Dear Bill, Tim, Waldek,
here are some patches for bug #355.
Index: special.spad.pamphlet
===================================================================
--- special.spad.pamphlet (revision 580)
+++ special.spad.pamphlet (working copy)
@@ -86,7 +86,7 @@
++ \spad{Y(v,x)}.
++ This function satisfies the differential equation:
++ \spad{x^2 w''(x) + x w'(x) + (x^2-v^2)w(x) = 0}.
- ++ Note: The default implmentation uses the relation
+ ++ Note: The default implementation uses the relation
++ \spad{Y(v,x) = (J(v,x) cos(v*%pi) - J(-v,x))/sin(v*%pi)}
++ so is not valid for integer values of v.
besselY: (C, C) -> C
@@ -94,7 +94,7 @@
++ \spad{Y(v,x)}.
++ This function satisfies the differential equation:
++ \spad{x^2 w''(x) + x w'(x) + (x^2-v^2)w(x) = 0}.
- ++ Note: The default implmentation uses the relation
+ ++ Note: The default implementation uses the relation
++ \spad{Y(v,x) = (J(v,x) cos(v*%pi) - J(-v,x))/sin(v*%pi)}
++ so is not valid for integer values of v.
@@ -110,19 +110,19 @@
++ \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
besselK: (R, R) -> R
- ++ besselK(v,x) is the modified Bessel function of the first kind,
+ ++ besselK(v,x) is the modified Bessel function of the second kind,
++ \spad{K(v,x)}.
++ This function satisfies the differential equation:
++ \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
- ++ Note: The default implmentation uses the relation
+ ++ Note: The default implementation uses the relation
++ \spad{K(v,x) = %pi/2*(I(-v,x) - I(v,x))/sin(v*%pi)}.
++ so is not valid for integer values of v.
besselK: (C, C) -> C
- ++ besselK(v,x) is the modified Bessel function of the first kind,
+ ++ besselK(v,x) is the modified Bessel function of the second kind,
++ \spad{K(v,x)}.
++ This function satisfies the differential equation:
++ \spad{x^2 w''(x) + x w'(x) - (x^2+v^2)w(x) = 0}.
- ++ Note: The default implmentation uses the relation
+ ++ Note: The default implementation uses the relation
++ \spad{K(v,x) = %pi/2*(I(-v,x) - I(v,x))/sin(v*%pi)}
++ so is not valid for integer values of v.
Index: combfunc.spad.pamphlet
===================================================================
--- combfunc.spad.pamphlet (revision 580)
+++ combfunc.spad.pamphlet (working copy)
@@ -41,11 +41,6 @@
++ formal product;
@
-The latest change allows Axiom to reduce
-\begin{verbatim}
- sum(1/i,i=1..n)-sum(1/i,i=1..n)
-\end{verbatim}
-to reduce to zero.
<<package COMBF CombinatorialFunction>>=
)abbrev package COMBF CombinatorialFunction
++ Provides the usual combinatorial functions
@@ -703,21 +698,33 @@
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
+ ++ 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.
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
+ ++ 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.
besselJ: (F,F) -> F
- ++ besselJ(x,y) returns the besselj function applied to x and y
+ ++ 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.
besselY: (F,F) -> F
- ++ besselY(x,y) returns the bessely function applied to x and y
+ ++ 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.
besselI: (F,F) -> F
- ++ besselI(x,y) returns the besseli function applied to x and y
+ ++ 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.
besselK: (F,F) -> F
- ++ besselK(x,y) returns the besselk function applied to x and y
+ ++ 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.
airyAi: F -> F
++ airyAi(x) returns the airyai function applied to x
airyBi: F -> F
@@ -837,7 +844,16 @@
ahalf * (besselI (n-1,x) + besselI (n+1,x))
iBesselKGrad(l: List F): F ==
n := first l; x := second l
- ahalf * (besselK (n-1,x) + besselK (n+1,x))
+ - 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)
+
+<<package FSPECF FunctionalSpecialFunction>>=
ipolygammaGrad(l: List F): F ==
n := first l; x := second l
polygamma(n+1, x)
Bill: I have been unable to upload the files to MathAction again, and could not
find the instructions. Maybe you could provide some help on MathAction
itsself?
Tim, Waldek: would you be so kind to include these patches in the respective
distributions. By the way, I just found the documented version of the patch to
STTAYLOR, it is on MathAction, #312, powern.patch. I include a version that
applies smoothly to wh-sandbox here, too. Please apply this one also.
Index: sttaylor.spad.pamphlet
===================================================================
--- sttaylor.spad.pamphlet (revision 580)
+++ sttaylor.spad.pamphlet (working copy)
@@ -401,6 +401,14 @@
cssa(frst ststa,mapsa((rst x) * #1,comps(rst ststa,x)))
if A has Algebra RN then
+@
+
+The following defines lazy integration on streams interpreted as Taylor series.
+I.e. if [[x]] is $c_0,c_1,c_2,\dots$, then [[integ x]] returns
+$c_0,\frac{1}{2}c_1,\frac{1}{3}c_2,\dots$. [[integrate]] prepends a given
+constant of integration.
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
integre: (ST A,I) -> ST A
integre(x,n) == delay
empty? x => zro()
@@ -412,6 +420,9 @@
integrate(a,x) == concat(a,integ x)
lazyIntegrate(s,xf) == concat(s,integ(delay xf))
+@
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
nldere:(ST ST A,ST A) -> ST A
nldere(lslsa,c) == lazyIntegrate(0,addiag(comps(lslsa,c)))
nlde lslsa == YS(nldere(lslsa,#1))
@@ -420,10 +431,45 @@
smult: (RN,ST A) -> ST A
smult(rn,x) == map(rn * #1,x)
+@
+
+The following helper function computes
+\begin{equation*}
+ 1+\int (rn+1) c x' dz - c (x-a_0),
+\end{equation*}
+where $a_0$ is the constant term of [[x]].
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
powerrn:(RN,ST A,ST A) -> ST A
powerrn(rn,x,c) == delay
concat(1,integ(smult(rn + 1,c * deriv x)) - rst x * c)
- powern(rn,x) ==
+@
+
+The following operation raises the power series [[x]] to a rational power
+[[rn]]. We first outline the general strategy.
+
+Suppose the constant term of [[x]] equals one. In this case, we have
+\begin{equation*}
+ x^{rn+1} = 1 + \int (rn+1) x^{rn} x' dz,
+\end{equation*}
+since $(x^{rn+1})'= (rn+1)x^{rn} x'$. Thus, $x^{rn}$ is the fixed point of %
+[[g +-> powerrn(rn, x, g)]]. We use [[Y$ParadoxicalCombinatorsForStreams(A)]]%$
+to compute this fixed point lazily.
+
+If the constant term of [[x]] is neither zero nor one, we use the identity
+\begin{equation*}
+ (c_0 + c_1*z + c_2 z^2\dots)^{rn} = c_0^{rn} (1 + \frac{c_1}{c_0}*z +\dots)^{rn}.
+\end{equation*}
+
+Finally, if the constant term of [[x]] is zero, we use the identity
+\begin{equation*}
+ (c_0*z^o + c_1*z^{o+1} +\dots)^{rn} = z^{o rn} (c_0 + c_1*z +\dots)^{rn}.
+\end{equation*}
+
+This implementation should be compared with [[cRationalPower$ISUPS]].%$
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
+ powern(rn, x) ==
order : I := 0
for n in 0.. repeat
empty? x => return zro()
@@ -431,17 +477,46 @@
x := rst x
n = 1000 =>
error "**: series with many leading zero coefficients"
+@
+First we determine the order of [[x]], i.e., the index of the first non-zero
+coefficient.
+
+Remarks:
+\begin{itemize}
+\item Note that usually [[leave]] takes no arguments. I don't know what it does
+ here.
+\item [[empty? x]] tests whether the stream [[x]] has no elements. This is
+ mathematically the same as the corresponding power series being zero.
+\end{itemize}
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
(ord := (order exquo denom(rn))) case "failed" =>
error "**: rational power does not exist"
- co := frst x
+@
+
+[[ord*numer(rn)]] will be the order of the new power series. If it is not an
+integer, we won't get a Taylor expansion and thus raise an error.
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
if ord > 0 and rn < 0 then
error "**: negative power does not exist"
+@
+
+If [[order]] was non-zero, we may not raise to a negative power. This test
+should really be done before the previous one.
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
+ co := frst x
(invCo := recip co) case "failed" =>
error "** rational power of coefficient undefined"
--- This error message is misleading, isn't it? see sups.spad/cRationalPower
+@
+
+We need to be able to invert the leading coefficient. The error message is
+slightly misleading, see [[sups.spad/cRationalPower]].
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
power :=
--- one? co => YS(powerrn(rn,x,#1))
- (co = 1) => YS(powerrn(rn,x,#1))
+ one? co => YS(powerrn(rn,x,#1))
(denom rn) = 1 =>
not negative?(num := numer rn) =>
-- It seems that this cannot happen, but I don't know why
@@ -450,8 +525,23 @@
RATPOWERS => co**rn * YS(powerrn(rn,(invCo :: A) * x,#1))
error "** rational power of coefficient undefined"
+@
+
+We now invoke the fixed point computation as explained above. We can do the
+computation if
+\begin{itemize}
+\item the constant term equals one, or
+\item [[rn]] is an integer, or
+\item we have rational powers in the coefficient ring.
+\end{itemize}
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
monom(1,(ord :: I) * numer(rn)) * power
+@
+Finally, we return the result with the correct order.
+
+<<package STTAYLOR StreamTaylorSeriesOperations>>=
if A has Field then
mapdiv(x,y) == delay
empty? y => error "stream division by zero"
Waldek: do you prefer if I commit to your branch myself, or is this form good
for you, too?
Many thanks,
Martin
_______________________________________________
Axiom-developer mailing list
[email protected]
http://lists.nongnu.org/mailman/listinfo/axiom-developer
