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

)abbrev domain PFR PartialFraction
++ Author: Robert S. Sutor
++ Date Created: 1986
++ Change History: 05/20/91 BMT Converted to the new library
++ Description:
++ The domain \spadtype{PartialFraction} implements partial fractions
++ over a euclidean domain \spad{R}.  This requirement on the
++ argument domain allows us to normalize the fractions.  Of
++ particular interest are the 2 forms for these fractions.  The
++ ``compact'' form has only one fractional term per prime in the
++ denominator, while the ``p-adic'' form expands each numerator
++ p-adically via the prime p in the denominator.  For computational
++ efficiency, the compact form is used, though the p-adic form may
++ be gotten by calling the function padicFraction}.  For a
++ general euclidean domain, it is not known how to factor the
++ denominator.  Thus the function partialFraction takes as its
++ second argument an element of \spadtype{Factored(R)}.

PartialFraction(R) : SIG == CODE where
  R : EuclideanDomain

  FRR  ==> Factored R
  SUPR ==> SparseUnivariatePolynomial R

  SIG ==> Join(Field, Algebra R) with

    coerce : % -> Fraction R
      ++ coerce(p) sums up the components of the partial fraction and
      ++ returns a single fraction.
      ++
      ++X a:=(13/74)::PFR(INT)
      ++X a::FRAC(INT)

    coerce : Fraction FRR -> %
      ++ coerce(f) takes a fraction with numerator and denominator in
      ++ factored form and creates a partial fraction.  It is
      ++ necessary for the parts to be factored because it is not
      ++ known in general how to factor elements of \spad{R} and
      ++ this is needed to decompose into partial fractions.
      ++
      ++X (13/74)::PFR(INT)

    compactFraction : % -> %
      ++ compactFraction(p) normalizes the partial fraction \spad{p}
      ++ to the compact representation. In this form, the partial
      ++ fraction has only one fractional term per prime in the
      ++ denominator.
      ++
      ++X a:=partialFraction(1,factorial 10)
      ++X b:=padicFraction(a)
      ++X compactFraction(b)

    firstDenom : % -> FRR
      ++ firstDenom(p) extracts the denominator of the first fractional
      ++ term. This returns 1 if there is no fractional part (use
      ++ wholePart from PartialFraction to get the whole part).
      ++
      ++X a:=partialFraction(1,factorial 10)
      ++X firstDenom(a)

    firstNumer : % -> R
      ++ firstNumer(p) extracts the numerator of the first fractional
      ++ term. This returns 0 if there is no fractional part (use
      ++ wholePart from PartialFraction to get the whole part).
      ++
      ++X a:=partialFraction(1,factorial 10)
      ++X firstNumer(a)

    nthFractionalTerm : (%,Integer) -> %
      ++ nthFractionalTerm(p,n) extracts the nth fractional term from
      ++ the partial fraction \spad{p}.  This returns 0 if the index
      ++ \spad{n} is out of range.
      ++
      ++X a:=partialFraction(1,factorial 10)
      ++X b:=padicFraction(a)
      ++X nthFractionalTerm(b,3)

    numberOfFractionalTerms : % -> Integer
      ++ numberOfFractionalTerms(p) computes the number of fractional
      ++ terms in \spad{p}. This returns 0 if there is no fractional
      ++ part.
      ++
      ++X a:=partialFraction(1,factorial 10)
      ++X b:=padicFraction(a)
      ++X numberOfFractionalTerms(b)

    padicallyExpand : (R,R) -> SUPR
      ++ padicallyExpand(p,x) is a utility function that expands
      ++ the second argument \spad{x} ``p-adically'' in
      ++ the first.

    padicFraction : % -> %
      ++ padicFraction(q) expands the fraction p-adically in the primes
      ++ \spad{p} in the denominator of \spad{q}. For example,
      ++ \spad{padicFraction(3/(2**2)) = 1/2 + 1/(2**2)}.
      ++ Use compactFraction from PartialFraction to 
      ++ return to compact form.
      ++
      ++X a:=partialFraction(1,factorial 10)
      ++X padicFraction(a)

    partialFraction : (R, FRR) -> %
      ++ partialFraction(numer,denom) is the main function for
      ++ constructing partial fractions. The second argument is the
      ++ denominator and should be factored.
      ++
      ++X partialFraction(1,factorial 10)

    wholePart : % -> R
      ++ wholePart(p) extracts the whole part of the partial fraction
      ++ \spad{p}.
      ++
      ++X a:=(74/13)::PFR(INT)
      ++X wholePart(a)

  CODE ==> add

    -- some constructor assignments and macros

    Ex     ==> OutputForm

    fTerm  ==> Record(num: R, den: FRR)           -- den should have
                                                  -- unit = 1 and only
                                                  -- 1 factor
    LfTerm ==> List Record(num: R, den: FRR)

    QR     ==> Record(quotient: R, remainder: R)

    Rep    := Record(whole:R, fract: LfTerm)

    -- private function signatures

    copypf: % -> %

    LessThan: (fTerm, fTerm) -> Boolean

    multiplyFracTerms: (fTerm, fTerm) -> %

    normalizeFracTerm: fTerm -> %

    partialFractionNormalized: (R, FRR) -> %

    -- declarations

    a,b: %

    n: Integer

    r: R

    -- private function definitions

    copypf(a: %): % == [a.whole,copy a.fract]$%

    LessThan(s: fTerm, t: fTerm) ==
      -- have to wait until FR has < operation
      if (GGREATERP(s.den,t.den)$Lisp : Boolean) then false
      else true

    multiplyFracTerms(s : fTerm, t : fTerm) ==
      nthFactor(s.den,1) = nthFactor(t.den,1) =>
        normalizeFracTerm([s.num * t.num, s.den * t.den]$fTerm) : Rep
      i : Union(Record(coef1: R, coef2: R),"failed")
      coefs : Record(coef1: R, coef2: R)
      i := extendedEuclidean(expand t.den, expand s.den,s.num * t.num)
      i case "failed" => error "PartialFraction: not in ideal"
      coefs := (i :: Record(coef1: R, coef2: R))
      c : % := copypf 0$%
      d : %
      if coefs.coef2 ^= 0$R then
        c := normalizeFracTerm ([coefs.coef2, t.den]$fTerm)
      if coefs.coef1 ^= 0$R then
        d := normalizeFracTerm ([coefs.coef1, s.den]$fTerm)
        c.whole := c.whole + d.whole
        not (null d.fract) => c.fract := append(d.fract,c.fract)
      c

    normalizeFracTerm(s : fTerm) ==
      -- makes sure num is "less than" den, whole may be non-zero
      qr : QR := divide(s.num, (expand s.den))
      qr.remainder = 0$R => [qr.quotient, nil()$LfTerm]
      -- now verify num and den are coprime
      f : R := nthFactor(s.den,1)
      nexpon : Integer := nthExponent(s.den,1)
      expon  : Integer := 0
      q : QR := divide(qr.remainder, f)
      while (q.remainder = 0$R) and (expon < nexpon) repeat
        expon := expon + 1
        qr.remainder := q.quotient
        q := divide(qr.remainder,f)
      expon = 0 => [qr.quotient,[[qr.remainder, s.den]$fTerm]$LfTerm]
      expon = nexpon => (qr.quotient + qr.remainder) :: %
      [qr.quotient,[[qr.remainder, nilFactor(f,nexpon-expon)]$fTerm]$LfTerm]

    partialFractionNormalized(nm: R, dn : FRR) ==
      -- assume unit dn = 1
      nm = 0$R   => 0$%
      dn = 1$FRR => nm :: %
      qr : QR := divide(nm, expand dn)
      c : % := [0$R,[[qr.remainder,
        nilFactor(nthFactor(dn,1), nthExponent(dn,1))]$fTerm]$LfTerm]
      d : %
      for i in 2..numberOfFactors(dn) repeat
        d :=
          [0$R,[[1$R,nilFactor(nthFactor(dn,i),_
            nthExponent(dn,i))]$fTerm]$LfTerm]
        c := c * d
      (qr.quotient :: %) + c

    -- public function definitions

    padicFraction(a : %) ==
      b: % := compactFraction a
      null b.fract => b
      l : LfTerm := nil
      s : fTerm
      f : R
      e,d: Integer
      for s in b.fract repeat
        e := nthExponent(s.den,1)
        e = 1 => l := cons(s,l)
        f := nthFactor(s.den,1)
        d := degree(sp := padicallyExpand(f,s.num))
        while (sp ^= 0$SUPR) repeat
          l := cons([leadingCoefficient sp,nilFactor(f,e-d)]$fTerm, l)
          d := degree(sp := reductum sp)
      [b.whole, sort(LessThan,l)]$%

    compactFraction(a : %) ==
      -- only one power for each distinct denom will remain
      2 > # a.fract => a
      af : LfTerm := reverse a.fract
      bf : LfTerm := nil
      bw : R := a.whole
      b : %
      s : fTerm := [(first af).num,(first af).den]$fTerm
      f : R := nthFactor(s.den,1)
      e : Integer := nthExponent(s.den,1)
      t : fTerm
      for t in rest af repeat
        f = nthFactor(t.den,1) =>
          s.num := s.num + (t.num *
            (f **$R ((e - nthExponent(t.den,1)) : NonNegativeInteger)))
        b := normalizeFracTerm s
        bw := bw + b.whole
        if not (null b.fract) then bf := cons(first b.fract,bf)
        s := [t.num, t.den]$fTerm
        f := nthFactor(s.den,1)
        e := nthExponent(s.den,1)
      b := normalizeFracTerm s
      [bw + b.whole,append(b.fract,bf)]$%

    0 == [0$R, nil()$LfTerm]

    1 == [1$R, nil()$LfTerm]

    characteristic() == characteristic()$R

    coerce(r): % == [r, nil()$LfTerm]

    coerce(n): % == [(n :: R), nil()$LfTerm]

    coerce(a): Fraction R ==
      q : Fraction R := (a.whole :: Fraction R)
      s : fTerm
      for s in a.fract repeat
        q := q + (s.num / (expand s.den))
      q

    coerce(q: Fraction FRR): % ==
      u : R := (recip unit denom q):: R
      r1 : R := u * expand numer q
      partialFractionNormalized(r1, u * denom q)

    a exquo b ==
      b = 0$% => "failed"
      b = 1$% => a
      br : Fraction R := inv (b :: Fraction R)
      a * partialFraction(numer br,(denom br) :: FRR)

    recip a == (1$% exquo a)

    firstDenom a ==         -- denominator of 1st fractional term
      null a.fract => 1$FRR
      (first a.fract).den

    firstNumer a ==         -- numerator of 1st fractional term
      null a.fract => 0$R
      (first a.fract).num

    numberOfFractionalTerms a == # a.fract

    nthFractionalTerm(a,n) ==
      l : LfTerm := a.fract
      (n < 1) or (n > # l) => 0$%
      [0$R,[l.n]$LfTerm]$%

    wholePart a == a.whole

    partialFraction(nm: R, dn : FRR) ==
      nm = 0$R => 0$%
      -- move inv unit of den to numerator
      u : R := unit dn
      u := (recip u) :: R
      partialFractionNormalized(u * nm,u * dn)

    padicallyExpand(p : R, r : R) ==
      -- expands r as a sum of powers of p, with coefficients
      -- r = HornerEval(padicallyExpand(p,r),p)
      qr : QR := divide(r, p)
      qr.quotient = 0$R => qr.remainder :: SUPR
      (qr.remainder :: SUPR) + monomial(1$R,1$NonNegativeInteger)$SUPR *
        padicallyExpand(p,qr.quotient)

    a = b ==
      a.whole ^= b.whole => false  -- must verify this
      (null a.fract) =>
        null b.fract => a.whole = b.whole
        false
      null b.fract => false
      -- oh, no! following is temporary
      (a :: Fraction R) = (b :: Fraction R)

    - a ==
      s: fTerm
      l: LfTerm := nil
      for s in reverse a.fract repeat l := cons([- s.num,s.den]$fTerm,l)
      [- a.whole,l]

    r * a ==
      r = 0$R => 0$%
      r = 1$R => a
      b : % := (r * a.whole) :: %
      c : %
      s : fTerm
      for s in reverse a.fract repeat
        c := normalizeFracTerm [r * s.num, s.den]$fTerm
        b.whole := b.whole + c.whole
        not (null c.fract) => b.fract := append(c.fract, b.fract)
      b

    n * a == (n :: R) * a

    a + b ==
      compactFraction
        [a.whole + b.whole,
          sort(LessThan,append(a.fract,copy b.fract))]$%

    a * b ==
      null a.fract => a.whole * b
      null b.fract => b.whole * a
      af : % := [0$R, a.fract]$%   --     a - a.whole
      c: % := (a.whole * b) + (b.whole * af)
      s,t : fTerm
      for s in a.fract repeat
        for t in b.fract repeat
          c := c + multiplyFracTerms(s,t)
      c

    coerce(a): Ex ==
      null a.fract => a.whole :: Ex
      s : fTerm
      l : List Ex
      if a.whole = 0 then l := nil else l := [a.whole :: Ex]
      for s in a.fract repeat
        s.den = 1$FRR => l := cons(s.num :: Ex, l)
        l := cons(s.num :: Ex /  s.den :: Ex, l)
      # l = 1 => first l
      reduce("+", reverse l)