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

)abbrev domain ULSCONS UnivariateLaurentSeriesConstructor
++ Authors: Bill Burge, Clifton J. Williamson
++ Date Created: August 1988
++ Date Last Updated: 17 June 1996
++ Description:
++ This package enables one to construct a univariate Laurent series
++ domain from a univariate Taylor series domain. Univariate
++ Laurent series are represented by a pair \spad{[n,f(x)]}, where n is
++ an arbitrary integer and \spad{f(x)} is a Taylor series.  This pair
++ represents the Laurent series \spad{x**n * f(x)}.

UnivariateLaurentSeriesConstructor(Coef,UTS) : SIG == CODE where
  Coef : Ring
  UTS : UnivariateTaylorSeriesCategory Coef

  I     ==> Integer
  L     ==> List
  NNI   ==> NonNegativeInteger
  OUT   ==> OutputForm
  P     ==> Polynomial Coef
  RF    ==> Fraction Polynomial Coef
  RN    ==> Fraction Integer
  ST    ==> Stream Coef
  TERM  ==> Record(k:I,c:Coef)
  monom ==> monomial$UTS
  EFULS ==> ElementaryFunctionsUnivariateLaurentSeries(Coef,UTS,%)
  STTAYLOR ==> StreamTaylorSeriesOperations Coef

  SIG ==> UnivariateLaurentSeriesConstructorCategory(Coef,UTS)

  CODE ==> add

    Rep := Record(expon:I,ps:UTS)

    getExpon : % -> I
    getUTS   : % -> UTS

    getExpon x == x.expon

    getUTS   x == x.ps

    laurent(n,psr) == [n,psr]

    taylorRep x    == getUTS x

    degree x       == getExpon x

    0 == laurent(0,0)

    1 == laurent(0,1)

    monomial(s,e) == laurent(e,s::UTS)

    coerce(uts:UTS):% == laurent(0,uts)

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

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

    taylorIfCan uls ==
      n := getExpon uls
      n < 0 =>
        uls := removeZeroes(-n,uls)
        getExpon(uls) < 0 => "failed"
        getUTS uls
      n = 0 => getUTS uls
      getUTS(uls) * monom(1,n :: NNI)

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

    termExpon: TERM -> I
    termExpon term == term.k

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

    rec: (I,Coef) -> TERM
    rec(exponent,coef) == [exponent,coef]

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

    terms x == recs(coefficients getUTS x,getExpon x)

    recsToCoefs: (Stream TERM,I) -> ST
    recsToCoefs(st,n) == delay
      empty? st => empty()
      term := frst st; ex := termExpon term
      n = ex => concat(termCoef term,recsToCoefs(rst st,n + 1))
      concat(0,recsToCoefs(rst st,n + 1))

    series st ==
      empty? st => 0
      ex := termExpon frst st
      laurent(ex,series recsToCoefs(st,ex))

--% normalizations

    removeZeroes x ==
      empty? coefficients(xUTS := getUTS x) => 0
      coefficient(xUTS,0) = 0 =>
        removeZeroes laurent(getExpon(x) + 1,quoByVar xUTS)
      x

    removeZeroes(n,x) ==
      n <= 0 => x
      empty? coefficients(xUTS := getUTS x) => 0
      coefficient(xUTS,0) = 0 =>
        removeZeroes(n - 1,laurent(getExpon(x) + 1,quoByVar xUTS))
      x

--% predicates

    x = y ==
      EQ(x,y)$Lisp => true
      (expDiff := getExpon(x) - getExpon(y)) = 0 =>
        getUTS(x) = getUTS(y)
      abs(expDiff) > _$streamCount$Lisp => false
      expDiff > 0 =>
        getUTS(x) * monom(1,expDiff :: NNI) = getUTS(y)
      getUTS(y) * monom(1,(- expDiff) :: NNI) = getUTS(x)

    pole? x ==
      (n := degree x) >= 0 => false
      x := removeZeroes(-n,x)
      degree x < 0

--% arithmetic

    x + y  ==
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        laurent(getExpon y,getUTS(y) + getUTS(x) * monom(1,n::NNI))
      laurent(getExpon x,getUTS(x) + getUTS(y) * monom(1,(-n)::NNI))

    x - y  ==
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        laurent(getExpon y,getUTS(x) * monom(1,n::NNI) - getUTS(y))
      laurent(getExpon x,getUTS(x) - getUTS(y) * monom(1,(-n)::NNI))

    x:% * y:% == laurent(getExpon x + getExpon y,getUTS x * getUTS y)

    x:% ** n:NNI ==
      zero? n =>
        zero? x => error "0 ** 0 is undefined"
        1
      laurent(n * getExpon(x),getUTS(x) ** n)

    recip x ==
      x := removeZeroes(1000,x)
      zero? coefficient(x,d := degree x) => "failed"
      (uts := recip getUTS x) case "failed" => "failed"
      laurent(-d,uts :: UTS)

    elt(uls1:%,uls2:%) ==
      (uts := taylorIfCan uls2) case "failed" =>
        error "elt: second argument must have positive order"
      uts2 := uts :: UTS
      not zero? coefficient(uts2,0) =>
        error "elt: second argument must have positive order"
      if (deg := getExpon uls1) < 0 then uls1 := removeZeroes(-deg,uls1)
      (deg := getExpon uls1) < 0 =>
        (recipr := recip(uts2 :: %)) case "failed" =>
          error "elt: second argument not invertible"
        uts1 := taylor(uls1 * monomial(1,-deg))
        (elt(uts1,uts2) :: %) * (recipr :: %) ** ((-deg) :: NNI)
      elt(taylor uls1,uts2) :: %

    eval(uls:%,r:Coef) ==
      if (n := getExpon uls) < 0 then uls := removeZeroes(-n,uls)
      uts := getUTS uls
      (n := getExpon uls) < 0 =>
        zero? r => error "eval: 0 raised to negative power"
        (recipr := recip r) case "failed" =>
          error "eval: non-unit raised to negative power"
        (recipr :: Coef) ** ((-n) :: NNI) *$STTAYLOR eval(uts,r)
      zero? n => eval(uts,r)
      r ** (n :: NNI) *$STTAYLOR eval(uts,r)

--% values

    variable x == variable getUTS x

    center   x == center   getUTS x

    coefficient(x,n) ==
      a := n - getExpon(x)
      a >= 0 => coefficient(getUTS x,a :: NNI)
      0

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

--% other functions

    order x == getExpon x + order getUTS x

    order(x,n) ==
      (m := n - (e := getExpon x)) < 0 => n
      e + order(getUTS x,m :: NNI)

    truncate(x,n) ==
      (m := n - (e := getExpon x)) < 0 => 0
      laurent(e,truncate(getUTS x,m :: NNI))

    truncate(x,n1,n2) ==
      if n2 < n1 then (n1,n2) := (n2,n1)
      (m1 := n1 - (e := getExpon x)) < 0 => truncate(x,n2)
      laurent(e,truncate(getUTS x,m1 :: NNI,(n2 - e) :: NNI))

    if Coef has IntegralDomain then

      rationalFunction(x,n) ==
        (m := n - (e := getExpon x)) < 0 => 0
        poly := polynomial(getUTS x,m :: NNI) :: RF
        zero? e => poly
        v := variable(x) :: RF; c := center(x) :: P :: RF
        positive? e => poly * (v - c) ** (e :: NNI)
        poly / (v - c) ** ((-e) :: NNI)

      rationalFunction(x,n1,n2) ==
        if n2 < n1 then (n1,n2) := (n2,n1)
        (m1 := n1 - (e := getExpon x)) < 0 => rationalFunction(x,n2)
        poly := polynomial(getUTS x,m1 :: NNI,(n2 - e) :: NNI) :: RF
        zero? e => poly
        v := variable(x) :: RF; c := center(x) :: P :: RF
        positive? e => poly * (v - c) ** (e :: NNI)
        poly / (v - c) ** ((-e) :: NNI)

      x exquo y ==
        x := removeZeroes(1000,x)
        y := removeZeroes(1000,y)
        zero? coefficient(y, d := degree y) => "failed"
        (uts := (getUTS x) exquo (getUTS y)) case "failed" => "failed"
        laurent(degree x-d,uts :: UTS)

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

        approximate(x,n) ==
          (m := n - (e := getExpon x)) < 0 => 0
          app := approximate(getUTS x,m :: NNI)
          zero? e => app
          app * ((variable(x) :: Coef) - center(x)) ** e

    complete x == laurent(getExpon x,complete getUTS x)

    extend(x,n) ==
      e := getExpon x
      (m := n - e) < 0 => x
      laurent(e,extend(getUTS x,m :: NNI))

    map(f:Coef -> Coef,x:%) == laurent(getExpon x,map(f,getUTS x))

    multiplyCoefficients(f,x) ==
      e := getExpon x
      laurent(e,multiplyCoefficients((z1:I):Coef +-> f(e + z1),getUTS x))

    multiplyExponents(x,n) ==
      laurent(n * getExpon x,multiplyExponents(getUTS x,n))

    differentiate x ==
      e := getExpon x
      laurent(e - 1,
              multiplyCoefficients((z1:I):Coef +-> (e + z1)::Coef,getUTS x))

    if Coef has PartialDifferentialRing(Symbol) then

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

    characteristic() == characteristic()$Coef

    if Coef has Field then

      retract(x:%):UTS                      == taylor x

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

      (x:%) ** (n:I) ==
        zero? n =>
          zero? x => error "0 ** 0 is undefined"
          1
        n > 0 => laurent(n * getExpon(x),getUTS(x) ** (n :: NNI))
        xInv := inv x; minusN := (-n) :: NNI
        laurent(minusN * getExpon(xInv),getUTS(xInv) ** minusN)

      (x:UTS) * (y:%) == (x :: %) * y

      (x:%) * (y:UTS) == x * (y :: %)

      inv x ==
        (xInv := recip x) case "failed" =>
          error "multiplicative inverse does not exist"
        xInv :: %

      (x:%) / (y:%) ==
        (yInv := recip y) case "failed" =>
          error "inv: multiplicative inverse does not exist"
        x * (yInv :: %)

      (x:UTS) / (y:UTS) == (x :: %) / (y :: %)

      numer x ==
        (n := degree x) >= 0 => taylor x
        x := removeZeroes(-n,x)
        (n := degree x) = 0 => taylor x
        getUTS x

      denom x ==
        (n := degree x) >= 0 => 1
        x := removeZeroes(-n,x)
        (n := degree x) = 0 => 1
        monom(1,(-n) :: NNI)

--% algebraic and transcendental functions

    if Coef has Algebra Fraction Integer then

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

      if Coef has Field then
         (x:%) ** (r:RN) == x **$EFULS r

      exp x   == exp(x)$EFULS

      log x   == log(x)$EFULS

      sin x   == sin(x)$EFULS

      cos x   == cos(x)$EFULS

      tan x   == tan(x)$EFULS

      cot x   == cot(x)$EFULS

      sec x   == sec(x)$EFULS

      csc x   == csc(x)$EFULS

      asin x  == asin(x)$EFULS

      acos x  == acos(x)$EFULS

      atan x  == atan(x)$EFULS

      acot x  == acot(x)$EFULS

      asec x  == asec(x)$EFULS

      acsc x  == acsc(x)$EFULS

      sinh x  == sinh(x)$EFULS

      cosh x  == cosh(x)$EFULS

      tanh x  == tanh(x)$EFULS

      coth x  == coth(x)$EFULS

      sech x  == sech(x)$EFULS

      csch x  == csch(x)$EFULS

      asinh x == asinh(x)$EFULS

      acosh x == acosh(x)$EFULS

      atanh x == atanh(x)$EFULS

      acoth x == acoth(x)$EFULS

      asech x == asech(x)$EFULS

      acsch x == acsch(x)$EFULS

      ratInv: I -> Coef
      ratInv n ==
        zero? n => 1
        inv(n :: RN) :: Coef

      integrate x ==
        not zero? coefficient(x,-1) =>
          error "integrate: series has term of order -1"
        e := getExpon x
        laurent(e+1,multiplyCoefficients((z:I):Coef+->ratInv(e+1+z),getUTS 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((z1:Coef):Coef+->integrate(z1,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((z1:Coef):Coef +-> integrateWithOneAnswer(z1,s),x)
          error "integrate: center is a function of variable of integration"

    termOutput:(I,Coef,OUT) -> OUT
    termOutput(k,c,vv) ==
    -- creates a term c * vv ** k
      k = 0 => c :: OUT
      mon :=
        k = 1 => vv
        vv ** (k :: OUT)
      c = 1 => mon
      c = -1 => -mon
      (c :: OUT) * mon

    -- check a global Lisp variable
    showAll?:() -> Boolean
    showAll?() == true

    termsToOutputForm:(I,ST,OUT) -> OUT
    termsToOutputForm(m,uu,xxx) ==
      l : L OUT := empty()
      empty? uu => (0$Coef) :: OUT
      n : NNI ; count : NNI := _$streamCount$Lisp
      for n in 0..count while not empty? uu repeat
        if frst(uu) ^= 0 then
          l := concat(termOutput((n :: I) + m,frst(uu),xxx),l)
        uu := rst uu
      if showAll?() then
        for n in (count + 1).. while explicitEntries? uu and _
               not eq?(uu,rst uu) repeat
          if frst(uu) ^= 0 then
            l := concat(termOutput((n::I) + m,frst(uu),xxx),l)
          uu := rst uu
      l :=
        explicitlyEmpty? uu => l
        eq?(uu,rst uu) and frst uu = 0 => l
        concat(prefix("O" :: OUT,[xxx ** ((n :: I) + m) :: OUT]),l)
      empty? l => (0$Coef) :: OUT
      reduce("+",reverse_! l)

    coerce(x:%):OUT ==
      x := removeZeroes(_$streamCount$Lisp,x)
      m := degree x
      uts := getUTS x
      p := coefficients uts
      var := variable uts; cen := center uts
      xxx :=
        zero? cen => var :: OUT
        paren(var :: OUT - cen :: OUT)
      termsToOutputForm(m,p,xxx)