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

)abbrev domain UPXSCONS UnivariatePuiseuxSeriesConstructor
++ Author: Clifton J. Williamson
++ Date Created: 9 May 1989
++ Date Last Updated: 30 November 1994
++ Description:
++ This package enables one to construct a univariate Puiseux series
++ domain from a univariate Laurent series domain. Univariate
++ Puiseux series are represented by a pair \spad{[r,f(x)]}, where r is
++ a positive rational number and \spad{f(x)} is a Laurent series.
++ This pair represents the Puiseux series \spad{f(x^r)}.

UnivariatePuiseuxSeriesConstructor(Coef,ULS) : SIG == CODE where
  Coef : Ring
  ULS : UnivariateLaurentSeriesCategory Coef

  I     ==> Integer
  L     ==> List
  NNI   ==> NonNegativeInteger
  OUT   ==> OutputForm
  PI    ==> PositiveInteger
  RN    ==> Fraction Integer
  ST    ==> Stream Coef
  LTerm ==> Record(k:I,c:Coef)
  PTerm ==> Record(k:RN,c:Coef)
  ST2LP ==> StreamFunctions2(LTerm,PTerm)
  ST2PL ==> StreamFunctions2(PTerm,LTerm)

  SIG ==> UnivariatePuiseuxSeriesConstructorCategory(Coef,ULS)

  CODE ==> add

    Rep := Record(expon:RN,lSeries:ULS)

    getExpon: % -> RN
    getULS  : % -> ULS

    getExpon pxs == pxs.expon

    getULS   pxs == pxs.lSeries

    puiseux(n,ls)   == [n,ls]

    laurentRep x    == getULS x

    rationalPower x == getExpon x

    degree x        == getExpon(x) * degree(getULS(x))

    0 == puiseux(1,0)

    1 == puiseux(1,1)

    monomial(c,k) ==
      k = 0 => c :: %
      k < 0 => puiseux(-k,monomial(c,-1))
      puiseux(k,monomial(c,1))

    coerce(ls:ULS) == puiseux(1,ls)

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

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

    laurentIfCan upxs ==
      r := getExpon upxs
      (denom r) = 1 =>
        multiplyExponents(getULS upxs,numer(r) :: PI)
      "failed"

    laurent upxs ==
      (uls := laurentIfCan upxs) case "failed" =>
        error "laurent: Puiseux series has fractional powers"
      uls :: ULS

    multExp: (RN,LTerm) -> PTerm
    multExp(r,lTerm) == [r * lTerm.k,lTerm.c]

    terms upxs == 
      map((t1:LTerm):PTerm+->multExp(getExpon upxs,t1),terms getULS upxs)$ST2LP

    clearDen: (I,PTerm) -> LTerm
    clearDen(n,lTerm) ==
      (int := retractIfCan(n * lTerm.k)@Union(I,"failed")) case "failed" =>
        error "series: inappropriate denominator"
      [int :: I,lTerm.c]

    series(n,stream) ==
      str := map((t1:PTerm):LTerm +-> clearDen(n,t1),stream)$ST2PL
      puiseux(1/n,series str)

--% normalizations

    rewrite:(%,PI) -> %
    rewrite(upxs,m) ==
      -- rewrites a series in x**r as a series in x**(r/m)
      puiseux((getExpon upxs)*(1/m),multiplyExponents(getULS upxs,m))

    ratGcd: (RN,RN) -> RN
    ratGcd(r1,r2) ==
      -- if r1 = prod(p prime,p ** ep(r1)) and
      -- if r2 = prod(p prime,p ** ep(r2)), then
      -- ratGcd(r1,r2) = prod(p prime,p ** min(ep(r1),ep(r2)))
      gcd(numer r1,numer r2) / lcm(denom r1,denom r2)

    withNewExpon:(%,RN) -> %
    withNewExpon(upxs,r) ==
      rewrite(upxs,numer(getExpon(upxs)/r) pretend PI)

--% predicates

    upxs1 = upxs2 ==
      r1 := getExpon upxs1; r2 := getExpon upxs2
      ls1 := getULS upxs1; ls2 := getULS upxs2
      (r1 = r2) => (ls1 = ls2)
      r := ratGcd(r1,r2)
      m1 := numer(getExpon(upxs1)/r) pretend PI
      m2 := numer(getExpon(upxs2)/r) pretend PI
      multiplyExponents(ls1,m1) = multiplyExponents(ls2,m2)

    pole? upxs == pole? getULS upxs

--% arithmetic

    applyFcn:((ULS,ULS) -> ULS,%,%) -> %
    applyFcn(op,pxs1,pxs2) ==
      r1 := getExpon pxs1; r2 := getExpon pxs2
      ls1 := getULS pxs1; ls2 := getULS pxs2
      (r1 = r2) => puiseux(r1,op(ls1,ls2))
      r := ratGcd(r1,r2)
      m1 := numer(getExpon(pxs1)/r) pretend PI
      m2 := numer(getExpon(pxs2)/r) pretend PI
      puiseux(r,op(multiplyExponents(ls1,m1),multiplyExponents(ls2,m2)))

    pxs1 + pxs2     == applyFcn((z1:ULS,z2:ULS):ULS+->z1 +$ULS z2,pxs1,pxs2)

    pxs1 - pxs2     == applyFcn((z1:ULS,z2:ULS):ULS+->z1 -$ULS z2,pxs1,pxs2)

    pxs1:% * pxs2:% == applyFcn((z1:ULS,z2:ULS):ULS+->z1 *$ULS z2,pxs1,pxs2)

    pxs:% ** n:NNI == puiseux(getExpon pxs,getULS(pxs)**n)

    recip pxs ==
      rec := recip getULS pxs
      rec case "failed" => "failed"
      puiseux(getExpon pxs,rec :: ULS)

    RATALG : Boolean := Coef has Algebra(Fraction Integer)

    elt(upxs1:%,upxs2:%) ==
      uls1 := laurentRep upxs1; uls2 := laurentRep upxs2
      r1 := rationalPower upxs1; r2 := rationalPower upxs2
      (n := retractIfCan(r1)@Union(Integer,"failed")) case Integer =>
        puiseux(r2,uls1(uls2 ** r1))
      RATALG =>
        if zero? (coef := coefficient(uls2,deg := degree uls2)) then
          deg := order(uls2,deg + 1000)
          zero? (coef := coefficient(uls2,deg)) =>
            error "elt: series with many leading zero coefficients"
        -- a fractional power of a Laurent series may not be defined:
        -- if f(x) = c * x**n + ..., then f(x) ** (p/q) will be defined
        -- only if q divides n
        b := lcm(denom r1,deg); c := b quo deg
        mon : ULS := monomial(1,c)
        uls2 := elt(uls2,mon) ** r1
        puiseux(r2*(1/c),elt(uls1,uls2))
      error "elt: rational powers not available for this coefficient domain"

    if Coef has "**": (Coef,Integer) -> Coef and
       Coef has "**": (Coef, RN) -> Coef then

         eval(upxs:%,a:Coef) == eval(getULS upxs,a ** getExpon(upxs))

    if Coef has Field then

      pxs1:% / pxs2:% == applyFcn((z1:ULS,z2:ULS):ULS+->z1 /$ULS z2,pxs1,pxs2)

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

--% values

    variable upxs == variable getULS upxs

    center   upxs == center   getULS upxs

    coefficient(upxs,rn) ==
      (denom(n := rn / getExpon upxs)) = 1 =>
        coefficient(getULS upxs,numer n)
      0

    elt(upxs:%,rn:RN) == coefficient(upxs,rn)

--% other functions

    roundDown: RN -> I
    roundDown rn ==
      -- returns the largest integer <= rn
      (den := denom rn) = 1 => numer rn
      n := (num := numer rn) quo den
      positive?(num) => n
      n - 1

    roundUp: RN -> I
    roundUp rn ==
      -- returns the smallest integer >= rn
      (den := denom rn) = 1 => numer rn
      n := (num := numer rn) quo den
      positive?(num) => n + 1
      n

    order upxs == getExpon upxs * order getULS upxs

    order(upxs,r) ==
      e := getExpon upxs
      ord := order(getULS upxs, n := roundDown(r / e))
      ord = n => r
      ord * e

    truncate(upxs,r) ==
      e := getExpon upxs
      puiseux(e,truncate(getULS upxs,roundDown(r / e)))

    truncate(upxs,r1,r2) ==
      e := getExpon upxs
      puiseux(e,truncate(getULS upxs,roundUp(r1 / e),roundDown(r2 / e)))

    complete upxs == puiseux(getExpon upxs,complete getULS upxs)

    extend(upxs,r) ==
      e := getExpon upxs
      puiseux(e,extend(getULS upxs,roundDown(r / e)))

    map(fcn,upxs) == puiseux(getExpon upxs,map(fcn,getULS upxs))

    characteristic() == characteristic()$Coef

    multiplyExponents(upxs:%,n:RN) ==
      puiseux(n * getExpon(upxs),getULS upxs)

    multiplyExponents(upxs:%,n:PI) ==
      puiseux(n * getExpon(upxs),getULS upxs)

    if Coef has "*": (Fraction Integer, Coef) -> Coef then

      differentiate upxs ==
        r := getExpon upxs
        puiseux(r,differentiate getULS upxs) * monomial(r :: Coef,r-1)

      if Coef has PartialDifferentialRing(Symbol) then

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

    if Coef has Algebra Fraction Integer then

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

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

      integrate upxs ==
        not zero? coefficient(upxs,-1) =>
          error "integrate: series has term of order -1"
        r := getExpon upxs
        uls := getULS upxs
        uls := multiplyCoefficients((z1:Integer):Coef+->ratInv(z1*r+1),uls)
        monomial(1,1) * puiseux(r,uls)

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

        integrate(upxs:%,s:Symbol) ==
          (s = variable(upxs)) => integrate upxs
          not entry?(s,variables center upxs) 
            => map((z1:Coef):Coef +-> integrate(z1,s),upxs)
          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(upxs:%,s:Symbol) ==
          (s = variable(upxs)) => integrate upxs
          not entry?(s,variables center upxs) =>
            map((z1:Coef):Coef +-> integrateWithOneAnswer(z1,s),upxs)
          error "integrate: center is a function of variable of integration"

      if Coef has Field then

         (upxs:%) ** (q:RN) ==
           num := numer q; den := denom q
           den = 1 => upxs ** num
           r := rationalPower upxs; uls := laurentRep upxs
           deg := degree uls
           if zero?(coef := coefficient(uls,deg)) then
             deg := order(uls,deg + 1000)
             zero?(coef := coefficient(uls,deg)) =>
               error "power of series with many leading zero coefficients"
           ulsPow := (uls * monomial(1,-deg)$ULS) ** q
           puiseux(r,ulsPow) * monomial(1,deg*q*r)

      applyUnary: (ULS -> ULS,%) -> %
      applyUnary(fcn,upxs) ==
        puiseux(rationalPower upxs,fcn laurentRep upxs)

      exp upxs   == applyUnary(exp,upxs)

      log upxs   == applyUnary(log,upxs)

      sin upxs   == applyUnary(sin,upxs)

      cos upxs   == applyUnary(cos,upxs)

      tan upxs   == applyUnary(tan,upxs)

      cot upxs   == applyUnary(cot,upxs)

      sec upxs   == applyUnary(sec,upxs)

      csc upxs   == applyUnary(csc,upxs)

      asin upxs  == applyUnary(asin,upxs)

      acos upxs  == applyUnary(acos,upxs)

      atan upxs  == applyUnary(atan,upxs)

      acot upxs  == applyUnary(acot,upxs)

      asec upxs  == applyUnary(asec,upxs)

      acsc upxs  == applyUnary(acsc,upxs)

      sinh upxs  == applyUnary(sinh,upxs)

      cosh upxs  == applyUnary(cosh,upxs)

      tanh upxs  == applyUnary(tanh,upxs)

      coth upxs  == applyUnary(coth,upxs)

      sech upxs  == applyUnary(sech,upxs)

      csch upxs  == applyUnary(csch,upxs)

      asinh upxs == applyUnary(asinh,upxs)

      acosh upxs == applyUnary(acosh,upxs)

      atanh upxs == applyUnary(atanh,upxs)

      acoth upxs == applyUnary(acoth,upxs)

      asech upxs == applyUnary(asech,upxs)

      acsch upxs == applyUnary(acsch,upxs)