--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.

)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 should be invertible.
    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) ==
      negative? z => 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 and not empty? rst x =>
        case recip rst x is
          y@ST(A) => lagrange y
          otherwise => error "revert: should start 0,a,... with invertible a"
      error "revert: argument is not reversible"

--% 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"
        if positive? ord and negative? rn then
           error "**: negative power does not exist"
        co := frst x
        (invCo := recip co) case "failed" =>
           error "** rational power of coefficient undefined"
        power :=
          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
              (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"
        monom(1, (ord :: I) * numer(rn)) * power
    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()
        not one? frst x => 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))

