[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] 20080103.01.tpd.patch (7090/355)
|
From: |
daly |
|
Subject: |
[Axiom-developer] 20080103.01.tpd.patch (7090/355) |
|
Date: |
Fri, 4 Jan 2008 21:51:11 -0600 |
Handle besselK differently. Apply Waldek's patches. Fix 7090/355.
(1) -> D(besselK(a,x),x)
- besselK(a + 1,x) - besselK(a - 1,x)
(1) -------------------------------------
2
Type: Expression Integer
(2) -> D(besselK(a,x),a)
(2) besselK (a,x)
,1
Type: Expression Integer
(3) -> integrate(D(besselK(a,x),a),a)
(3) besselK(a,x)
Type: Union(Expression Integer,...)
(4) -> limit(D(besselK(a,x),a),a=1/2)
(4) "failed"
Type: Union("failed",...)
=======================================================================
diff --git a/changelog b/changelog
index daecc5d..8acaf43 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,6 @@
+20080103 wxh src/algebra/sf.spad handle besselK (7090/355)
+20080103 wxh src/algebra/op.spad handle besselK (7090/355)
+20080103 wxh src/algebra/combfunc.spad handle besselK (7090/355)
20080102 tpd src/hyper/Makefile.pamphlet fix typo for axbook
20071230 tpd src/hyper/Makefile prevent spurious remake of axbook (7052)
20071230 tpd src/input/summation.input update tests with new mathml output
diff --git a/src/algebra/combfunc.spad.pamphlet
b/src/algebra/combfunc.spad.pamphlet
index c7d54b7..9da05e9 100644
--- a/src/algebra/combfunc.spad.pamphlet
+++ b/src/algebra/combfunc.spad.pamphlet
@@ -83,6 +83,53 @@ CombinatorialFunction(R, F): Exports == Implementation where
binomial : (F, F) -> F
++ binomial(n, r) returns the number of subsets of r objects
++ taken among n objects, i.e. n!/(r! * (n-r)!);
+@
+
+We currently simplify binomial coefficients only for non-negative integral
+second argument, using the formula
+$$ \binom{n}{k}=\frac{1}{k!}\prod_{i=0..k-1} (n-i),$$
+except if the second argument is symbolic: in this case [[binomial(n,n)]] is
+simplified to one.
+
+Note that there are at least two different ways to define binomial coefficients
+for negative integral second argument. One way, particular suitable for
+combinatorics, is to set the binomial coefficient equal to zero for negative
+second argument. This is, partially, also the approach taken in
+[[combinat.spad]], where we find
+
+\begin{verbatim}
+ binomial(n, m) ==
+ n < 0 or m < 0 or m > n => 0
+ m = 0 => 1
+\end{verbatim}
+
+Of course, here [[n]] and [[m]] are integers. This definition agrees with the
+recurrence
+
+$$\binom{n}{k}+\binom{n}{k+1}=\binom{n+1}{k+1}.$$
+
+Alternatively, one can use the formula
+$$ \binom{n}{k}=\frac{\Gamma(n+1)}{\Gamma(k+1)\Gamma(n-k+1)}, $$
+and leave the case where $k\in\mathbb Z$, $n\in\mathbb Z$ and $k \leq n < 0$
+undefined, since the limit does not exist in this case:
+
+Since we then have that $n-k+1\geq 1$, $\Gamma(n-k+1)$ is finite. So it is
+sufficient to consider $\frac{\Gamma(n+1)}{\Gamma(k+1)}$. On the one hand, we
+have
+$$\lim_{n_0\to n} \lim_{k_0\to k}\frac{\Gamma(n_0+1)}{\Gamma(k_0+1)} = 0,$$
+since for any non-integral $n_0$, $\Gamma(n_0+1)$ is finite. On the other
+hand,
+$$\lim_{k_0\to k} \lim_{n_0\to n}\frac{\Gamma(n_0+1)}{\Gamma(k_0+1)}$$
+does not exist, since for non-integral $k_0$, $\Gamma(k_0+1)$ is finite while
+$\Gamma(n_0+1)$ is unbounded.
+
+However, since for $k\in\mathbb Z$, $n\in\mathbb Z$ and $0 < k < n$ both
+definitions agree, one could also combine them. This is what, for example,
+Mathematica does. It seems that MuPAD sets [[binomial(n,n)=1]] for all
+arguments [[n]], and returns [[binomial(-2, n)]] unevaluated. Provisos may help
+here.
+
+<<package COMBF CombinatorialFunction>>=
permutation: (F, F) -> F
++ permutation(n, r) returns the number of permutations of
++ n objects taken r at a time, i.e. n!/(n-r)!;
@@ -517,17 +564,33 @@ dummy variable is introduced to make the indexing
variable \lq local\rq.
if R has RetractableTo(Z) and F has Algebra(Fraction(Z)) then
iibinom l ==
+ (s:=retractIfCan(second l)@Union(R,"failed")) case R and
+ (t:=retractIfCan(s)@Union(Z,"failed")) case Z and t>0 =>
+ ans:=1::F
+ for i in 0..t-1 repeat
+ ans:=ans*(first l - i::R::F)
+ (1/factorial t) * ans
(s:=retractIfCan(first l-second l)@Union(R,"failed")) case R and
- (t:=retractIfCan(s)@Union(Z,"failed")) case Z and s>0=>
- ans:=1::F
- for i in 1..t repeat
- ans:=ans*(second l+i::R::F)
- (1/factorial t) * ans
+ (t:=retractIfCan(s)@Union(Z,"failed")) case Z and t>0 =>
+ ans:=1::F
+ for i in 1..t repeat
+ ans:=ans*(second l+i::R::F)
+ (1/factorial t) * ans
(r1 := retractIfCan(first l)@Union(R,"failed")) case "failed" or
(r2 := retractIfCan(second l)@Union(R,"failed")) case "failed"
=> ibinom l
binomial(r1::R, r2::R)::F
+@
+
+[[iibinom]] checks those cases in which the binomial coefficient may be
+evaluated explicitly. Note that up to [[patch--51]], the case where the second
+argument is a positive integer was not checked.(Issue~\#336) Currently, the
+naive iterative algorithm is used to calculate the coefficient, there is room
+for improvement here.
+
+<<package COMBF CombinatorialFunction>>=
+
else
iibinom l ==
(r1 := retractIfCan(first l)@Union(R,"failed")) case "failed" or
@@ -622,6 +685,7 @@ FunctionalSpecialFunction(R, F): Exports == Implementation
where
OP ==> BasicOperator
K ==> Kernel F
SE ==> Symbol
+ SPECIALDIFF ==> "%specialDiff"
Exports ==> with
belong? : OP -> Boolean
@@ -635,46 +699,73 @@ FunctionalSpecialFunction(R, F): Exports ==
Implementation where
Gamma : F -> F
++ Gamma(f) returns the formal Gamma function applied to f
Gamma : (F,F) -> F
- ++ Gamma(a,x) returns the incomplete Gamma function applied to a and x.
- ++ Concerning differentiation, it is regarded as a function in the
- ++ second argument only.
+ ++ Gamma(a,x) returns the incomplete Gamma function applied to a and x
Beta: (F,F) -> F
++ Beta(x,y) returns the beta function applied to x and y
digamma: F->F
++ digamma(x) returns the digamma function applied to x
polygamma: (F,F) ->F
- ++ polygamma(x,y) returns the polygamma function applied to x and y.
- ++ Concerning differentiation, it is regarded as a function in the
- ++ second argument only.
+ ++ polygamma(x,y) returns the polygamma function applied to x and y
besselJ: (F,F) -> F
- ++ besselJ(x,y) returns the besselj function applied to x and y.
- ++ Concerning differentiation, it is regarded as a function in the
- ++ second argument only.
+ ++ besselJ(x,y) returns the besselj function applied to x and y
besselY: (F,F) -> F
- ++ besselY(x,y) returns the bessely function applied to x and y.
- ++ Concerning differentiation, it is regarded as a function in the
- ++ second argument only.
+ ++ besselY(x,y) returns the bessely function applied to x and y
besselI: (F,F) -> F
- ++ besselI(x,y) returns the besseli function applied to x and y.
- ++ Concerning differentiation, it is regarded as a function in the
- ++ second argument only.
+ ++ besselI(x,y) returns the besseli function applied to x and y
besselK: (F,F) -> F
- ++ besselK(x,y) returns the besselk function applied to x and y.
- ++ Concerning differentiation, it is regarded as a function in the
- ++ second argument only.
+ ++ besselK(x,y) returns the besselk function applied to x and y
airyAi: F -> F
++ airyAi(x) returns the airyai function applied to x
airyBi: F -> F
++ airyBi(x) returns the airybi function applied to x
+@
+
+In case we want to have more special function operators here, do not forget to
+add them to the list [[specop]] in [[CommonOperators]]. Otherwise they will
+not have the \lq special\rq\ attribute and will not be recognised here. One
+effect could be that
+\begin{verbatim}
+myNewSpecOp(1::Expression Integer)::Expression DoubleFloat
+\end{verbatim}
+might not re-evaluate the operator.
+
+<<package FSPECF FunctionalSpecialFunction>>=
iiGamma : F -> F
++ iiGamma(x) should be local but conditional;
iiabs : F -> F
++ iiabs(x) should be local but conditional;
+ iiBeta : List F -> F
+ ++ iiGamma(x) should be local but conditional;
+ iidigamma : F -> F
+ ++ iidigamma(x) should be local but conditional;
+ iipolygamma: List F -> F
+ ++ iipolygamma(x) should be local but conditional;
+ iiBesselJ : List F -> F
+ ++ iiBesselJ(x) should be local but conditional;
+ iiBesselY : List F -> F
+ ++ iiBesselY(x) should be local but conditional;
+ iiBesselI : List F -> F
+ ++ iiBesselI(x) should be local but conditional;
+ iiBesselK : List F -> F
+ ++ iiBesselK(x) should be local but conditional;
+ iiAiryAi : F -> F
+ ++ iiAiryAi(x) should be local but conditional;
+ iiAiryBi : F -> F
+ ++ iiAiryBi(x) should be local but conditional;
Implementation ==> add
- iabs : F -> F
- iGamma: F -> F
+ iabs : F -> F
+ iGamma : F -> F
+ iBeta : (F, F) -> F
+ idigamma : F -> F
+ iiipolygamma: (F, F) -> F
+ iiiBesselJ : (F, F) -> F
+ iiiBesselY : (F, F) -> F
+ iiiBesselI : (F, F) -> F
+ iiiBesselK : (F, F) -> F
+ iAiryAi : F -> F
+ iAiryBi : F -> F
opabs := operator("abs"::Symbol)$CommonOperators
opGamma := operator("Gamma"::Symbol)$CommonOperators
@@ -732,6 +823,17 @@ FunctionalSpecialFunction(R, F): Exports == Implementation
where
x < 0 => kernel(opabs, -x)
kernel(opabs, x)
+ iBeta(x, y) == kernel(opBeta, [x, y])
+ idigamma x == kernel(opdigamma, x)
+ iiipolygamma(n, x) == kernel(oppolygamma, [n, x])
+ iiiBesselJ(x, y) == kernel(opBesselJ, [x, y])
+ iiiBesselY(x, y) == kernel(opBesselY, [x, y])
+ iiiBesselI(x, y) == kernel(opBesselI, [x, y])
+ iiiBesselK(x, y) == kernel(opBesselK, [x, y])
+ iAiryAi x == kernel(opAiryAi, x)
+ iAiryBi x == kernel(opAiryBi, x)
+
+
-- Could put more conditional special rules for other functions here
if R has abs : R -> R then
@@ -750,6 +852,54 @@ FunctionalSpecialFunction(R, F): Exports == Implementation
where
(r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iGamma x
Gamma(r::R)::F
+ iiBeta l ==
+ (r:=retractIfCan(first l)@Union(R,"failed")) case "failed" or _
+ (s:=retractIfCan(second l)@Union(R,"failed")) case "failed" _
+ => iBeta(first l, second l)
+ Beta(r::R, s::R)::F
+
+ iidigamma x ==
+ (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => idigamma x
+ digamma(r::R)::F
+
+ iipolygamma l ==
+ (s:=retractIfCan(first l)@Union(R,"failed")) case "failed" or _
+ (r:=retractIfCan(second l)@Union(R,"failed")) case "failed" _
+ => iiipolygamma(first l, second l)
+ polygamma(s::R, r::R)::F
+
+ iiBesselJ l ==
+ (r:=retractIfCan(first l)@Union(R,"failed")) case "failed" or _
+ (s:=retractIfCan(second l)@Union(R,"failed")) case "failed" _
+ => iiiBesselJ(first l, second l)
+ besselJ(r::R, s::R)::F
+
+ iiBesselY l ==
+ (r:=retractIfCan(first l)@Union(R,"failed")) case "failed" or _
+ (s:=retractIfCan(second l)@Union(R,"failed")) case "failed" _
+ => iiiBesselY(first l, second l)
+ besselY(r::R, s::R)::F
+
+ iiBesselI l ==
+ (r:=retractIfCan(first l)@Union(R,"failed")) case "failed" or _
+ (s:=retractIfCan(second l)@Union(R,"failed")) case "failed" _
+ => iiiBesselI(first l, second l)
+ besselI(r::R, s::R)::F
+
+ iiBesselK l ==
+ (r:=retractIfCan(first l)@Union(R,"failed")) case "failed" or _
+ (s:=retractIfCan(second l)@Union(R,"failed")) case "failed" _
+ => iiiBesselK(first l, second l)
+ besselK(r::R, s::R)::F
+
+ iiAiryAi x ==
+ (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iAiryAi x
+ airyAi(r::R)::F
+
+ iiAiryBi x ==
+ (r:=retractIfCan(x)@Union(R,"failed")) case "failed" => iAiryBi x
+ airyBi(r::R)::F
+
else
if R has RetractableTo Integer then
iiGamma x ==
@@ -759,39 +909,113 @@ FunctionalSpecialFunction(R, F): Exports ==
Implementation where
else
iiGamma x == iGamma x
+ iiBeta l == iBeta(first l, second l)
+ iidigamma x == idigamma x
+ iipolygamma l == iiipolygamma(first l, second l)
+ iiBesselJ l == iiiBesselJ(first l, second l)
+ iiBesselY l == iiiBesselY(first l, second l)
+ iiBesselI l == iiiBesselI(first l, second l)
+ iiBesselK l == iiiBesselK(first l, second l)
+ iiAiryAi x == iAiryAi x
+ iiAiryBi x == iAiryBi x
+
-- Default behaviour is to build a kernel
evaluate(opGamma, iiGamma)$BasicOperatorFunctions1(F)
evaluate(opabs, iiabs)$BasicOperatorFunctions1(F)
+-- evaluate(opGamma2 ,iiGamma2 )$BasicOperatorFunctions1(F)
+ evaluate(opBeta ,iiBeta )$BasicOperatorFunctions1(F)
+ evaluate(opdigamma ,iidigamma )$BasicOperatorFunctions1(F)
+ evaluate(oppolygamma ,iipolygamma)$BasicOperatorFunctions1(F)
+ evaluate(opBesselJ ,iiBesselJ )$BasicOperatorFunctions1(F)
+ evaluate(opBesselY ,iiBesselY )$BasicOperatorFunctions1(F)
+ evaluate(opBesselI ,iiBesselI )$BasicOperatorFunctions1(F)
+ evaluate(opBesselK ,iiBesselK )$BasicOperatorFunctions1(F)
+ evaluate(opAiryAi ,iiAiryAi )$BasicOperatorFunctions1(F)
+ evaluate(opAiryBi ,iiAiryBi )$BasicOperatorFunctions1(F)
+@
+
+\subsection{differentiation of special functions}
+
+In the following we define the symbolic derivatives of the special functions we
+provide. The formulas we use for the Bessel functions can be found in Milton
+Abramowitz and Irene A. Stegun, eds. (1965). Handbook of Mathematical
+Functions with Formulas, Graphs, and Mathematical Tables. New York: Dover. ISBN
+0-486-61272-4, Equations~9.1.27 and 9.6.26. Up to [[patch--50]] the formula
+for $K$ missed the minus sign. (Issue~\#355)
+
+We do not attempt to provide formulas for the derivative with respect to the
+first argument currently. Instead, we leave such derivatives unevaluated.
+<<package FSPECF FunctionalSpecialFunction>>=
import Fraction Integer
ahalf: F := recip(2::F)::F
athird: F := recip(2::F)::F
twothirds: F := 2*recip(3::F)::F
+@
+
+We need to get hold of the differentiation operator as modified by
+[[FunctionSpace]]. Otherwise, for example, display will be ugly. We accomplish
+that by differentiating an operator, which will certainly result in a single
+kernel only.
+
+<<package FSPECF FunctionalSpecialFunction>>=
+ dummyArg: SE := new()$SE
+ opdiff := operator first kernels D((operator(new()$SE)$BasicOperator)
+ (dummyArg::F), dummyArg)
+@
+
+The differentiation operator [[opdiff]] takes three arguments corresponding to
+$$
+F_{,i}(a_1,a_2,\dots,a_n):
+$$
+\begin{enumerate}
+\item $F(a_1,...,dm,...a_n)$, where the $i$\textsuperscript{th} argument is a
+ dummy variable,
+\item $dm$, the dummy variable, and
+\item $a_i$, the point at which the differential is evaluated.
+\end{enumerate}
+
+In the following, it seems to be safe to use the same dummy variable
+troughout. At least, this is done also in [[FunctionSpace]], and did not cause
+problems.
- lzero(l: List F): F == 0
+The operation [[symbolicGrad]] returns the first component of the gradient of
+[[op l]].
+
+<<package FSPECF FunctionalSpecialFunction>>=
+ dm := new()$SE :: F
- iBesselJGrad(l: List F): F ==
+ iBesselJ(l: List F, t: SE): F ==
n := first l; x := second l
- ahalf * (besselJ (n-1,x) - besselJ (n+1,x))
- iBesselYGrad(l: List F): F ==
+ differentiate(n, t)*kernel(opdiff, [opBesselJ [dm, x], dm, n])
+ + differentiate(x, t) * ahalf * (besselJ (n-1,x) - besselJ (n+1,x))
+
+ iBesselY(l: List F, t: SE): F ==
n := first l; x := second l
- ahalf * (besselY (n-1,x) - besselY (n+1,x))
- iBesselIGrad(l: List F): F ==
+ differentiate(n, t)*kernel(opdiff, [opBesselY [dm, x], dm, n])
+ + differentiate(x, t) * ahalf * (besselY (n-1,x) - besselY (n+1,x))
+
+ iBesselI(l: List F, t: SE): F ==
n := first l; x := second l
- ahalf * (besselI (n-1,x) + besselI (n+1,x))
- iBesselKGrad(l: List F): F ==
+ differentiate(n, t)*kernel(opdiff, [opBesselI [dm, x], dm, n])
+ + differentiate(x, t)* ahalf * (besselI (n-1,x) + besselI (n+1,x))
+
+ iBesselK(l: List F, t: SE): F ==
n := first l; x := second l
- - ahalf * (besselK (n-1,x) + besselK (n+1,x))
+ differentiate(n, t)*kernel(opdiff, [opBesselK [dm, x], dm, n])
+ - differentiate(x, t)* ahalf * (besselK (n-1,x) + besselK (n+1,x))
+
@
-The formulas above for the Bessel functions can be found in Milton Abramowitz
-and Irene A. Stegun, eds. (1965). Handbook of Mathematical Functions with
-Formulas, Graphs, and Mathematical Tables. New York: Dover. ISBN 0-486-61272-4,
-Equations~9.1.27 and 9.6.26. Up to [[patch--50]] the formula for $K$ missed
-the minus sign. (Issue~\#355)
+
+For the moment we throw an error if we try to differentiate [[polygamma]] with
+respect to the first argument.
+
<<package FSPECF FunctionalSpecialFunction>>=
- ipolygammaGrad(l: List F): F ==
- n := first l; x := second l
- polygamma(n+1, x)
+ ipolygamma(l: List F, x: SE): F ==
+ member?(x, variables first l) =>
+ error "cannot differentiate polygamma with respect to the first
argument"
+ n := first l; y := second l
+ differentiate(y, x)*polygamma(n+1, y)
iBetaGrad1(l: List F): F ==
x := first l; y := second l
Beta(x,y)*(digamma x - digamma(x+y))
@@ -800,20 +1024,36 @@ the minus sign. (Issue~\#355)
Beta(x,y)*(digamma y - digamma(x+y))
if F has ElementaryFunctionCategory then
- iGamma2Grad(l: List F):F ==
+ iGamma2(l: List F, t: SE): F ==
a := first l; x := second l
- - x ** (a - 1) * exp(-x)
- derivative(opGamma2, [lzero, iGamma2Grad])
+ differentiate(a, t)*kernel(opdiff, [opGamma2 [dm, x], dm, a])
+ - differentiate(x, t)* x ** (a - 1) * exp(-x)
+ setProperty(opGamma2, SPECIALDIFF, iGamma2@((List F, SE)->F)
+ pretend None)
+@
+
+Finally, we tell Axiom to use these functions for differentiation. Note that
+up to [[patch--50]], the properties for the Bessel functions were set using
+[[derivative(oppolygamma, [lzero, ipolygammaGrad])]], where [[lzero]] returned
+zero always. Trying to replace [[lzero]] by a function that returns the first
+component of the gradient failed, it resulted in an infinite loop for
+[[integrate(D(besselJ(a,x),a),a)]].
+<<package FSPECF FunctionalSpecialFunction>>=
derivative(opabs, abs(#1) * inv(#1))
derivative(opGamma, digamma #1 * Gamma #1)
derivative(opBeta, [iBetaGrad1, iBetaGrad2])
derivative(opdigamma, polygamma(1, #1))
- derivative(oppolygamma, [lzero, ipolygammaGrad])
- derivative(opBesselJ, [lzero, iBesselJGrad])
- derivative(opBesselY, [lzero, iBesselYGrad])
- derivative(opBesselI, [lzero, iBesselIGrad])
- derivative(opBesselK, [lzero, iBesselKGrad])
+ setProperty(oppolygamma, SPECIALDIFF, ipolygamma@((List F, SE)->F)
+ pretend None)
+ setProperty(opBesselJ, SPECIALDIFF, iBesselJ@((List F, SE)->F)
+ pretend None)
+ setProperty(opBesselY, SPECIALDIFF, iBesselY@((List F, SE)->F)
+ pretend None)
+ setProperty(opBesselI, SPECIALDIFF, iBesselI@((List F, SE)->F)
+ pretend None)
+ setProperty(opBesselK, SPECIALDIFF, iBesselK@((List F, SE)->F)
+ pretend None)
@
\section{package SUMFS FunctionSpaceSum}
diff --git a/src/algebra/op.spad.pamphlet b/src/algebra/op.spad.pamphlet
index e443d3c..51d4d6c 100644
--- a/src/algebra/op.spad.pamphlet
+++ b/src/algebra/op.spad.pamphlet
@@ -460,9 +460,11 @@ BasicOperator(): Exports == Implementation where
-- property EQUAL? contains a function f: (BOP, BOP) -> Boolean
-- such that f(o1, o2) is true iff o1 = o2
op1 = op2 ==
+ (EQ$Lisp)(op1, op2) => true
name(op1) ^= name(op2) => false
op1.narg ^= op2.narg => false
- brace(keys properties op1)^=$Set(String) brace(keys properties op2) =>
false
+ brace(keys properties op1)^=$Set(String) _
+ brace(keys properties op2) => false
(func := property(op1, EQUAL?)) case None =>
((func::None) pretend (($, $) -> Boolean)) (op1, op2)
true
@@ -721,7 +723,8 @@ CommonOperators(): Exports == Implementation where
combop := [opfact, opperm, opbinom, oppow,
opsum, opdsum, opprod, opdprod]
specop := [opGamma, opGamma2, opBeta, opdigamma, oppolygamma, opabs,
- opBesselJ, opBesselY, opBesselI, opBesselK]
+ opBesselJ, opBesselY, opBesselI, opBesselK, opAiryAi,
+ opAiryBi]
anyop := [oppren, opdiff, opbox, opquote]
allop := concat(concat(concat(concat(concat(
algop,elemop),primop),combop),specop),anyop)
@@ -746,7 +749,23 @@ CommonOperators(): Exports == Implementation where
dpi l == "%pi"::Symbol::O
dfact x == postfix("!"::Symbol::O, (ATOM(x)$Lisp => x; paren x))
dquote l == prefix(quote(first(l)::O), rest l)
+@
+It is certainly an abuse of OutputForm to produce a Gamma as done below.
+Originally, it was even worse (Issue~\#6):
+\begin{verbatim}
dgamma l == prefix(hconcat("|"::Symbol::O, overbar(" "::Symbol::O)), l)
+\end{verbatim}
+which was TeXed to
+$${|{\overline{\ }}}
+\left(
+{x}
+\right).
+$$
+
+The right thing would be to introduce Greek letters in OutputForm, but
+that should be coordinated with the new mathml package
+<<package COMMONOP CommonOperators>>=
+ dgamma l == prefix(super("|"::Symbol::O, "-"::Symbol::O), l)
setDummyVar(op, n) == setProperty(op, DUMMYVAR, n pretend None)
dexp x ==
diff --git a/src/algebra/sf.spad.pamphlet b/src/algebra/sf.spad.pamphlet
index 8ed563d..571ccd0 100644
--- a/src/algebra/sf.spad.pamphlet
+++ b/src/algebra/sf.spad.pamphlet
@@ -993,7 +993,8 @@ o $AXIOM/doc/src/algebra/sf.spad.dvi
-- I've put some timing comparisons in the notes for the Float
-- domain about the difference in speed between the two domains.
DoubleFloat(): Join(FloatingPointSystem, DifferentialRing, OpenMath,
- TranscendentalFunctionCategory, ConvertibleTo InputForm) with
+ TranscendentalFunctionCategory, SpecialFunctionCategory, _
+ ConvertibleTo InputForm) with
_/ : (%, Integer) -> %
++ x / i computes the division from x by an integer i.
_*_* : (%,%) -> %
@@ -1082,16 +1083,18 @@ DoubleFloat(): Join(FloatingPointSystem,
DifferentialRing, OpenMath,
base() = 2 => precision()
base() = 16 => 4*precision()
wholePart(precision()*log2(base()::%))::PositiveInteger
- max() == MOST_-POSITIVE_-LONG_-FLOAT$Lisp
- min() == MOST_-NEGATIVE_-LONG_-FLOAT$Lisp
+ max() == MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp
+ min() == MOST_-NEGATIVE_-DOUBLE_-FLOAT$Lisp
order(a) == precision() + exponent a - 1
- 0 == FLOAT(0$Lisp,MOST_-POSITIVE_-LONG_-FLOAT$Lisp)$Lisp
- 1 == FLOAT(1$Lisp,MOST_-POSITIVE_-LONG_-FLOAT$Lisp)$Lisp
+ 0 == FLOAT(0$Lisp,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp
+ 1 == FLOAT(1$Lisp,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp
-- rational approximation to e accurate to 23 digits
- exp1() ==
FLOAT(534625820200,MOST_-POSITIVE_-LONG_-FLOAT$Lisp)$Lisp /
FLOAT(196677847971,MOST_-POSITIVE_-LONG_-FLOAT$Lisp)$Lisp
- pi() == PI$Lisp
+ exp1() == FLOAT(534625820200,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp / _
+ FLOAT(196677847971,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp
+ pi() == FLOAT(PI$Lisp,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp
coerce(x:%):OutputForm ==
- outputForm(FORMAT(NIL$Lisp,format,x)$Lisp pretend DoubleFloat)
+ x >= 0 => message(FORMAT(NIL$Lisp,format,x)$Lisp pretend String)
+ - (message(FORMAT(NIL$Lisp,format,-x)$Lisp pretend String))
convert(x:%):InputForm == convert(x pretend DoubleFloat)$InputForm
x < y == (x<y)$Lisp
- x == (-x)$Lisp
@@ -1107,7 +1110,7 @@ DoubleFloat(): Join(FloatingPointSystem,
DifferentialRing, OpenMath,
log10 x == checkComplex log(x)$Lisp
x:% ** i:Integer == EXPT(x,i)$Lisp
x:% ** y:% == checkComplex EXPT(x,y)$Lisp
- coerce(i:Integer):% == FLOAT(i,MOST_-POSITIVE_-LONG_-FLOAT$Lisp)$Lisp
+ coerce(i:Integer):% == FLOAT(i,MOST_-POSITIVE_-DOUBLE_-FLOAT$Lisp)$Lisp
exp x == EXP(x)$Lisp
log x == checkComplex LN(x)$Lisp
log2 x == checkComplex LOG2(x)$Lisp
@@ -1145,8 +1148,22 @@ DoubleFloat(): Join(FloatingPointSystem,
DifferentialRing, OpenMath,
SFSFUN ==> DoubleFloatSpecialFunctions()
sfx ==> x pretend DoubleFloat
sfy ==> y pretend DoubleFloat
- Gamma x == Gamma(sfx)$SFSFUN pretend %
+ airyAi x == airyAi(sfx)$SFSFUN pretend %
+ airyBi x == airyBi(sfx)$SFSFUN pretend %
+ besselI(x,y) == besselI(sfx,sfy)$SFSFUN pretend %
+ besselJ(x,y) == besselJ(sfx,sfy)$SFSFUN pretend %
+ besselK(x,y) == besselK(sfx,sfy)$SFSFUN pretend %
+ besselY(x,y) == besselY(sfx,sfy)$SFSFUN pretend %
Beta(x,y) == Beta(sfx,sfy)$SFSFUN pretend %
+ digamma x == digamma(sfx)$SFSFUN pretend %
+ Gamma x == Gamma(sfx)$SFSFUN pretend %
+-- not implemented in SFSFUN
+-- Gamma(x,y) == Gamma(sfx,sfy)$SFSFUN pretend %
+ polygamma(x,y) ==
+ if (n := retractIfCan(x:%):Union(Integer, "failed")) case Integer _
+ and n >= 0
+ then polygamma(n::Integer::NonNegativeInteger,sfy)$SFSFUN pretend %
+ else error "polygamma: first argument should be a nonnegative integer"
wholePart x == FIX(x)$Lisp
float(ma,ex,b) == ma*(b::%)**ex
- [Axiom-developer] 20080103.01.tpd.patch (7090/355),
daly <=