path: root/src/algebra/taylor.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra taylor.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain ITAYLOR InnerTaylorSeries}
<<domain ITAYLOR InnerTaylorSeries>>=
)abbrev domain ITAYLOR InnerTaylorSeries
++ Author: Clifton J. Williamson
++ Date Created: 21 December 1989
++ Date Last Updated: October 26, 2009
++ Basic Operations:
++ Related Domains: UnivariateTaylorSeries(Coef,var,cen)
++ Also See:
++ AMS Classifications:
++ Keywords: stream, dense Taylor series
++ Examples:
++ References:
++ Description: Internal package for dense Taylor series.
++ This is an internal Taylor series type in which Taylor series
++ are represented by a \spadtype{Stream} of \spadtype{Ring} elements.
++ For univariate series, the \spad{Stream} elements are the Taylor
++ coefficients. For multivariate series, the \spad{n}th Stream element
++ is a form of degree n in the power series variables.

InnerTaylorSeries(Coef): Exports == Implementation where
  Coef  : Ring
  I   ==> Integer
  NNI ==> NonNegativeInteger
  ST  ==> Stream Coef
  STT ==> StreamTaylorSeriesOperations Coef

  Exports == Join(Ring, BiModule(Coef,Coef)) with
    coefficients: % -> Stream Coef
      ++\spad{coefficients(x)} returns a stream of ring elements.
      ++ When x is a univariate series, this is a stream of Taylor
      ++ coefficients. When x is a multivariate series, the
      ++ \spad{n}th element of the stream is a form of
      ++ degree n in the power series variables.
    series: Stream Coef -> %
      ++\spad{series(s)} creates a power series from a stream of
      ++ ring elements.
      ++ For univariate series types, the stream s should be a stream
      ++ of Taylor coefficients. For multivariate series types, the
      ++ stream s should be a stream of forms the \spad{n}th element
      ++ of which is a
      ++ form of degree n in the power series variables.
    pole?: % -> Boolean
      ++\spad{pole?(x)} tests if the series x has a pole.
      ++ Note: this is false when x is a Taylor series.
    order: % -> NNI
      ++\spad{order(x)} returns the order of a power series x,
      ++ i.e. the degree of the first non-zero term of the series.
    order: (%,NNI) -> NNI
      ++\spad{order(x,n)} returns the minimum of n and the order of x.
    * : (%,Integer)->%
      ++\spad{x*i} returns the product of integer i and the series x.
    if Coef has IntegralDomain then IntegralDomain
      --++ An IntegralDomain provides 'exquo'

  Implementation == add

    Rep == ST

    -- In what follows, we will be calling operations on Streams
    -- which are NOT defined in the package Stream.  Thus, it is
    -- necessary to explicitly pass back and forth between Rep and %.

    series st == per st

    0 == per coerce(0)$STT
    1 == per coerce(1)$STT

    x = y ==
      -- tests if two power series are equal
      -- difference must be a finite stream of zeroes of length <= n + 1,
      -- where n = $streamCount$Lisp
      st : ST := rep(x - y)
      n : I := _$streamCount$Lisp
      for i in 0..n repeat
        empty? st => return true
        frst st ~= 0 => return false
        st := rst st
      empty? st

    coefficients x == rep x

    x + y            == per(rep x +$STT rep y)
    x - y            == per(rep x -$STT rep y)
    (x:%) * (y:%)    == per(rep x *$STT rep y)
    - x              == per(-$STT rep x)
    (i:I) * (x:%)    == per((i::Coef) *$STT rep x)
    (x:%) * (i:I)    == per(rep x *$STT (i::Coef))
    (c:Coef) * (x:%) == per(c *$STT rep x)
    (x:%) * (c:Coef) == per(rep x *$STT c)

    recip x ==
      (rec := recip$STT rep x) case "failed" => "failed"
      per(rec :: ST)

    if Coef has IntegralDomain then

      x exquo y ==
        (quot := rep x exquo$STT rep y) case "failed" => "failed"
        per(quot :: ST)

    x:% ** n:NNI ==
      n = 0 => 1
      expt(x,n :: PositiveInteger)$RepeatedSquaring(%)

    characteristic == characteristic$Coef
    pole? x == false

    iOrder: (ST,NNI,NNI) -> NNI
    iOrder(st,n,n0) ==
      (n = n0) or (empty? st) => n0
      zero? frst st => iOrder(rst st,n + 1,n0)
      n

    order(x,n) == iOrder(rep x,0,n)

    iOrder2: (ST,NNI) -> NNI
    iOrder2(st,n) ==
      empty? st => error "order: series has infinite order"
      zero? frst st => iOrder2(rst st,n + 1)
      n

    order x == iOrder2(rep x,0)

    zero? x == empty? rep x

@
\section{domain UTS UnivariateTaylorSeries}
<<domain UTS UnivariateTaylorSeries>>=
)abbrev domain UTS UnivariateTaylorSeries
++ Author: Clifton J. Williamson
++ Date Created: 21 December 1989
++ Date Last Updated: June 18, 2010
++ Basic Operations:
++ Related Domains: UnivariateLaurentSeries(Coef,var,cen), UnivariatePuiseuxSeries(Coef,var,cen)
++ Also See:
++ AMS Classifications:
++ Keywords: dense, Taylor series
++ Examples:
++ References:
++ 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): Exports == Implementation 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)
  Exports ==> Join(UnivariateTaylorSeriesCategory(Coef),_
                   PartialDifferentialDomain(%,Variable var)) 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.
    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 invertible 1st order coefficient.
    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.

  Implementation ==> 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(differentiate(#1,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(integrate(#1,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(integrateWithOneAnswer(#1,s),x)
          error "integrate: center is a function of variable of integration"

--% OutputForms
--  We use the default coerce: % -> OutputForm in UTSCAT&

@
\section{package UTS2 UnivariateTaylorSeriesFunctions2}
<<package UTS2 UnivariateTaylorSeriesFunctions2>>=
)abbrev package UTS2 UnivariateTaylorSeriesFunctions2
++ Author: Clifton J. Williamson
++ Date Created: 9 February 1990
++ Date Last Updated: 9 February 1990
++ Basic Operations:
++ Related Domains: UnivariateTaylorSeries(Coef1,var,cen)
++ Also See:
++ AMS Classifications:
++ Keywords: Taylor series, map
++ Examples:
++ References:
++ Description: Mapping package for univariate Taylor series.
++   This package allows one to apply a function to the coefficients of
++   a univariate Taylor series.
UnivariateTaylorSeriesFunctions2(Coef1,Coef2,UTS1,UTS2):_
 Exports == Implementation where
  Coef1 : Ring
  Coef2 : Ring
  UTS1  : UnivariateTaylorSeriesCategory Coef1
  UTS2  : UnivariateTaylorSeriesCategory Coef2
  ST2 ==> StreamFunctions2(Coef1,Coef2)

  Exports ==> with
    map: (Coef1 -> Coef2,UTS1) -> UTS2
      ++\spad{map(f,g(x))} applies the map f to the coefficients of
      ++ the Taylor series \spad{g(x)}.

  Implementation ==> add

    map(f,uts) == series map(f,coefficients uts)$ST2

@
\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>>

<<domain ITAYLOR InnerTaylorSeries>>
<<domain UTS UnivariateTaylorSeries>>
<<package UTS2 UnivariateTaylorSeriesFunctions2>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}