diff options
Diffstat (limited to 'src/algebra/sttaylor.spad.pamphlet')
| -rw-r--r-- | src/algebra/sttaylor.spad.pamphlet | 515 |
1 files changed, 515 insertions, 0 deletions
diff --git a/src/algebra/sttaylor.spad.pamphlet b/src/algebra/sttaylor.spad.pamphlet new file mode 100644 index 00000000..27a16ec0 --- /dev/null +++ b/src/algebra/sttaylor.spad.pamphlet @@ -0,0 +1,515 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra sttaylor.spad} +\author{William Burge, Stephen Watt, Clifton J. Williamson} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{package STTAYLOR StreamTaylorSeriesOperations} +Problems raising a UTS to a negative integer power. + +The code in [[powern(rn,x)]] which raises an unnecessary error +where no distinction between rational and integer powers are made. + +The fix is easy. Since the problem does not exist in SUPS we can +just take the definition there. + +<<package STTAYLOR StreamTaylorSeriesOperations>>= +)abbrev package STTAYLOR StreamTaylorSeriesOperations +++ Author: William Burge, Stephen Watt, Clifton J. Williamson +++ Date Created: 1986 +++ Date Last Updated: 26 May 1994 +++ Basic Operations: +++ Related Domains: Stream(A), ParadoxicalCombinatorsForStreams(A), +++ StreamTranscendentalFunctions(A), +++ StreamTranscendentalFunctionsNonCommutative(A) +++ Also See: +++ AMS Classifications: +++ Keywords: stream, Taylor series +++ Examples: +++ References: +++ Description: +++ StreamTaylorSeriesOperations implements Taylor series arithmetic, +++ where a Taylor series is represented by a stream of its coefficients. +StreamTaylorSeriesOperations(A): Exports == Implementation where + A : Ring + RN ==> Fraction Integer + I ==> Integer + NNI ==> NonNegativeInteger + ST ==> Stream + SP2 ==> StreamFunctions2 + SP3 ==> StreamFunctions3 + L ==> List + LA ==> List A + YS ==> Y$ParadoxicalCombinatorsForStreams(A) + UN ==> Union(ST A,"failed") + Exports ==> with + "+" : (ST A,ST A) -> ST A + ++ a + b returns the power series sum of \spad{a} and \spad{b}: + ++ \spad{[a0,a1,..] + [b0,b1,..] = [a0 + b0,a1 + b1,..]} + "-" : (ST A,ST A) -> ST A + ++ a - b returns the power series difference of \spad{a} and + ++ \spad{b}: \spad{[a0,a1,..] - [b0,b1,..] = [a0 - b0,a1 - b1,..]} + "-" : ST A -> ST A + ++ - a returns the power series negative of \spad{a}: + ++ \spad{- [a0,a1,...] = [- a0,- a1,...]} + "*" : (ST A,ST A) -> ST A + ++ a * b returns the power series (Cauchy) product of \spad{a} and b: + ++ \spad{[a0,a1,...] * [b0,b1,...] = [c0,c1,...]} where + ++ \spad{ck = sum(i + j = k,ai * bk)}. + "*" : (A,ST A) -> ST A + ++ r * a returns the power series scalar multiplication of r by \spad{a}: + ++ \spad{r * [a0,a1,...] = [r * a0,r * a1,...]} + "*" : (ST A,A) -> ST A + ++ a * r returns the power series scalar multiplication of \spad{a} by r: + ++ \spad{[a0,a1,...] * r = [a0 * r,a1 * r,...]} + "exquo" : (ST A,ST A) -> Union(ST A,"failed") + ++ exquo(a,b) returns the power series quotient of \spad{a} by b, + ++ if the quotient exists, and "failed" otherwise + "/" : (ST A,ST A) -> ST A + ++ a / b returns the power series quotient of \spad{a} by b. + ++ An error message is returned if \spad{b} is not invertible. + ++ This function is used in fixed point computations. + recip : ST A -> UN + ++ recip(a) returns the power series reciprocal of \spad{a}, or + ++ "failed" if not possible. + monom : (A,I) -> ST A + ++ monom(deg,coef) is a monomial of degree deg with coefficient + ++ coef. + integers : I -> ST I + ++ integers(n) returns \spad{[n,n+1,n+2,...]}. + oddintegers : I -> ST I + ++ oddintegers(n) returns \spad{[n,n+2,n+4,...]}. + int : A -> ST A + ++ int(r) returns [r,r+1,r+2,...], where r is a ring element. + mapmult : (ST A,ST A) -> ST A + ++ mapmult([a0,a1,..],[b0,b1,..]) + ++ returns \spad{[a0*b0,a1*b1,..]}. + deriv : ST A -> ST A + ++ deriv(a) returns the derivative of the power series with + ++ respect to the power series variable. Thus + ++ \spad{deriv([a0,a1,a2,...])} returns \spad{[a1,2 a2,3 a3,...]}. + gderiv : (I -> A,ST A) -> ST A + ++ gderiv(f,[a0,a1,a2,..]) returns + ++ \spad{[f(0)*a0,f(1)*a1,f(2)*a2,..]}. + coerce : A -> ST A + ++ coerce(r) converts a ring element r to a stream with one element. + eval : (ST A,A) -> ST A + ++ eval(a,r) returns a stream of partial sums of the power series + ++ \spad{a} evaluated at the power series variable equal to r. + compose : (ST A,ST A) -> ST A + ++ compose(a,b) composes the power series \spad{a} with + ++ the power series b. + lagrange : ST A -> ST A + ++ lagrange(g) produces the power series for f where f is + ++ implicitly defined as \spad{f(z) = z*g(f(z))}. + revert : ST A -> ST A + ++ revert(a) computes the inverse of a power series \spad{a} + ++ with respect to composition. + ++ the series should have constant coefficient 0 and first + ++ order coefficient 1. + addiag : ST ST A -> ST A + ++ addiag(x) performs diagonal addition of a stream of streams. if x = + ++ \spad{[[a<0,0>,a<0,1>,..],[a<1,0>,a<1,1>,..],[a<2,0>,a<2,1>,..],..]} + ++ and \spad{addiag(x) = [b<0,b<1>,...], then b<k> = sum(i+j=k,a<i,j>)}. + lambert : ST A -> ST A + ++ lambert(st) computes \spad{f(x) + f(x**2) + f(x**3) + ...} + ++ if st is a stream representing \spad{f(x)}. + ++ This function is used for computing infinite products. + ++ If \spad{f(x)} is a power series with constant coefficient 1 then + ++ \spad{prod(f(x**n),n = 1..infinity) = exp(lambert(log(f(x))))}. + oddlambert : ST A -> ST A + ++ oddlambert(st) computes \spad{f(x) + f(x**3) + f(x**5) + ...} + ++ if st is a stream representing \spad{f(x)}. + ++ This function is used for computing infinite products. + ++ If f(x) is a power series with constant coefficient 1 then + ++ \spad{prod(f(x**(2*n-1)),n=1..infinity) = exp(oddlambert(log(f(x))))}. + evenlambert : ST A -> ST A + ++ evenlambert(st) computes \spad{f(x**2) + f(x**4) + f(x**6) + ...} + ++ if st is a stream representing \spad{f(x)}. + ++ This function is used for computing infinite products. + ++ If \spad{f(x)} is a power series with constant coefficient 1, then + ++ \spad{prod(f(x**(2*n)),n=1..infinity) = exp(evenlambert(log(f(x))))}. + generalLambert : (ST A,I,I) -> ST A + ++ generalLambert(f(x),a,d) returns + ++ \spad{f(x**a) + f(x**(a + d)) + f(x**(a + 2 d)) + ...}. + ++ \spad{f(x)} should have zero constant + ++ coefficient and \spad{a} and d should be positive. + multisect : (I,I,ST A) -> ST A + ++ multisect(a,b,st) + ++ selects the coefficients of \spad{x**((a+b)*n+a)}, + ++ and changes them to \spad{x**n}. + invmultisect : (I,I,ST A) -> ST A + ++ invmultisect(a,b,st) substitutes \spad{x**((a+b)*n)} for \spad{x**n} + ++ and multiplies by \spad{x**b}. + if A has Algebra RN then + integrate : (A,ST A) -> ST A + ++ integrate(r,a) returns the integral of the power series \spad{a} + ++ with respect to the power series variableintegration where + ++ r denotes the constant of integration. Thus + ++ \spad{integrate(a,[a0,a1,a2,...]) = [a,a0,a1/2,a2/3,...]}. + lazyIntegrate : (A,() -> ST A) -> ST A + ++ lazyIntegrate(r,f) is a local function + ++ used for fixed point computations. + nlde : ST ST A -> ST A + ++ nlde(u) solves a + ++ first order non-linear differential equation described by u of the + ++ form \spad{[[b<0,0>,b<0,1>,...],[b<1,0>,b<1,1>,.],...]}. + ++ the differential equation has the form + ++ \spad{y' = sum(i=0 to infinity,j=0 to infinity,b<i,j>*(x**i)*(y**j))}. + powern : (RN,ST A) -> ST A + ++ powern(r,f) raises power series f to the power r. + if A has Field then + mapdiv : (ST A,ST A) -> ST A + ++ mapdiv([a0,a1,..],[b0,b1,..]) returns + ++ \spad{[a0/b0,a1/b1,..]}. + lazyGintegrate : (I -> A,A,() -> ST A) -> ST A + ++ lazyGintegrate(f,r,g) is used for fixed point computations. + power : (A,ST A) -> ST A + ++ power(a,f) returns the power series f raised to the power \spad{a}. + + Implementation ==> add + +--% definitions + + zro: () -> ST A + -- returns a zero power series + zro() == empty()$ST(A) + +--% arithmetic + + x + y == delay + empty? y => x + empty? x => y + eq?(x,rst x) => map(frst x + #1,y) + eq?(y,rst y) => map(frst y + #1,x) + concat(frst x + frst y,rst x + rst y) + + x - y == delay + empty? y => x + empty? x => -y + eq?(x,rst x) => map(frst x - #1,y) + eq?(y,rst y) => map(#1 - frst y,x) + concat(frst x - frst y,rst x - rst y) + + -y == map(_-#1,y) + + (x:ST A) * (y:ST A) == delay + empty? y => zro() + empty? x => zro() + concat(frst x * frst y,frst x * rst y + rst x * y) + + (s:A) * (x:ST A) == + zero? s => zro() + map(s * #1,x) + + (x:ST A) * (s:A) == + zero? s => zro() + map(#1 * s,x) + + iDiv: (ST A,ST A,A) -> ST A + iDiv(x,y,ry0) == delay + empty? x => empty() + c0 := frst x * ry0 + concat(c0,iDiv(rst x - c0 * rst y,y,ry0)) + + x exquo y == + for n in 1.. repeat + n > 1000 => return "failed" + empty? y => return "failed" + empty? x => return empty() + frst y = 0 => + frst x = 0 => (x := rst x; y := rst y) + return "failed" + leave "first entry in y is non-zero" + (ry0 := recip frst y) case "failed" => "failed" + empty? rst y => map(#1 * (ry0 :: A),x) + iDiv(x,y,ry0 :: A) + + (x:ST A) / (y:ST A) == delay + empty? y => error "/: division by zero" + empty? x => empty() + (ry0 := recip frst y) case "failed" => + error "/: second argument is not invertible" + empty? rst y => map(#1 * (ry0 :: A),x) + iDiv(x,y,ry0 :: A) + + recip x == + empty? x => "failed" + rh1 := recip frst x + rh1 case "failed" => "failed" + rh := rh1 :: A + delay + concat(rh,iDiv(- rh * rst x,x,rh)) + +--% coefficients + + rp: (I,A) -> L A + -- rp(z,s) is a list of length z each of whose entries is s. + rp(z,s) == + z <= 0 => empty() + concat(s,rp(z-1,s)) + + rpSt: (I,A) -> ST A + -- rpSt(z,s) is a stream of length z each of whose entries is s. + rpSt(z,s) == delay + z <= 0 => empty() + concat(s,rpSt(z-1,s)) + + monom(s,z) == + z < 0 => error "monom: cannot create monomial of negative degree" + concat(rpSt(z,0),concat(s,zro())) + +--% some streams of integers + nnintegers: NNI -> ST NNI + nnintegers zz == generate(#1 + 1,zz) + integers z == generate(#1 + 1,z) + oddintegers z == generate(#1 + 2,z) + int s == generate(#1 + 1,s) + +--% derivatives + + mapmult(x,y) == delay + empty? y => zro() + empty? x => zro() + concat(frst x * frst y,mapmult(rst x,rst y)) + + deriv x == + empty? x => zro() + mapmult(int 1,rest x) + + gderiv(f,x) == + empty? x => zro() + mapmult(map(f,integers 0)$SP2(I,A),x) + +--% coercions + + coerce(s:A) == + zero? s => zro() + concat(s,zro()) + +--% evaluations and compositions + + eval(x,at) == scan(0,#1 + #2,mapmult(x,generate(at * #1,1)))$SP2(A,A) + + compose(x,y) == delay + empty? y => concat(frst x,zro()) + not zero? frst y => + error "compose: 2nd argument should have 0 constant coefficient" + empty? x => zro() + concat(frst x,compose(rst x,y) * rst(y)) + +--% reversion + + lagrangere:(ST A,ST A) -> ST A + lagrangere(x,c) == delay(concat(0,compose(x,c))) + lagrange x == YS(lagrangere(x,#1)) + + revert x == + empty? x => error "revert should start 0,1,..." + zero? frst x => + empty? rst x => error "revert: should start 0,1,..." +-- one? frst rst x => lagrange(recip(rst x) :: (ST A)) + (frst rst x) = 1 => lagrange(recip(rst x) :: (ST A)) + error "revert:should start 0,1,..." + +--% lambert functions + + addiag(ststa:ST ST A) == delay + empty? ststa => zro() + empty? frst ststa => concat(0,addiag rst ststa) + concat(frst(frst ststa),rst(frst ststa) + addiag(rst ststa)) + +-- lambert operates on a series +/[a[i]x**i for i in 1..] , and produces +-- the series +/[a[i](x**i/(1-x**i)) for i in 1..] i.e. forms the +-- coefficients A[n] which is the sum of a[i] for all divisors i of n +-- (including 1 and n) + + rptg1:(I,A) -> ST A + -- --------- + -- returns the repeating stream [s,0,...,0]; (there are z zeroes) + rptg1(z,s) == repeating concat(s,rp(z,0)) + + rptg2:(I,A) -> ST A + -- --------- + -- returns the repeating stream [0,...,0,s,0,...,0] + -- there are z leading zeroes and z-1 in the period + rptg2(z,s) == repeating concat(rp(z,0),concat(s,rp(z-1,0))) + + rptg3:(I,I,I,A) -> ST A + rptg3(a,d,n,s) == + concat(rpSt(n*(a-1),0),repeating(concat(s,rp(d*n-1,0)))) + + lambert x == delay + empty? x => zro() + zero? frst x => + concat(0,addiag(map(rptg1,integers 0,rst x)$SP3(I,A,ST A))) + error "lambert:constant coefficient should be zero" + + oddlambert x == delay + empty? x => zro() + zero? frst x => + concat(0,addiag(map(rptg1,oddintegers 1,rst x)$SP3(I,A,ST A))) + error "oddlambert: constant coefficient should be zero" + + evenlambert x == delay + empty? x => zro() + zero? frst x => + concat(0,addiag(map(rptg2,integers 1,rst x)$SP3(I,A,ST A))) + error "evenlambert: constant coefficient should be zero" + + generalLambert(st,a,d) == delay + a < 1 or d < 1 => + error "generalLambert: both integer arguments must be positive" + empty? st => zro() + zero? frst st => + concat(0,addiag(map(rptg3(a,d,#1,#2),_ + integers 1,rst st)$SP3(I,A,ST A))) + error "generalLambert: constant coefficient should be zero" + +--% misc. functions + + ms: (I,I,ST A) -> ST A + ms(m,n,s) == delay + empty? s => zro() + zero? n => concat(frst s,ms(m,m-1,rst s)) + ms(m,n-1,rst s) + + multisect(b,a,x) == ms(a+b,0,rest(x,a :: NNI)) + + altn: (ST A,ST A) -> ST A + altn(zs,s) == delay + empty? s => zro() + concat(frst s,concat(zs,altn(zs,rst s))) + + invmultisect(a,b,x) == + concat(rpSt(b,0),altn(rpSt(a + b - 1,0),x)) + +-- comps(ststa,y) forms the composition of +/b[i,j]*y**i*x**j +-- where y is a power series in y. + + cssa ==> concat$(ST ST A) + mapsa ==> map$SP2(ST A,ST A) + comps: (ST ST A,ST A) -> ST ST A + comps(ststa,x) == delay$(ST ST A) + empty? ststa => empty()$(ST ST A) + empty? x => cssa(frst ststa,empty()$(ST ST A)) + cssa(frst ststa,mapsa((rst x) * #1,comps(rst ststa,x))) + + if A has Algebra RN then + integre: (ST A,I) -> ST A + integre(x,n) == delay + empty? x => zro() + concat((1$I/n) * frst(x),integre(rst x,n + 1)) + + integ: ST A -> ST A + integ x == integre(x,1) + + integrate(a,x) == concat(a,integ x) + lazyIntegrate(s,xf) == concat(s,integ(delay xf)) + + nldere:(ST ST A,ST A) -> ST A + nldere(lslsa,c) == lazyIntegrate(0,addiag(comps(lslsa,c))) + nlde lslsa == YS(nldere(lslsa,#1)) + + RATPOWERS : Boolean := A has "**": (A,RN) -> A + + smult: (RN,ST A) -> ST A + smult(rn,x) == map(rn * #1,x) + 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) == + order : I := 0 + for n in 0.. repeat + empty? x => return zro() + not zero? frst x => (order := n; leave x) + x := rst x + n = 1000 => + error "**: series with many leading zero coefficients" + (ord := (order exquo denom(rn))) case "failed" => + error "**: rational power does not exist" + 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 + power := +-- one? co => YS(powerrn(rn,x,#1)) + (co = 1) => 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 + (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)) + error "** rational power of coefficient undefined" + + if A has Field then + mapdiv(x,y) == delay + empty? y => error "stream division by zero" + empty? x => zro() + concat(frst x/frst y,mapdiv(rst x,rst y)) + + ginteg: (I -> A,ST A) -> ST A + ginteg(f,x) == mapdiv(x,map(f,integers 1)$SP2(I,A)) + + lazyGintegrate(fntoa,s,xf) == concat(s,ginteg(fntoa,delay xf)) + + finteg: ST A -> ST A + finteg x == mapdiv(x,int 1) + powerre: (A,ST A,ST A) -> ST A + powerre(s,x,c) == delay + empty? x => zro() + frst x^=1 => error "**:constant coefficient should be 1" + concat(frst x,finteg((s+1)*(c*deriv x))-rst x * c) + power(s,x) == YS(powerre(s,x,#1)) + +@ +\section{License} +<<license>>= +--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +--All rights reserved. +-- +--Redistribution and use in source and binary forms, with or without +--modification, are permitted provided that the following conditions are +--met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +@ +<<*>>= +<<license>> + +<<package STTAYLOR StreamTaylorSeriesOperations>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |