axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / STTAYLOR.spad

)abbrev package STTAYLOR StreamTaylorSeriesOperations
++ Author: William Burge, Stephen Watt, Clifton J. Williamson
++ Date Created: 1986
++ Date Last Updated: 26 May 1994
++ Description:
++ StreamTaylorSeriesOperations implements Taylor series arithmetic,
++ where a Taylor series is represented by a stream of its coefficients.

StreamTaylorSeriesOperations(A) : SIG == CODE 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")

  SIG ==> 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}.

  CODE ==> 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(z +-> frst x+z, y)
      eq?(y,rst y) => map(z +-> frst y+z, x)
      concat(frst x + frst y,rst x + rst y)

    x - y == delay
      empty? y => x
      empty? x => -y
      eq?(x,rst x) => map(z +-> frst x-z, y)
      eq?(y,rst y) => map(z +-> z-frst y, x)
      concat(frst x - frst y,rst x - rst y)

    -y == map(z +-> -z, 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(z +-> s*z, x)

    (x:ST A) * (s:A) ==
      zero? s => zro()
      map(z +-> z*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(z +-> z*(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(z +-> z*(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(y +-> y+1, zz)

    integers z    == generate(y +-> y+1, z)

    oddintegers z == generate(y +-> y+2, z)

    int s         == generate(y +-> y+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,(y,z) +-> y+z,mapmult(x,generate(y +-> at*y,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(y +-> lagrangere(x,y))

    revert x ==
      empty? x => error "revert should start 0,1,..."
      zero? frst x =>
        empty? rst x => error "revert: should start 0,1,..."
        (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..] forms the
-- coefficients A[n] which is the sum of a[i] for all divisors i of n
-- (including 1 and n)

    --                               ---------
    -- returns the repeating stream [s,0,...,0]; (there are z zeroes)
    rptg1:(I,A) -> ST A
    rptg1(z,s) == repeating concat(s,rp(z,0))

    --                                       ---------
    -- returns the repeating stream [0,...,0,s,0,...,0]
    -- there are z leading zeroes and z-1 in the period
    rptg2:(I,A) -> ST A
    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((x,y) +-> rptg3(a,d,x,y),
                 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(y +-> (rst x)*y,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(y +-> nldere(lslsa,y))

      RATPOWERS : Boolean := A has "**": (A,RN) -> A

      smult: (RN,ST A) -> ST A
      smult(rn,x) == map(y +-> rn*y, 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"
        power :=
          (co = 1) => YS(y +-> powerrn(rn,x,y))
          (denom rn) = 1 =>
            not negative?(num := numer rn) => 
              (co**num::NNI) * YS(y +-> powerrn(rn,(invCo :: A) * x, y))
            (invCo::A)**((-num)::NNI) * YS(y +-> powerrn(rn,(invCo :: A)*x, y))
          RATPOWERS => co**rn * YS(y +-> powerrn(rn,(invCo :: A)*x, y))
          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(y +-> powerre(s,x,y))