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

)abbrev domain SULS SparseUnivariateLaurentSeries
++ Author: Clifton J. Williamson
++ Date Created: 11 November 1994
++ Date Last Updated: 10 March 1995
++ Description:
++ Sparse Laurent series in one variable
++ \spadtype{SparseUnivariateLaurentSeries} is a domain representing Laurent
++ 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,
++ \spad{SparseUnivariateLaurentSeries(Integer,x,3)} represents Laurent
++ series in \spad{(x - 3)} with integer coefficients.

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

  I     ==> Integer
  NNI   ==> NonNegativeInteger
  OUT   ==> OutputForm
  P     ==> Polynomial Coef
  RF    ==> Fraction Polynomial Coef
  RN    ==> Fraction Integer
  S     ==> String
  SUTS  ==> SparseUnivariateTaylorSeries(Coef,var,cen)
  EFULS ==> ElementaryFunctionsUnivariateLaurentSeries(Coef,SUTS,%)

  SIG ==> UnivariateLaurentSeriesConstructorCategory(Coef,SUTS) with

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

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

    if Coef has Algebra Fraction Integer then

      integrate : (%,Variable(var)) -> %
        ++ \spad{integrate(f(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 ==> InnerSparseUnivariatePowerSeries(Coef) add

    Rep := InnerSparseUnivariatePowerSeries(Coef)

    variable x == var
    center   x == cen

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

    pole? x == negative? order(x,0)

--% operations with Taylor series

    coerce(uts:SUTS) == uts pretend %

    taylorIfCan uls ==
      pole? uls => "failed"
      uls pretend SUTS

    taylor uls ==
      (uts := taylorIfCan uls) case "failed" =>
        error "taylor: Laurent series has a pole"
      uts :: SUTS

    retractIfCan(x:%):Union(SUTS,"failed") == taylorIfCan x

    laurent(n,uts) == monomial(1,n) * (uts :: %)

    removeZeroes uls    == uls
    removeZeroes(n,uls) == uls

    taylorRep uls == taylor(monomial(1,-order(uls,0)) * uls)
    degree uls    == order(uls,0)

    numer uls == taylorRep uls
    denom uls == monomial(1,(-order(uls,0)) :: NNI)$SUTS

    (uts:SUTS) * (uls:%) == (uts :: %) * uls
    (uls:%) * (uts:SUTS) == uls * (uts :: %)

    if Coef has Field then
      (uts1:SUTS) / (uts2:SUTS) == (uts1 :: %) / (uts2 :: %)

    recip(uls) == iExquo(1,uls,false)

    if Coef has IntegralDomain then
      uls1 exquo uls2 == iExquo(uls1,uls2,false)

    if Coef has Field then
      uls1:% / uls2:% ==
        (q := uls1 exquo uls2) case "failed" =>
          error "quotient cannot be computed"
        q :: %

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

    elt(uls1:%,uls2:%) ==
      order(uls2,1) < 1 =>
        error "elt: second argument must have positive order"
      negative?(ord := order(uls1,0)) =>
        (recipr := recip uls2) case "failed" =>
          error "elt: second argument not invertible"
        uls3 := uls1 * monomial(1,-ord)
        iCompose(uls3,uls2) * (recipr :: %) ** ((-ord) :: NNI)
      iCompose(uls1,uls2)

    if Coef has IntegralDomain then
      rationalFunction(uls,n) ==
        zero?(e := order(uls,0)) =>
          negative? n => 0
          polynomial(taylor uls,n :: NNI) :: RF
        negative?(m := n - e) => 0
        poly := polynomial(taylor(monomial(1,-e) * uls),m :: NNI) :: RF
        v := variable(uls) :: RF; c := center(uls) :: P :: RF
        poly / (v - c) ** ((-e) :: NNI)

      rationalFunction(uls,n1,n2) == rationalFunction(truncate(uls,n1,n2),n2)

    if Coef has Algebra Fraction Integer then

      integrate uls ==
        zero? coefficient(uls,-1) =>
          error "integrate: series has term of order -1"
        integrate(uls)$Rep

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

      (uls1:%) ** (uls2:%) == exp(log(uls1) * uls2)

      exp uls   == exp(uls)$EFULS

      log uls   == log(uls)$EFULS

      sin uls   == sin(uls)$EFULS

      cos uls   == cos(uls)$EFULS

      tan uls   == tan(uls)$EFULS

      cot uls   == cot(uls)$EFULS

      sec uls   == sec(uls)$EFULS

      csc uls   == csc(uls)$EFULS

      asin uls  == asin(uls)$EFULS

      acos uls  == acos(uls)$EFULS

      atan uls  == atan(uls)$EFULS

      acot uls  == acot(uls)$EFULS

      asec uls  == asec(uls)$EFULS

      acsc uls  == acsc(uls)$EFULS

      sinh uls  == sinh(uls)$EFULS

      cosh uls  == cosh(uls)$EFULS

      tanh uls  == tanh(uls)$EFULS

      coth uls  == coth(uls)$EFULS

      sech uls  == sech(uls)$EFULS

      csch uls  == csch(uls)$EFULS

      asinh uls == asinh(uls)$EFULS

      acosh uls == acosh(uls)$EFULS

      atanh uls == atanh(uls)$EFULS

      acoth uls == acoth(uls)$EFULS

      asech uls == asech(uls)$EFULS

      acsch uls == acsch(uls)$EFULS

      if Coef has CommutativeRing then

        (uls:%) ** (r:RN) == cRationalPower(uls,r)

      else

        (uls:%) ** (r:RN) ==
          negative?(ord0 := order(uls,0)) =>
            order := ord0 :: I
            (n := order exquo denom(r)) case "failed" =>
              error "**: rational power does not exist"
            uts := retract(uls * monomial(1,-order))@SUTS
            utsPow := (uts ** r) :: %
            monomial(1,(n :: I) * numer(r)) * utsPow
          uts := retract(uls)@SUTS
          (uts ** r) :: %

--% OutputForms

    coerce(uls:%): OUT ==
      st := getStream uls
      if not(explicitlyEmpty? st or explicitEntries? st) _
        and (nx := retractIfCan(elt getRef uls))@Union(I,"failed") case I then
        count : NNI := _$streamCount$Lisp
        degr := min(count,(nx :: I) + count + 1)
        extend(uls,degr)
      seriesToOutputForm(st,getRef uls,variable uls,center uls,1)