aboutsummaryrefslogtreecommitdiff
path: root/src/algebra
diff options
context:
space:
mode:
Diffstat (limited to 'src/algebra')
-rw-r--r--src/algebra/sttaylor.spad.pamphlet109
1 files changed, 101 insertions, 8 deletions
diff --git a/src/algebra/sttaylor.spad.pamphlet b/src/algebra/sttaylor.spad.pamphlet
index a16e32ce..6d2a2afd 100644
--- a/src/algebra/sttaylor.spad.pamphlet
+++ b/src/algebra/sttaylor.spad.pamphlet
@@ -401,6 +401,14 @@ StreamTaylorSeriesOperations(A): Exports == Implementation where
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 @@ StreamTaylorSeriesOperations(A): Exports == Implementation where
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 @@ StreamTaylorSeriesOperations(A): Exports == Implementation where
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,24 +477,71 @@ StreamTaylorSeriesOperations(A): Exports == Implementation where
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"
+@
+
+[[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) =>
+ not negative?(num := numer rn) =>
-- It seems that this cannot happen, but I don't know why
- (co**num::NNI) * YS(powerrn(rn,(invCo :: A) * x,#1))
- (invCo :: A)**((-num)::NNI) * YS(powerrn(rn,(invCo :: A) * x,#1))
+ (co**num::NNI) * YS(powerrn(rn, (invCo :: A) * x, #1))
+ (invCo::A)**((-num)::NNI) * YS(powerrn(rn, (invCo :: A) * x, #1))
- RATPOWERS => co**rn * YS(powerrn(rn,(invCo :: A) * x,#1))
+ 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"