--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.
 
import NonNegativeInteger
import List
import Stream
)abbrev domain SMTS SparseMultivariateTaylorSeries
++ This domain provides multivariate Taylor series
++ Authors: William Burge, Stephen Watt, Clifton Williamson
++ Date Created: 15 August 1988
++ Date Last Updated: 18 May 1991
++ Basic Operations:
++ Related Domains:
++ Also See: UnivariateTaylorSeries
++ AMS Classifications:
++ Keywords: multivariate, Taylor, series
++ Examples:
++ References:
++ Description:
++   This domain provides multivariate Taylor series with variables
++   from an arbitrary ordered set.  A Taylor series is represented
++   by a stream of polynomials from the polynomial domain SMP.
++   The nth element of the stream is a form of degree n.  SMTS is an
++   internal domain.
SparseMultivariateTaylorSeries(Coef,Var,SMP):_
 Exports == Implementation where
  Coef : Ring
  Var  : OrderedSet
  SMP  : PolynomialCategory(Coef,IndexedExponents Var,Var)
  I   ==> Integer
  L   ==> List
  NNI ==> NonNegativeInteger
  OUT ==> OutputForm
  PS  ==> InnerTaylorSeries SMP
  RN  ==> Fraction Integer
  ST  ==> Stream
  StS ==> Stream SMP
  STT ==> StreamTaylorSeriesOperations SMP
  STF ==> StreamTranscendentalFunctions SMP
  ST2 ==> StreamFunctions2
  ST3 ==> StreamFunctions3
 
  Exports ==> MultivariateTaylorSeriesCategory(Coef,Var) with
    coefficient: (%,NNI) -> SMP
      ++ \spad{coefficient(s, n)} gives the terms of total degree n.
    coerce: Var -> %
      ++ \spad{coerce(var)} converts a variable to a Taylor series
    coerce: SMP -> %
      ++ \spad{coerce(poly)} regroups the terms by total degree and forms
      ++ a series.
    *:(SMP,%)->%
      ++\spad{smp*ts} multiplies a TaylorSeries by a monomial SMP.
    csubst:(L Var,L StS) -> (SMP -> StS)
      ++\spad{csubst(a,b)} is for internal use only
 
    if Coef has Algebra Fraction Integer then
      integrate: (%,Var,Coef) -> %
        ++\spad{integrate(s,v,c)} is the integral of s with respect
        ++ to v and having c as the constant of integration.
      fintegrate: (() -> %,Var,Coef) -> %
        ++\spad{fintegrate(f,v,c)} is the integral of \spad{f()} with respect
        ++ to v and having c as the constant of integration.
        ++ The evaluation of \spad{f()} is delayed.
 
  Implementation ==> PS add
 
    Rep := StS -- Below we use the fact that Rep of PS is Stream SMP.
    -- FIXME: The next two functions come free if assignment to Rep is removed.
    rep(x: %): Rep == x pretend Rep
    per(x: Rep): % == x pretend %
    extend(x,n) == extend(x,n + 1)$Rep
    complete x == complete(x)$Rep

    evalstream:(%,L Var,L SMP) -> StS
    evalstream(s:%,lv:(L Var),lsmp:(L SMP)):(ST SMP) ==
      scan(0,_+$SMP,map(eval(#1,lv,lsmp),s pretend StS))$ST2(SMP,SMP)
 
    addvariable:(Var,InnerTaylorSeries Coef) -> %
    addvariable(v,s) ==
      ints := integers(0)$STT pretend ST NNI
      map(monomial(#2 :: SMP,v,#1)$SMP,ints,s pretend ST Coef)$ST3(NNI,Coef,SMP)
 
    coefficient(s:%,n: NNI) == 
      elt(s,n + 1)$Rep  -- 1-based indexing for streams
 
--% creation of series
 
    coerce(r:Coef) == monom(r::SMP,0)$STT
    smp:SMP * p:% == per((smp *  rep p)$STT)
    r:Coef * p:% == per(((r::SMP) *  rep p )$STT)
    p:% * r:Coef == per(((r::SMP) * rep p)$STT)
    mts(p:SMP):% ==
      (uv := mainVariable p) case "failed" => monom(p,0)$STT
      v := uv :: Var
      s : % := 0
      up := univariate(p,v)
      while not zero? up repeat
        s := s + monomial(1,v,degree up) * mts(leadingCoefficient up)
        up := reductum up
      s
 
    coerce(p:SMP) == mts p
    coerce(v:Var) == v :: SMP :: %
 
    monomial(r:%,v:Var,n:NNI) ==
      r * monom(monomial(1,v,n)$SMP,n)$STT
 
--% evaluation
 
    substvar: (SMP,L Var,L %) -> %
    substvar(p,vl,q) ==
      null vl => monom(p,0)$STT
      (uv := mainVariable p) case "failed" => monom(p,0)$STT
      v := uv :: Var
      v = first vl =>
        s : % := 0
        up := univariate(p,v)
        while not zero? up repeat
          c := leadingCoefficient up
          s := s + first q ** degree up * substvar(c,rest vl,rest q)
          up := reductum up
        s
      substvar(p,rest vl,rest q)
 
    sortmfirst:(SMP,L Var,L %) -> %
    sortmfirst(p,vl,q) ==
      nlv : L Var := sort(#1 > #2,vl)
      nq : L % := [q position$(L Var) (i,vl) for i in nlv]
      substvar(p,nlv,nq)
 
    csubst(vl,q) == sortmfirst(#1,vl,q pretend L(%)) pretend StS
 
    restCheck(s:StS):StS ==
      -- checks that stream is null or first element is 0
      -- returns empty() or rest of stream
      empty? s => s
      not zero? frst s =>
        error "eval: constant coefficient should be 0"
      rst s
 
    eval(s:%,v:L Var,q:L %) ==
      #v ~= #q =>
        error "eval: number of variables should equal number of values"
      nq : L StS := [restCheck(i pretend StS) for i in q]
      addiag(map(csubst(v,nq),s pretend StS)$ST2(SMP,StS))$STT pretend %
 
    substmts(v:Var,p:SMP,q:%):% ==
      up := univariate(p,v)
      ss : % := 0
      while not zero? up repeat
        d:=degree up
        c:SMP:=leadingCoefficient up
        ss := ss + c* q ** d
        up := reductum up
      ss
 
    subststream(v:Var,p:SMP,q:StS):StS==
      substmts(v,p,q pretend %) pretend StS
 
    comp1:(Var,StS,StS) -> StS
    comp1(v,r,t)== addiag(map(subststream(v,#1,t),r)$ST2(SMP,StS))$STT
 
    comp(v:Var,s:StS,t:StS):StS == delay
      empty? s => s
      f := frst s; r : StS := rst s;
      empty? r => s
      empty? t => concat(f,comp1(v,r,empty()$StS))
      not zero? frst t =>
        error "eval: constant coefficient should be zero"
      concat(f,comp1(v,r,rst t))
 
    eval(s:%,v:Var,t:%) == comp(v,s pretend StS,t pretend StS)
 
--% differentiation and integration
 
    differentiate(s:%,v:Var):% ==
      empty? s => 0
      map(differentiate(#1,v),rst s)
 
    if Coef has Algebra Fraction Integer then
      (x:%) ** (r:RN) == powern(r,rep x)$STT
      (r:RN) * (x:%)  == per map(r * #1, rep x)$ST2(SMP,SMP)
      (x:%) * (r:RN)  == per map(#1 * r,rep x )$ST2(SMP,SMP)
 
      exp x == exp(rep x)$STF
      log x == log(rep x)$STF
 
      sin x == sin(rep x)$STF
      cos x == cos(rep x)$STF
      tan x == tan(rep x)$STF
      cot x == cot(rep x)$STF
      sec x == sec(rep x)$STF
      csc x == csc(rep x)$STF
 
      asin x == asin(rep x)$STF
      acos x == acos(rep x)$STF
      atan x == atan(rep x)$STF
      acot x == acot(rep x)$STF
      asec x == asec(rep x)$STF
      acsc x == acsc(rep x)$STF
 
      sinh x == sinh(rep x)$STF
      cosh x == cosh(rep x)$STF
      tanh x == tanh(rep x)$STF
      coth x == coth(rep x)$STF
      sech x == sech(rep x)$STF
      csch x == csch(rep x)$STF
 
      asinh x == asinh(rep x)$STF
      acosh x == acosh(rep x)$STF
      atanh x == atanh(rep x)$STF
      acoth x == acoth(rep x)$STF
      asech x == asech(rep x)$STF
      acsch x == acsch(rep x)$STF
 
      intsmp(v:Var,p: SMP): SMP ==
        up := univariate(p,v)
        ss : SMP := 0
        while not zero? up repeat
          d := degree up
          c := leadingCoefficient up
          ss := ss + inv((d+1) :: RN) * monomial(c,v,d+1)$SMP
          up := reductum up
        ss
 
      fintegrate(f,v,r) ==
        concat(r::SMP,delay map(intsmp(v,#1),f() pretend StS))
      integrate(s,v,r) ==
        concat(r::SMP,map(intsmp(v,#1),s pretend StS))
 
    -- If there is more than one term of the same order, group them.
    tout(p:SMP):OUT ==
      pe := p :: OUT
      monomial? p => pe
      paren pe
 
    showAll?: () -> Boolean
    -- check a global Lisp variable
    showAll?() == true
 
    coerce(s:%):OUT ==
      uu := s pretend Stream(SMP)
      empty? uu => (0$SMP) :: OUT
      count : NNI := _$streamCount$Lisp
      l : List OUT := empty()
      n : NNI := 0
      while n <= count and not empty? uu repeat
        if frst(uu) ~= 0 then l := concat(tout frst uu,l)
        uu := rst uu
        n := n + 1
      if showAll?() then
        while explicitEntries? uu and not eq?(uu,rst uu) repeat
          if frst(uu) ~= 0 then l := concat(tout frst uu,l)
          uu := rst uu
          n := n + 1
      l :=
        explicitlyEmpty? uu => l
        eq?(uu,rst uu) and frst uu = 0 => l
        concat(prefix("O" :: OUT,[n :: OUT]),l)
      empty? l => (0$SMP) :: OUT
      reduce("+",reverse! l)
    if Coef has Field then
      p:% / r:Coef == per(map(#1/$SMP r,rep p)$StreamFunctions2(SMP,SMP))

)abbrev domain TS TaylorSeries
++ Authors: Burge, Watt, Williamson
++ Date Created: 15 August 1988
++ Date Last Updated: 18 May 1991
++ Basic Operations:
++ Related Domains: SparseMultivariateTaylorSeries
++ Also See: UnivariateTaylorSeries
++ AMS Classifications:
++ Keywords: multivariate, Taylor, series
++ Examples:
++ References:
++ Description:
++   \spadtype{TaylorSeries} is a general multivariate Taylor series domain
++   over the ring Coef and with variables of type Symbol.
TaylorSeries(Coef): Exports == Implementation where
  Coef  : Ring
  L   ==> List
  NNI ==> NonNegativeInteger
  SMP ==> Polynomial Coef
  StS ==> Stream SMP
 
  Exports ==> MultivariateTaylorSeriesCategory(Coef,Symbol) with
    coefficient: (%,NNI) -> SMP
      ++\spad{coefficient(s, n)} gives the terms of total degree n.
    coerce: Symbol -> %
      ++\spad{coerce(s)} converts a variable to a Taylor series
    coerce: SMP -> %
      ++\spad{coerce(s)} regroups terms of s by total degree
      ++ and forms a series.
 
    if Coef has Algebra Fraction Integer then
      integrate: (%,Symbol,Coef) -> %
        ++\spad{integrate(s,v,c)} is the integral of s with respect
        ++ to v and having c as the constant of integration.
      fintegrate: (() -> %,Symbol,Coef) -> %
        ++\spad{fintegrate(f,v,c)} is the integral of \spad{f()} with respect
        ++ to v and having c as the constant of integration.
        ++ The evaluation of \spad{f()} is delayed.
 
  Implementation ==> SparseMultivariateTaylorSeries(Coef,Symbol,SMP) add
    Rep := StS -- Below we use the fact that Rep of PS is Stream SMP.
 
    polynomial(s,n) ==
      sum : SMP := 0
      for i in 0..n while not empty? s repeat
        sum := sum + frst s
        s:= rst s
      sum

