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

)abbrev domain UTS UnivariateTaylorSeries
++ Author: Clifton J. Williamson
++ Date Created: 21 December 1989
++ Date Last Updated: 21 September 1993
++ Description:
++ Dense Taylor series in one variable
++ \spadtype{UnivariateTaylorSeries} is a domain representing Taylor
++ series in
++ one variable with coefficients in an arbitrary ring.  The parameters
++ of the type specify the coefficient ring, the power series variable,
++ and the center of the power series expansion.  For example,
++ \spadtype{UnivariateTaylorSeries}(Integer,x,3) represents
++ Taylor series in
++ \spad{(x - 3)} with \spadtype{Integer} coefficients.

UnivariateTaylorSeries(Coef,var,cen) : SIG == CODE where
  Coef : Ring
  var : Symbol
  cen : Coef

  I    ==> Integer
  NNI  ==> NonNegativeInteger
  P    ==> Polynomial Coef
  RN   ==> Fraction Integer
  ST   ==> Stream
  STT  ==> StreamTaylorSeriesOperations Coef
  TERM ==> Record(k:NNI,c:Coef)
  UP   ==> UnivariatePolynomial(var,Coef)

  SIG ==> UnivariateTaylorSeriesCategory(Coef) with

    coerce : UP -> %
      ++\spad{coerce(p)} converts a univariate polynomial p in the variable
      ++\spad{var} to a univariate Taylor series in \spad{var}.

    univariatePolynomial : (%,NNI) -> UP
      ++\spad{univariatePolynomial(f,k)} returns a univariate polynomial
      ++ consisting of the sum of all terms of f of degree \spad{<= k}.

    coerce : Variable(var) -> %
      ++\spad{coerce(var)} converts the series variable \spad{var} into a
      ++ Taylor series.

    differentiate : (%,Variable(var)) -> %
      ++ \spad{differentiate(f(x),x)} computes the derivative of
      ++ \spad{f(x)} with respect to \spad{x}.

    lagrange : % -> %
      ++\spad{lagrange(g(x))} produces the Taylor series for \spad{f(x)}
      ++ where \spad{f(x)} is implicitly defined as \spad{f(x) = x*g(f(x))}.

    lambert : % -> %
      ++\spad{lambert(f(x))} returns \spad{f(x) + f(x^2) + f(x^3) + ...}.
      ++ This function is used for computing infinite products.
      ++ \spad{f(x)} should have zero constant coefficient.
      ++ If \spad{f(x)} is a Taylor series with constant term 1, then
      ++ \spad{product(n = 1..infinity,f(x^n)) = exp(log(lambert(f(x))))}.

    oddlambert : % -> %
      ++\spad{oddlambert(f(x))} returns \spad{f(x) + f(x^3) + f(x^5) + ...}.
      ++ \spad{f(x)} should have a zero constant coefficient.
      ++ This function is used for computing infinite products.
      ++ If \spad{f(x)} is a Taylor series with constant term 1, then
      ++ \spad{product(n=1..infinity,f(x^(2*n-1)))=exp(log(oddlambert(f(x))))}.

    evenlambert : % -> %
      ++\spad{evenlambert(f(x))} returns \spad{f(x^2) + f(x^4) + f(x^6) + ...}.
      ++ \spad{f(x)} should have a zero constant coefficient.
      ++ This function is used for computing infinite products.
      ++ If \spad{f(x)} is a Taylor series with constant term 1, then
      ++ \spad{product(n=1..infinity,f(x^(2*n))) = exp(log(evenlambert(f(x))))}

    generalLambert : (%,I,I) -> %
      ++\spad{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.

    revert : % -> %
      ++ \spad{revert(f(x))} returns a Taylor series \spad{g(x)} such that
      ++ \spad{f(g(x)) = g(f(x)) = x}. Series \spad{f(x)} should have constant
      ++ coefficient 0 and 1st order coefficient 1.

    multisect : (I,I,%) -> %
      ++\spad{multisect(a,b,f(x))} selects the coefficients of
      ++ \spad{x^((a+b)*n+a)}, and changes this monomial to \spad{x^n}.

    invmultisect : (I,I,%) -> %
      ++\spad{invmultisect(a,b,f(x))} substitutes \spad{x^((a+b)*n)}
      ++ for \spad{x^n} and multiples by \spad{x^b}.

    if Coef has Algebra Fraction Integer then

      integrate : (%,Variable(var)) -> %
        ++ \spad{integrate(f(x),x)} returns an anti-derivative of the power
        ++ series \spad{f(x)} with constant coefficient 0.
        ++ We may integrate a series when we can divide coefficients
        ++ by integers.

  CODE ==> InnerTaylorSeries(Coef) add

    Rep := Stream Coef

--% creation and destruction of series

    stream: % -> Stream Coef
    stream x  == x pretend Stream(Coef)

    coerce(v:Variable(var)) ==
      zero? cen => monomial(1,1)
      monomial(1,1) + monomial(cen,0)

    coerce(n:I)    == n :: Coef :: %

    coerce(r:Coef) == coerce(r)$STT

    monomial(c,n)  == monom(c,n)$STT

    getExpon: TERM -> NNI
    getExpon term == term.k

    getCoef: TERM -> Coef
    getCoef term == term.c

    rec: (NNI,Coef) -> TERM
    rec(expon,coef) == [expon,coef]

    recs: (ST Coef,NNI) -> ST TERM
    recs(st,n) == delay$ST(TERM)
      empty? st => empty()
      zero? (coef := frst st) => recs(rst st,n + 1)
      concat(rec(n,coef),recs(rst st,n + 1))

    terms x == recs(stream x,0)

    recsToCoefs: (ST TERM,NNI) -> ST Coef
    recsToCoefs(st,n) == delay
      empty? st => empty()
      term := frst st; expon := getExpon term
      n = expon => concat(getCoef term,recsToCoefs(rst st,n + 1))
      concat(0,recsToCoefs(st,n + 1))

    series(st: ST TERM) == recsToCoefs(st,0)

    stToPoly: (ST Coef,P,NNI,NNI) -> P
    stToPoly(st,term,n,n0) ==
      (n > n0) or (empty? st) => 0
      frst(st) * term ** n + stToPoly(rst st,term,n + 1,n0)

    polynomial(x,n) == stToPoly(stream x,(var :: P) - (cen :: P),0,n)

    polynomial(x,n1,n2) ==
      if n1 > n2 then (n1,n2) := (n2,n1)
      stToPoly(rest(stream x,n1),(var :: P) - (cen :: P),n1,n2)

    stToUPoly: (ST Coef,UP,NNI,NNI) -> UP
    stToUPoly(st,term,n,n0) ==
      (n > n0) or (empty? st) => 0
      frst(st) * term ** n + stToUPoly(rst st,term,n + 1,n0)

    univariatePolynomial(x,n) ==
      stToUPoly(stream x,monomial(1,1)$UP - monomial(cen,0)$UP,0,n)

    coerce(p:UP) ==
      zero? p => 0
      if not zero? cen then
        p := p(monomial(1,1)$UP + monomial(cen,0)$UP)
      st : ST Coef := empty()
      oldDeg : NNI := degree(p) + 1
      while not zero? p repeat
        deg := degree p
        delta := (oldDeg - deg - 1) :: NNI
        for i in 1..delta repeat st := concat(0$Coef,st)
        st := concat(leadingCoefficient p,st)
        oldDeg := deg; p := reductum p
      for i in 1..oldDeg repeat st := concat(0$Coef,st)
      st

    if Coef has coerce: Symbol -> Coef then
      if Coef has "**": (Coef,NNI) -> Coef then

        stToCoef: (ST Coef,Coef,NNI,NNI) -> Coef
        stToCoef(st,term,n,n0) ==
          (n > n0) or (empty? st) => 0
          frst(st) * term ** n + stToCoef(rst st,term,n + 1,n0)

        approximate(x,n) ==
          stToCoef(stream x,(var :: Coef) - cen,0,n)

--% values

    variable x == var

    center   s == cen

    coefficient(x,n) ==
       -- Cannot use elt!  Should return 0 if stream doesn't have it.
       u := stream x
       while not empty? u and n > 0 repeat
         u := rst u
         n := (n - 1) :: NNI
       empty? u or n ^= 0 => 0
       frst u

    elt(x:%,n:NNI) == coefficient(x,n)

--% functions

    map(f,x) == map(f,x)$Rep

    eval(x:%,r:Coef) == eval(stream x,r-cen)$STT

    differentiate x == deriv(stream x)$STT

    differentiate(x:%,v:Variable(var)) == differentiate x

    if Coef has PartialDifferentialRing(Symbol) then

      differentiate(x:%,s:Symbol) ==
        (s = variable(x)) => differentiate x
        map(y +-> differentiate(y,s),x) 
              - differentiate(center x,s)*differentiate(x)

    multiplyCoefficients(f,x) == gderiv(f,stream x)$STT

    lagrange x == lagrange(stream x)$STT

    lambert x == lambert(stream x)$STT

    oddlambert x == oddlambert(stream x)$STT

    evenlambert x == evenlambert(stream x)$STT

    generalLambert(x:%,a:I,d:I) == generalLambert(stream x,a,d)$STT

    extend(x,n) == extend(x,n+1)$Rep

    complete x == complete(x)$Rep

    truncate(x,n) == first(stream x,n + 1)$Rep

    truncate(x,n1,n2) ==
      if n2 < n1 then (n1,n2) := (n2,n1)
      m := (n2 - n1) :: NNI
      st := first(rest(stream x,n1)$Rep,m + 1)$Rep
      for i in 1..n1 repeat st := concat(0$Coef,st)
      st

    elt(x:%,y:%) == compose(stream x,stream y)$STT

    revert x == revert(stream x)$STT

    multisect(a,b,x) == multisect(a,b,stream x)$STT

    invmultisect(a,b,x) == invmultisect(a,b,stream x)$STT

    multiplyExponents(x,n) == invmultisect(n,0,x)

    quoByVar x == (empty? x => 0; rst x)

    if Coef has IntegralDomain then
      unit? x == unit? coefficient(x,0)
    if Coef has Field then
      if Coef is RN then

        (x:%) ** (s:Coef) == powern(s,stream x)$STT

      else

        (x:%) ** (s:Coef) == power(s,stream x)$STT

    if Coef has Algebra Fraction Integer then

      coerce(r:RN) == r :: Coef :: %

      integrate x == integrate(0,stream x)$STT

      integrate(x:%,v:Variable(var)) == integrate x

      if Coef has integrate: (Coef,Symbol) -> Coef and _
         Coef has variables: Coef -> List Symbol then

        integrate(x:%,s:Symbol) ==
          (s = variable(x)) => integrate x
          not entry?(s,variables center x) => map(y +-> integrate(y,s),x)
          error "integrate: center is a function of variable of integration"

      if Coef has TranscendentalFunctionCategory and _
         Coef has PrimitiveFunctionCategory and _
         Coef has AlgebraicallyClosedFunctionSpace Integer then

        integrateWithOneAnswer: (Coef,Symbol) -> Coef
        integrateWithOneAnswer(f,s) ==
          res := integrate(f,s)$FunctionSpaceIntegration(I,Coef)
          res case Coef => res :: Coef
          first(res :: List Coef)

        integrate(x:%,s:Symbol) ==
          (s = variable(x)) => integrate x
          not entry?(s,variables center x) =>
            map(y +-> integrateWithOneAnswer(y,s),x)
          error "integrate: center is a function of variable of integration"