path: root/src/algebra/pfr.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra pfr.spad}
\author{Robert Sutor, Barry Trager}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain PFR PartialFraction}
<<domain PFR PartialFraction>>=
)abbrev domain PFR PartialFraction
++ Author: Robert S. Sutor
++ Date Created: 1986
++ Change History:
++   05/20/91 BMT Converted to the new library
++ Basic Operations: (Field), (Algebra),
++   coerce, compactFraction, firstDenom, firstNumer,
++   nthFractionalTerm, numberOfFractionalTerms, padicallyExpand,
++   padicFraction, partialFraction, wholePart
++ Related Constructors:
++ Also See: ContinuedFraction
++ AMS Classifications:
++ Keywords: partial fraction, factorization, euclidean domain
++ References:
++ 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 \spadfunFrom{padicFraction}{PartialFraction}.  For a
++   general euclidean domain, it is not known how to factor the
++   denominator.  Thus the function \spadfunFrom{partialFraction}{PartialFraction} takes as its
++   second argument an element of \spadtype{Factored(R)}.

PartialFraction(R: EuclideanDomain): Cat == Capsule where
  FRR  ==> Factored R
  SUPR ==> SparseUnivariatePolynomial R

  Cat == Join(Field, Algebra R,CoercibleTo Fraction R) with
    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.

    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.

    firstDenom: % -> FRR
      ++ firstDenom(p) extracts the denominator of the first fractional
      ++ term. This returns 1 if there is no fractional part (use
      ++ \spadfunFrom{wholePart}{PartialFraction} to get the whole part).

    firstNumer:  % -> R
      ++ firstNumer(p) extracts the numerator of the first fractional
      ++ term. This returns 0 if there is no fractional part (use
      ++ \spadfunFrom{wholePart}{PartialFraction} to get the whole part).

    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.

    numberOfFractionalTerms: % -> Integer
      ++ numberOfFractionalTerms(p) computes the number of fractional
      ++ terms in \spad{p}. This returns 0 if there is no fractional
      ++ part.

    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 \spadfunFrom{compactFraction}{PartialFraction} to return to compact form.

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

    wholePart: % -> R
      ++ wholePart(p) extracts the whole part of the partial fraction
      ++ \spad{p}.

  Capsule == 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)

    import %before?: (FRR,FRR) -> Boolean from Foreign Builtin

    -- private function signatures

    copypf: % -> %
    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): Boolean ==
      -- have to wait until FR has < operation
      %before?(s.den,t.den)

    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
      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)
      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)
      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 ==
      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 : %
      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)
      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
      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)

@
\section{package PFRPAC PartialFractionPackage}
<<package PFRPAC PartialFractionPackage>>=
)abbrev package PFRPAC PartialFractionPackage
++ Author: Barry M. Trager
++ Date Created: 1992
++ BasicOperations:
++ Related Constructors: PartialFraction
++ Also See:
++ AMS Classifications:
++ Keywords: partial fraction, factorization, euclidean domain
++ References:
++ Description:
++   The package \spadtype{PartialFractionPackage} gives an easier
++   to use interfact the domain \spadtype{PartialFraction}.
++   The user gives a fraction of polynomials, and a variable and
++   the package converts it to the proper datatype for the
++   \spadtype{PartialFraction} domain.

PartialFractionPackage(R): Cat == Capsule where
--  R : UniqueFactorizationDomain -- not yet supported
   R : Join(EuclideanDomain, CharacteristicZero)
   FPR ==> Fraction Polynomial R
   INDE ==> IndexedExponents Symbol
   PR ==> Polynomial R
   SUP ==> SparseUnivariatePolynomial
   Cat == with
      partialFraction: (FPR, Symbol) -> Any
         ++ partialFraction(rf, var) returns the partial fraction decomposition
         ++ of the rational function rf with respect to the variable var.
      partialFraction: (PR, Factored PR, Symbol) -> Any
         ++ partialFraction(num, facdenom, var) returns the partial fraction
         ++ decomposition of the rational function whose numerator is num and
         ++ whose factored denominator is facdenom with respect to the variable var.
   Capsule == add
      partialFraction(rf, v) ==
         df := factor(denom rf)$MultivariateFactorize(Symbol, INDE,R,PR)
         partialFraction(numer rf, df, v)

      makeSup(p:Polynomial R, v:Symbol) : SparseUnivariatePolynomial FPR ==
         up := univariate(p,v)
         map(#1::FPR,up)$UnivariatePolynomialCategoryFunctions2(PR, SUP PR, FPR, SUP FPR)

      partialFraction(p, facq, v) ==
         up := UnivariatePolynomial(v, Fraction Polynomial R)
         fup := Factored up
         ffact : List(Record(irr:up,pow:Integer))
         ffact:=[[makeSup(u.factor,v) pretend up,u.exponent]
                      for u in factors facq]
         fcont:=makeSup(unit facq,v) pretend up
         nflist:fup := fcont*(*/[primeFactor(ff.irr,ff.pow) for ff in ffact])
         pfup:=partialFraction(makeSup(p,v) pretend up, nflist)$PartialFraction(up)
         coerce(pfup)$AnyFunctions1(PartialFraction up)

@
\section{License}
<<license>>=
--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.
@
<<*>>=
<<license>>

<<domain PFR PartialFraction>>
<<package PFRPAC PartialFractionPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}