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. + +<>= 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. + +<>= + 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. + +<>= 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. +<>= 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. + +<>= + 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]]. + +<>= + 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. + <>= - 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)]]. +<>= 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 +<>= + 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 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