Ondrej,
Most of this patch concerns support for 'conjugate' is FriCAS, but
specifically the following code implements your Wirtinger-based chain
rule for abs:
+ dvabs(arg, z) ==
+ f := first arg
+
(conjugate(f)*differentiate(f,z)+f*differentiate(conjugate(f),z))*inv(2*abs(f))
+
+ --derivative(opabs, (x : F) : F +-> conjugate(x)*inv(2*abs(x)))
+ setProperty(opabs, SPECIALDIFF, dvabs@((List F, SE)->F) pretend None)
On 19 November 2014 12:29, Bill Page <[email protected]> wrote:
> Since this mostly concerns FriCAS I am cross posting to that group. I will
> also post the patch there. For FriCAS list reference the original email
> thread is here:
>
> https://groups.google.com/forum/#!topic/sage-devel/6j-LcC6tpkE
> ...
>> >
>> > Would you mind posting your patch to FriCAS somewhere? I would be
>> > interested in how you implemented it.
>>
>> I'll try to compile FriCAS myself and apply your patch, so that I can
>> play with it. Can you also try:
>>
>> abs(I*x)
>>
>> 1/abs(x)
>> 1/abs(x)^2
>> x/abs(x)^3
>> abs(x)^2
>>
>> The x/abs(x)^3 is a Coulomb's law in 1D.
>>
>> Ondrej
>>
>
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.
M src/algebra/combfunc.spad
M src/algebra/expr.spad
M src/algebra/sf.spad
M src/algebra/trigcat.spad
Index: src/algebra/combfunc.spad
===================================================================
--- src/algebra/combfunc.spad (revision 1825)
+++ src/algebra/combfunc.spad (working copy)
@@ -458,15 +458,16 @@
setProperty(opdprod, SPECIALEQUAL, equaldsumprod@((K, K) -> Boolean) pretend None)
-
)abbrev package FSPECF FunctionalSpecialFunction
)boot $tryRecompileArguments := nil
++ Provides the special functions
++ Author: Manuel Bronstein
++ Date Created: 18 Apr 1989
-++ Date Last Updated: 4 October 1993
++ Description: Provides some special functions over an integral domain.
++ Keywords: special, function.
+++ Date Last Updated: 2 October 2014
+++ Added: 'conjugate'
+++ Modied: 'abs'
FunctionalSpecialFunction(R, F) : Exports == Implementation where
R : Join(Comparable, IntegralDomain)
F : FunctionSpace R
@@ -485,6 +486,8 @@
++ error if op is not a special function operator
abs : F -> F
++ abs(f) returns the absolute value operator applied to f
+ conjugate: F -> F
+ ++ conjugate(f) returns the conjugate value operator applied to f
Gamma : F -> F
++ Gamma(f) returns the formal Gamma function applied to f
Gamma : (F, F) -> F
@@ -639,6 +642,8 @@
++ iiGamma(x) should be local but conditional;
iiabs : F -> F
++ iiabs(x) should be local but conditional;
+ iiconjugate: F -> F
+ ++ iiconjugate(x) should be local but conditional;
iiBeta : List F -> F
++ iiBeta(x) should be local but conditional;
iidigamma : F -> F
@@ -685,6 +690,8 @@
SPECIALINPUT ==> '%specialInput
iabs : F -> F
+ dvabs : (List F, SE) -> F
+ iconjugate: F -> F
iGamma : F -> F
iBeta : (F, F) -> F
idigamma : F -> F
@@ -708,6 +715,7 @@
opabs := operator('abs)$CommonOperators
+ opconjugate := operator('conjugate)$CommonOperators
opGamma := operator('Gamma)$CommonOperators
opGamma2 := operator('Gamma2)$CommonOperators
opBeta := operator('Beta)$CommonOperators
@@ -740,6 +748,7 @@
op_erfis := operator('%erfis)$CommonOperators
abs x == opabs x
+ conjugate x == opconjugate x
Gamma(x) == opGamma(x)
Gamma(a, x) == opGamma2(a, x)
Beta(x, y) == opBeta(x, y)
@@ -1685,6 +1694,7 @@
operator op ==
is?(op, 'abs) => opabs
+ is?(op, 'conjugate)=> opconjugate
is?(op, 'Gamma) => opGamma
is?(op, 'Gamma2) => opGamma2
is?(op, 'Beta) => opBeta
@@ -1757,10 +1767,18 @@
iabs x ==
zero? x => 0
+ one? x => 1
is?(x, opabs) => x
+ is?(x, opconjugate) => kernel(opabs, argument(retract(x)@K)(1))
smaller?(x, 0) => kernel(opabs, -x)
kernel(opabs, x)
+ iconjugate x ==
+ zero? x => 0
+ is?(x, opconjugate) => argument(retract(x)@K)(1)
+ is?(x, opabs) => x
+ kernel(opconjugate, x)
+
iBeta(x, y) == kernel(opBeta, [x, y])
idigamma x == kernel(opdigamma, x)
iiipolygamma(n, x) == kernel(oppolygamma, [n, x])
@@ -1836,8 +1854,26 @@
(a := retractIfCan(numer f)@Union(R, "failed")) case "failed" or
(b := retractIfCan(denom f)@Union(R,"failed")) case "failed" => iabs x
abs(a::R)::F / abs(b::R)::F
+ else
+ if F has RadicalCategory and R has conjugate : R -> R then
+ iiabs x ==
+ (r := retractIfCan(x)@Union(R, "failed"))
+ case "failed" => iabs x
+ sqrt( (r::R*conjugate(r::R))::F)
+ else iiabs x == iabs x
- else iiabs x == iabs x
+ iiconjugate(x:F):F ==
+ is?(x, opconjugate) => argument(retract(x)@K)(1)
+ is?(x, opabs) => x
+ retractIfCan(x)@Union(Symbol, "failed") case Symbol => iconjugate(x)
+ x:=eval(x,kernels(x), _
+ map((k:Kernel F):F +-> _
+ ( (height(k)=0 or height(k)=1 and retractIfCan(k::F)@Union(Symbol, "failed") case Symbol) =>iconjugate(k::F);map(iiconjugate,k)), _
+ kernels(x))$ListFunctions2(Kernel F,F))
+ if R has conjugate : R -> R then
+ x:=map(conjugate$R,numer x)::F / _
+ map(conjugate$R,denom x)::F
+ return x
if R has SpecialFunctionCategory then
iiGamma x ==
@@ -1929,6 +1965,7 @@
-- Default behaviour is to build a kernel
evaluate(opGamma, iiGamma)$BasicOperatorFunctions1(F)
evaluate(opabs, iiabs)$BasicOperatorFunctions1(F)
+ evaluate(opconjugate, iiconjugate)$BasicOperatorFunctions1(F)
-- evaluate(opGamma2 , iiGamma2 )$BasicOperatorFunctions1(F)
evaluate(opBeta , iiBeta )$BasicOperatorFunctions1(F)
evaluate(opdigamma , iidigamma )$BasicOperatorFunctions1(F)
@@ -2119,6 +2156,9 @@
SEX ==> SExpression
NNI ==> NonNegativeInteger
+ dconjugate(lo : List OF) : OF == overbar lo.1
+ display(opconjugate,dconjugate)
+
if F has RetractableTo(Integer) then
get_int_listf : List F -> List Integer
get_int_listo : (Integer, List OF) -> List Integer
@@ -2369,7 +2409,14 @@
derivative(op_erfs, d_erfs)
derivative(op_erfis, d_erfis)
- derivative(opabs, (x : F) : F +-> abs(x)*inv(x))
+ dvabs(arg, z) ==
+ f := first arg
+ (conjugate(f)*differentiate(f,z)+f*differentiate(conjugate(f),z))*inv(2*abs(f))
+
+ --derivative(opabs, (x : F) : F +-> conjugate(x)*inv(2*abs(x)))
+ setProperty(opabs, SPECIALDIFF, dvabs@((List F, SE)->F) pretend None)
+
+ derivative(opconjugate, (x : F) : F +-> 1)
derivative(opGamma, (x : F) : F +-> digamma(x)*Gamma(x))
derivative(op_log_gamma, (x : F) : F +-> digamma(x))
derivative(opBeta, [iBetaGrad1, iBetaGrad2])
@@ -2402,6 +2449,7 @@
setProperty(opPolylog, SPECIALDIFF, dPolylog@((List F, SE)->F)
pretend None)
+
)abbrev package SUMFS FunctionSpaceSum
++ Top-level sum function
++ Author: Manuel Bronstein
Index: src/algebra/expr.spad
===================================================================
--- src/algebra/expr.spad (revision 1825)
+++ src/algebra/expr.spad (working copy)
@@ -302,6 +302,7 @@
acsch x == acsch(x)$EF
abs x == abs(x)$FSF
+ conjugate x == conjugate(x)$FSF
Gamma x == Gamma(x)$FSF
Gamma(a, x) == Gamma(a, x)$FSF
Beta(x, y) == Beta(x, y)$FSF
Index: src/algebra/sf.spad
===================================================================
--- src/algebra/sf.spad (revision 1825)
+++ src/algebra/sf.spad (working copy)
@@ -572,9 +572,10 @@
sign(x) == retract FLOAT_-SIGN(x, 1)$Lisp
abs x == abs_DF(x)$Lisp
+ -- floats are real
+ conjugate x == x
-
manexp(x) ==
zero? x => [0, 0]
s := sign x; x := abs x
Index: src/algebra/trigcat.spad
===================================================================
--- src/algebra/trigcat.spad (revision 1825)
+++ src/algebra/trigcat.spad (working copy)
@@ -261,6 +261,8 @@
SpecialFunctionCategory() : Category == with
abs : % -> %
++ abs(x) returns the absolute value of x.
+ conjugate : % -> %
+ ++ conjugate(x) returns the conjugate of x.
Gamma : % -> %
++ Gamma(x) is the Euler Gamma function.
Beta : (%, %)->%
