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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra sups.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain ISUPS InnerSparseUnivariatePowerSeries}
<<domain ISUPS InnerSparseUnivariatePowerSeries>>=
)abbrev domain ISUPS InnerSparseUnivariatePowerSeries
++ Author: Clifton J. Williamson
++ Date Created: 28 October 1994
++ Date Last Updated: 9 March 1995
++ Basic Operations:
++ Related Domains: SparseUnivariateTaylorSeries, SparseUnivariateLaurentSeries
++   SparseUnivariatePuiseuxSeries
++ Also See:
++ AMS Classifications:
++ Keywords: sparse, series
++ Examples:
++ References:
++ Description: InnerSparseUnivariatePowerSeries is an internal domain
++   used for creating sparse Taylor and Laurent series.
InnerSparseUnivariatePowerSeries(Coef): Exports == Implementation where
  Coef  : Ring
  B    ==> Boolean
  COM  ==> OrderedCompletion Integer
  I    ==> Integer
  L    ==> List
  NNI  ==> NonNegativeInteger
  OUT  ==> OutputForm
  PI   ==> PositiveInteger
  REF  ==> Reference OrderedCompletion Integer
  RN   ==> Fraction Integer
  Term ==> Record(k:Integer,c:Coef)
  SG   ==> String
  ST   ==> Stream Term

  Exports ==> UnivariatePowerSeriesCategory(Coef,Integer) with
    makeSeries: (REF,ST) -> %
      ++ makeSeries(refer,str) creates a power series from the reference
      ++ \spad{refer} and the stream \spad{str}.
    getRef: % -> REF
      ++ getRef(f) returns a reference containing the order to which the
      ++ terms of f have been computed.
    getStream: % -> ST
      ++ getStream(f) returns the stream of terms representing the series f.
    series: ST -> %
      ++ series(st) creates a series from a stream of non-zero terms,
      ++ where a term is an exponent-coefficient pair.  The terms in the
      ++ stream should be ordered by increasing order of exponents.
    monomial?: % -> B
      ++ monomial?(f) tests if f is a single monomial.
    multiplyCoefficients: (I -> Coef,%) -> %
      ++ multiplyCoefficients(fn,f) returns the series
      ++ \spad{sum(fn(n) * an * x^n,n = n0..)},
      ++ where f is the series \spad{sum(an * x^n,n = n0..)}.
    iExquo: (%,%,B) -> Union(%,"failed")
      ++ iExquo(f,g,taylor?) is the quotient of the power series f and g.
      ++ If \spad{taylor?} is \spad{true}, then we must have
      ++ \spad{order(f) >= order(g)}.
    taylorQuoByVar: % -> %
      ++ taylorQuoByVar(a0 + a1 x + a2 x**2 + ...)
      ++ returns \spad{a1 + a2 x + a3 x**2 + ...}
    iCompose: (%,%) -> %
      ++ iCompose(f,g) returns \spad{f(g(x))}.  This is an internal function
      ++ which should only be called for Taylor series \spad{f(x)} and
      ++ \spad{g(x)} such that the constant coefficient of \spad{g(x)} is zero.
    seriesToOutputForm: (ST,REF,Symbol,Coef,RN) -> OutputForm
      ++ seriesToOutputForm(st,refer,var,cen,r) prints the series
      ++ \spad{f((var - cen)^r)}.
    if Coef has Algebra Fraction Integer then
      integrate: % -> %
        ++ integrate(f(x)) returns an anti-derivative of the power series
        ++ \spad{f(x)} with constant coefficient 0.
        ++ Warning: function does not check for a term of degree -1.
      cPower: (%,Coef) -> %
        ++ cPower(f,r) computes \spad{f^r}, where f has constant coefficient 1.
        ++ For use when the coefficient ring is commutative.
      cRationalPower: (%,RN) -> %
        ++ cRationalPower(f,r) computes \spad{f^r}.
        ++ For use when the coefficient ring is commutative.
      cExp: % -> %
        ++ cExp(f) computes the exponential of the power series f.
        ++ For use when the coefficient ring is commutative.
      cLog: % -> %
        ++ cLog(f) computes the logarithm of the power series f.
        ++ For use when the coefficient ring is commutative.
      cSin: % -> %
        ++ cSin(f) computes the sine of the power series f.
        ++ For use when the coefficient ring is commutative.
      cCos: % -> %
        ++ cCos(f) computes the cosine of the power series f.
        ++ For use when the coefficient ring is commutative.
      cTan: % -> %
        ++ cTan(f) computes the tangent of the power series f.
        ++ For use when the coefficient ring is commutative.
      cCot: % -> %
        ++ cCot(f) computes the cotangent of the power series f.
        ++ For use when the coefficient ring is commutative.
      cSec: % -> %
        ++ cSec(f) computes the secant of the power series f.
        ++ For use when the coefficient ring is commutative.
      cCsc: % -> %
        ++ cCsc(f) computes the cosecant of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAsin: % -> %
        ++ cAsin(f) computes the arcsine of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAcos: % -> %
        ++ cAcos(f) computes the arccosine of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAtan: % -> %
        ++ cAtan(f) computes the arctangent of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAcot: % -> %
        ++ cAcot(f) computes the arccotangent of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAsec: % -> %
        ++ cAsec(f) computes the arcsecant of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAcsc: % -> %
        ++ cAcsc(f) computes the arccosecant of the power series f.
        ++ For use when the coefficient ring is commutative.
      cSinh: % -> %
        ++ cSinh(f) computes the hyperbolic sine of the power series f.
        ++ For use when the coefficient ring is commutative.
      cCosh: % -> %
        ++ cCosh(f) computes the hyperbolic cosine of the power series f.
        ++ For use when the coefficient ring is commutative.
      cTanh: % -> %
        ++ cTanh(f) computes the hyperbolic tangent of the power series f.
        ++ For use when the coefficient ring is commutative.
      cCoth: % -> %
        ++ cCoth(f) computes the hyperbolic cotangent of the power series f.
        ++ For use when the coefficient ring is commutative.
      cSech: % -> %
        ++ cSech(f) computes the hyperbolic secant of the power series f.
        ++ For use when the coefficient ring is commutative.
      cCsch: % -> %
        ++ cCsch(f) computes the hyperbolic cosecant of the power series f.
        ++ For use when the coefficient ring is commutative.
      cAsinh: % -> %
        ++ cAsinh(f) computes the inverse hyperbolic sine of the power
        ++ series f.  For use when the coefficient ring is commutative.
      cAcosh: % -> %
        ++ cAcosh(f) computes the inverse hyperbolic cosine of the power
        ++ series f.  For use when the coefficient ring is commutative.
      cAtanh: % -> %
        ++ cAtanh(f) computes the inverse hyperbolic tangent of the power
        ++ series f.  For use when the coefficient ring is commutative.
      cAcoth: % -> %
        ++ cAcoth(f) computes the inverse hyperbolic cotangent of the power
        ++ series f.  For use when the coefficient ring is commutative.
      cAsech: % -> %
        ++ cAsech(f) computes the inverse hyperbolic secant of the power
        ++ series f.  For use when the coefficient ring is commutative.
      cAcsch: % -> %
        ++ cAcsch(f) computes the inverse hyperbolic cosecant of the power
        ++ series f.  For use when the coefficient ring is commutative.

  Implementation ==> add
    import REF

    Rep := Record(%ord: REF,%str: Stream Term)
    -- when the value of 'ord' is n, this indicates that all non-zero
    -- terms of order up to and including n have been computed;
    -- when 'ord' is plusInfinity, all terms have been computed;
    -- lazy evaluation of 'str' has the side-effect of modifying the value
    -- of 'ord'

--% Local functions

    makeTerm:        (Integer,Coef) -> Term
    getCoef:         Term -> Coef
    getExpon:        Term -> Integer
    iSeries:         (ST,REF) -> ST
    iExtend:         (ST,COM,REF) -> ST
    iTruncate0:      (ST,REF,REF,COM,I,I) -> ST
    iTruncate:       (%,COM,I) -> %
    iCoefficient:    (ST,Integer) -> Coef
    iOrder:          (ST,COM,REF) -> I
    iMap1:           ((Coef,I) -> Coef,I -> I,B,ST,REF,REF,Integer) -> ST
    iMap2:           ((Coef,I) -> Coef,I -> I,B,%) -> %
    iPlus1:          ((Coef,Coef) -> Coef,ST,REF,ST,REF,REF,I) -> ST
    iPlus2:          ((Coef,Coef) -> Coef,%,%) -> %
    productByTerm:   (Coef,I,ST,REF,REF,I) -> ST
    productLazyEval: (ST,REF,ST,REF,COM) -> Void
    iTimes:          (ST,REF,ST,REF,REF,I) -> ST
    iDivide:         (ST,REF,ST,REF,Coef,I,REF,I) -> ST
    divide:          (%,I,%,I,Coef) -> %
    compose0:        (ST,REF,ST,REF,I,%,%,I,REF,I) -> ST
    factorials?:     () -> Boolean
    termOutput:      (RN,Coef,OUT) -> OUT
    showAll?:        () -> Boolean

--% macros

    makeTerm(exp,coef) == [exp,coef]
    getCoef term == term.c
    getExpon term == term.k

    makeSeries(refer,x) == [refer,x]
    getRef ups == ups.%ord
    getStream ups == ups.%str

--% creation and destruction of series

    monomial(coef,expon) ==
      nix : ST := empty()
      st :=
        zero? coef => nix
        concat(makeTerm(expon,coef),nix)
      makeSeries(ref plusInfinity(),st)

    monomial? ups == (not empty? getStream ups) and (empty? rst getStream ups)

    coerce(n:I)    == n :: Coef :: %
    coerce(r:Coef) == monomial(r,0)

    iSeries(x,refer) ==
      empty? x => (setref(refer,plusInfinity()); empty())
      setref(refer,(getExpon frst x) :: COM)
      concat(frst x,iSeries(rst x,refer))

    series(x:ST) ==
      empty? x => 0
      n := getExpon frst x; refer := ref(n :: COM)
      makeSeries(refer,iSeries(x,refer))

--% values

    characteristic == characteristic$Coef

    0 == monomial(0,0)
    1 == monomial(1,0)

    iExtend(st,n,refer) ==
      (deref refer) < n =>
        explicitlyEmpty? st => (setref(refer,plusInfinity()); st)
        explicitEntries? st => iExtend(rst st,n,refer)
        iExtend(lazyEvaluate st,n,refer)
      st

    extend(x,n) == (iExtend(getStream x,n :: COM,getRef x); x)
    complete x  == (iExtend(getStream x,plusInfinity(),getRef x); x)

    iTruncate0(x,xRefer,refer,minExp,maxExp,n) == delay
      explicitlyEmpty? x => (setref(refer,plusInfinity()); empty())
      nn := n :: COM
      while (deref xRefer) < nn repeat lazyEvaluate x
      explicitEntries? x =>
        (nx := getExpon(xTerm := frst x)) > maxExp =>
          (setref(refer,plusInfinity()); empty())
        setref(refer,nx :: COM)
        (nx :: COM) >= minExp =>
          concat(makeTerm(nx,getCoef xTerm),_
                 iTruncate0(rst x,xRefer,refer,minExp,maxExp,nx + 1))
        iTruncate0(rst x,xRefer,refer,minExp,maxExp,nx + 1)
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := retract(deref xRefer)@I
      setref(refer,degr :: COM)
      iTruncate0(x,xRefer,refer,minExp,maxExp,degr + 1)

    iTruncate(ups,minExp,maxExp) ==
      x := getStream ups; xRefer := getRef ups
      explicitlyEmpty? x => 0
      explicitEntries? x =>
        deg := getExpon frst x
        refer := ref((deg - 1) :: COM)
        makeSeries(refer,iTruncate0(x,xRefer,refer,minExp,maxExp,deg))
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := retract(deref xRefer)@I
      refer := ref(degr :: COM)
      makeSeries(refer,iTruncate0(x,xRefer,refer,minExp,maxExp,degr + 1))

    truncate(ups,n) == iTruncate(ups,minusInfinity(),n)
    truncate(ups,n1,n2) ==
      if n1 > n2 then (n1,n2) := (n2,n1)
      iTruncate(ups,n1 :: COM,n2)

    iCoefficient(st,n) ==
      explicitEntries? st =>
        term := frst st
        (expon := getExpon term) > n => 0
        expon = n => getCoef term
        iCoefficient(rst st,n)
      0

    coefficient(x,n)   == (extend(x,n); iCoefficient(getStream x,n))
    elt(x:%,n:Integer) == coefficient(x,n)

    iOrder(st,n,refer) ==
      explicitlyEmpty? st => error "order: series has infinite order"
      explicitEntries? st =>
        ((r := getExpon frst st) :: COM) >= n => retract(n)@Integer
        r
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := retract(deref refer)@I
      (degr :: COM) >= n => retract(n)@Integer
      iOrder(lazyEvaluate st,n,refer)

    order x    == iOrder(getStream x,plusInfinity(),getRef x)
    order(x,n) == iOrder(getStream x,n :: COM,getRef x)

    terms x    == getStream x

--% predicates

    zero? ups ==
      x := getStream ups; ref := getRef ups
      whatInfinity(n := deref ref) = 1 => explicitlyEmpty? x
      count : NNI := _$streamCount$Lisp
      for i in 1..count repeat
        explicitlyEmpty? x => return true
        explicitEntries? x => return false
        lazyEvaluate x
      false

    ups1 = ups2 == zero?(ups1 - ups2)

--% arithmetic

    iMap1(cFcn,eFcn,check?,x,xRefer,refer,n) == delay
      -- when this function is called, all terms in 'x' of order < n have been
      -- computed and we compute the eFcn(n)th order coefficient of the result
      explicitlyEmpty? x => (setref(refer,plusInfinity()); empty())
      -- if terms in 'x' up to order n have not been computed,
      -- apply lazy evaluation
      nn := n :: COM
      while (deref xRefer) < nn repeat lazyEvaluate x
      -- 'x' may now be empty: retest
      explicitlyEmpty? x => (setref(refer,plusInfinity()); empty())
      -- must have nx >= n
      explicitEntries? x =>
        xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm
        newCoef := cFcn(xCoef,nx); m := eFcn nx
        setref(refer,m :: COM)
        not check? =>
          concat(makeTerm(m,newCoef),_
                 iMap1(cFcn,eFcn,check?,rst x,xRefer,refer,nx + 1))
        zero? newCoef => iMap1(cFcn,eFcn,check?,rst x,xRefer,refer,nx + 1)
        concat(makeTerm(m,newCoef),_
               iMap1(cFcn,eFcn,check?,rst x,xRefer,refer,nx + 1))
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := retract(deref xRefer)@I
      setref(refer,eFcn(degr) :: COM)
      iMap1(cFcn,eFcn,check?,x,xRefer,refer,degr + 1)

    iMap2(cFcn,eFcn,check?,ups) ==
      -- 'eFcn' must be a strictly increasing function,
      -- i.e. i < j => eFcn(i) < eFcn(j)
      xRefer := getRef ups; x := getStream ups
      explicitlyEmpty? x => 0
      explicitEntries? x =>
        deg := getExpon frst x
        refer := ref(eFcn(deg - 1) :: COM)
        makeSeries(refer,iMap1(cFcn,eFcn,check?,x,xRefer,refer,deg))
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := retract(deref xRefer)@I
      refer := ref(eFcn(degr) :: COM)
      makeSeries(refer,iMap1(cFcn,eFcn,check?,x,xRefer,refer,degr + 1))

    map(fcn,x) == iMap2(fcn(#1),#1,true,x)
    differentiate x == iMap2(#2 * #1,#1 - 1,true,x)
    multiplyCoefficients(f,x) == iMap2(f(#2) * #1,#1,true,x)
    multiplyExponents(x,n) == iMap2(#1,n * #1,false,x)

    iPlus1(op,x,xRefer,y,yRefer,refer,n) == delay
      -- when this function is called, all terms in 'x' and 'y' of order < n
      -- have been computed and we are computing the nth order coefficient of
      -- the result; note the 'op' is either '+' or '-'
      explicitlyEmpty? x => iMap1(op(0,#1),#1,false,y,yRefer,refer,n)
      explicitlyEmpty? y => iMap1(op(#1,0),#1,false,x,xRefer,refer,n)
      -- if terms up to order n have not been computed,
      -- apply lazy evaluation
      nn := n :: COM
      while (deref xRefer) < nn repeat lazyEvaluate x
      while (deref yRefer) < nn repeat lazyEvaluate y
      -- 'x' or 'y' may now be empty: retest
      explicitlyEmpty? x => iMap1(op(0,#1),#1,false,y,yRefer,refer,n)
      explicitlyEmpty? y => iMap1(op(#1,0),#1,false,x,xRefer,refer,n)
      -- must have nx >= n, ny >= n
      -- both x and y have explicit terms
      explicitEntries?(x) and explicitEntries?(y) =>
        xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm
        yCoef := getCoef(yTerm := frst y); ny := getExpon yTerm
        nx = ny =>
          setref(refer,nx :: COM)
          zero? (coef := op(xCoef,yCoef)) =>
            iPlus1(op,rst x,xRefer,rst y,yRefer,refer,nx + 1)
          concat(makeTerm(nx,coef),_
                 iPlus1(op,rst x,xRefer,rst y,yRefer,refer,nx + 1))
        nx < ny =>
          setref(refer,nx :: COM)
          concat(makeTerm(nx,op(xCoef,0)),_
                 iPlus1(op,rst x,xRefer,y,yRefer,refer,nx + 1))
        setref(refer,ny :: COM)
        concat(makeTerm(ny,op(0,yCoef)),_
               iPlus1(op,x,xRefer,rst y,yRefer,refer,ny + 1))
      -- y has no term of degree n
      explicitEntries? x =>
        xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm
        -- can't have deref(yRefer) = infty unless all terms have been computed
        (degr := retract(deref yRefer)@I) < nx =>
          setref(refer,deref yRefer)
          iPlus1(op,x,xRefer,y,yRefer,refer,degr + 1)
        setref(refer,nx :: COM)
        concat(makeTerm(nx,op(xCoef,0)),_
               iPlus1(op,rst x,xRefer,y,yRefer,refer,nx + 1))
      -- x has no term of degree n
      explicitEntries? y =>
        yCoef := getCoef(yTerm := frst y); ny := getExpon yTerm
        -- can't have deref(xRefer) = infty unless all terms have been computed
        (degr := retract(deref xRefer)@I) < ny =>
          setref(refer,deref xRefer)
          iPlus1(op,x,xRefer,y,yRefer,refer,degr + 1)
        setref(refer,ny :: COM)
        concat(makeTerm(ny,op(0,yCoef)),_
               iPlus1(op,x,xRefer,rst y,yRefer,refer,ny + 1))
      -- neither x nor y has a term of degree n
      setref(refer,xyRef := min(deref xRefer,deref yRefer))
      -- can't have xyRef = infty unless all terms have been computed
      iPlus1(op,x,xRefer,y,yRefer,refer,retract(xyRef)@I + 1)

    iPlus2(op,ups1,ups2) ==
      xRefer := getRef ups1; x := getStream ups1
      xDeg :=
        explicitlyEmpty? x => return map(op(0$Coef,#1),ups2)
        explicitEntries? x => (getExpon frst x) - 1
        -- can't have deref(xRefer) = infty unless all terms have been computed
        retract(deref xRefer)@I
      yRefer := getRef ups2; y := getStream ups2
      yDeg :=
        explicitlyEmpty? y => return map(op(#1,0$Coef),ups1)
        explicitEntries? y => (getExpon frst y) - 1
        -- can't have deref(yRefer) = infty unless all terms have been computed
        retract(deref yRefer)@I
      deg := min(xDeg,yDeg); refer := ref(deg :: COM)
      makeSeries(refer,iPlus1(op,x,xRefer,y,yRefer,refer,deg + 1))

    x + y == iPlus2(#1 + #2,x,y)
    x - y == iPlus2(#1 - #2,x,y)
    - y   == iMap2(_-#1,#1,false,y)

    -- gives correct defaults for I, NNI and PI
    n:I   * x:% == (zero? n => 0; map(n * #1,x))
    n:NNI * x:% == (zero? n => 0; map(n * #1,x))
    n:PI  * x:% == (zero? n => 0; map(n * #1,x))

    productByTerm(coef,expon,x,xRefer,refer,n) ==
      iMap1(coef * #1,#1 + expon,true,x,xRefer,refer,n)

    productLazyEval(x,xRefer,y,yRefer,nn) ==
      explicitlyEmpty?(x) or explicitlyEmpty?(y) => void()
      explicitEntries? x =>
        explicitEntries? y => void()
        xDeg := (getExpon frst x) :: COM
        while (xDeg + deref(yRefer)) < nn repeat lazyEvaluate y

      explicitEntries? y =>
        yDeg := (getExpon frst y) :: COM
        while (yDeg + deref(xRefer)) < nn repeat lazyEvaluate x

      lazyEvaluate x
      -- if x = y, then y may now have explicit entries
      if lazy? y then lazyEvaluate y
      productLazyEval(x,xRefer,y,yRefer,nn)

    iTimes(x,xRefer,y,yRefer,refer,n) == delay
      -- when this function is called, we are computing the nth order
      -- coefficient of the product
      productLazyEval(x,xRefer,y,yRefer,n :: COM)
      explicitlyEmpty?(x) or explicitlyEmpty?(y) =>
        (setref(refer,plusInfinity()); empty())
      -- must have nx + ny >= n
      explicitEntries?(x) and explicitEntries?(y) =>
        xCoef := getCoef(xTerm := frst x); xExpon := getExpon xTerm
        yCoef := getCoef(yTerm := frst y); yExpon := getExpon yTerm
        expon := xExpon + yExpon
        setref(refer,expon :: COM)
        scRefer := ref(expon :: COM)
        scMult := productByTerm(xCoef,xExpon,rst y,yRefer,scRefer,yExpon + 1)
        prRefer := ref(expon :: COM)
        pr := iTimes(rst x,xRefer,y,yRefer,prRefer,expon + 1)
        sm := iPlus1(#1 + #2,scMult,scRefer,pr,prRefer,refer,expon + 1)
        zero?(coef := xCoef * yCoef) => sm
        concat(makeTerm(expon,coef),sm)
      explicitEntries? x =>
        xExpon := getExpon frst x
        -- can't have deref(yRefer) = infty unless all terms have been computed
        degr := retract(deref yRefer)@I
        setref(refer,(xExpon + degr) :: COM)
        iTimes(x,xRefer,y,yRefer,refer,xExpon + degr + 1)
      explicitEntries? y =>
        yExpon := getExpon frst y
        -- can't have deref(xRefer) = infty unless all terms have been computed
        degr := retract(deref xRefer)@I
        setref(refer,(yExpon + degr) :: COM)
        iTimes(x,xRefer,y,yRefer,refer,yExpon + degr + 1)
      -- can't have deref(xRefer) = infty unless all terms have been computed
      xDegr := retract(deref xRefer)@I
      yDegr := retract(deref yRefer)@I
      setref(refer,(xDegr + yDegr) :: COM)
      iTimes(x,xRefer,y,yRefer,refer,xDegr + yDegr + 1)

    ups1:% * ups2:% ==
      xRefer := getRef ups1; x := getStream ups1
      xDeg :=
        explicitlyEmpty? x => return 0
        explicitEntries? x => (getExpon frst x) - 1
        -- can't have deref(xRefer) = infty unless all terms have been computed
        retract(deref xRefer)@I
      yRefer := getRef ups2; y := getStream ups2
      yDeg :=
        explicitlyEmpty? y => return 0
        explicitEntries? y => (getExpon frst y) - 1
        -- can't have deref(yRefer) = infty unless all terms have been computed
        retract(deref yRefer)@I
      deg := xDeg + yDeg + 1; refer := ref(deg :: COM)
      makeSeries(refer,iTimes(x,xRefer,y,yRefer,refer,deg + 1))

    iDivide(x,xRefer,y,yRefer,rym,m,refer,n) == delay
      -- when this function is called, we are computing the nth order
      -- coefficient of the result
      explicitlyEmpty? x => (setref(refer,plusInfinity()); empty())
      -- if terms up to order n - m have not been computed,
      -- apply lazy evaluation
      nm := (n + m) :: COM
      while (deref xRefer) < nm repeat lazyEvaluate x
      -- 'x' may now be empty: retest
      explicitlyEmpty? x => (setref(refer,plusInfinity()); empty())
      -- must have nx >= n + m
      explicitEntries? x =>
        newCoef := getCoef(xTerm := frst x) * rym; nx := getExpon xTerm
        prodRefer := ref(nx :: COM)
        prod := productByTerm(-newCoef,nx - m,rst y,yRefer,prodRefer,1)
        sumRefer := ref(nx :: COM)
        sum := iPlus1(#1 + #2,rst x,xRefer,prod,prodRefer,sumRefer,nx + 1)
        setref(refer,(nx - m) :: COM); term := makeTerm(nx - m,newCoef)
        concat(term,iDivide(sum,sumRefer,y,yRefer,rym,m,refer,nx - m + 1))
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := retract(deref xRefer)@I
      setref(refer,(degr - m) :: COM)
      iDivide(x,xRefer,y,yRefer,rym,m,refer,degr - m + 1)

    divide(ups1,deg1,ups2,deg2,r) ==
      xRefer := getRef ups1; x := getStream ups1
      yRefer := getRef ups2; y := getStream ups2
      refer := ref((deg1 - deg2) :: COM)
      makeSeries(refer,iDivide(x,xRefer,y,yRefer,r,deg2,refer,deg1 - deg2 + 1))

    iExquo(ups1,ups2,taylor?) ==
      xRefer := getRef ups1; x := getStream ups1
      yRefer := getRef ups2; y := getStream ups2
      n : I := 0
      -- try to find first non-zero term in y
      -- give up after 1000 lazy evaluations
      while not explicitEntries? y repeat
        explicitlyEmpty? y => return "failed"
        lazyEvaluate y
        (n := n + 1) > 1000 => return "failed"
      yCoef := getCoef(yTerm := frst y); ny := getExpon yTerm
      (ry := recip yCoef) case "failed" => "failed"
      nn := ny :: COM
      if taylor? then
        while (deref(xRefer) < nn) repeat
          explicitlyEmpty? x => return 0
          explicitEntries? x => return "failed"
          lazyEvaluate x
      -- check if ups2 is a monomial
      empty? rst y => iMap2(#1 * (ry :: Coef),#1 - ny,false,ups1)
      explicitlyEmpty? x => 0
      nx :=
        explicitEntries? x =>
          ((deg := getExpon frst x) < ny) and taylor? => return "failed"
          deg - 1
        -- can't have deref(xRefer) = infty unless all terms have been computed
        retract(deref xRefer)@I
      divide(ups1,nx,ups2,ny,ry :: Coef)

    taylorQuoByVar ups ==
      iMap2(#1,#1 - 1,false,ups - monomial(coefficient(ups,0),0))

    compose0(x,xRefer,y,yRefer,yOrd,y1,yn0,n0,refer,n) == delay
      -- when this function is called, we are computing the nth order
      -- coefficient of the composite
      explicitlyEmpty? x => (setref(refer,plusInfinity()); empty())
      -- if terms in 'x' up to order n have not been computed,
      -- apply lazy evaluation
      nn := n :: COM; yyOrd := yOrd :: COM
      while (yyOrd * deref(xRefer)) < nn repeat lazyEvaluate x
      explicitEntries? x =>
        xCoef := getCoef(xTerm := frst x); n1 := getExpon xTerm
        zero? n1 =>
          setref(refer,n1 :: COM)
          concat(makeTerm(n1,xCoef),_
                 compose0(rst x,xRefer,y,yRefer,yOrd,y1,yn0,n0,refer,n1 + 1))
        yn1 := yn0 * y1 ** ((n1 - n0) :: NNI)
        z := getStream yn1; zRefer := getRef yn1
        degr := yOrd * n1; prodRefer := ref((degr - 1) :: COM)
        prod := iMap1(xCoef * #1,#1,true,z,zRefer,prodRefer,degr)
        coRefer := ref((degr + yOrd - 1) :: COM)
        co := compose0(rst x,xRefer,y,yRefer,yOrd,y1,yn1,n1,coRefer,degr + yOrd)
        setref(refer,(degr - 1) :: COM)
        iPlus1(#1 + #2,prod,prodRefer,co,coRefer,refer,degr)
      -- can't have deref(xRefer) = infty unless all terms have been computed
      degr := yOrd * (retract(deref xRefer)@I + 1)
      setref(refer,(degr - 1) :: COM)
      compose0(x,xRefer,y,yRefer,yOrd,y1,yn0,n0,refer,degr)

    iCompose(ups1,ups2) ==
      x := getStream ups1; xRefer := getRef ups1
      y := getStream ups2; yRefer := getRef ups2
      -- try to compute the order of 'ups2'
      n : I := _$streamCount$Lisp
      for i in 1..n while not explicitEntries? y repeat
        explicitlyEmpty? y => return coefficient(ups1,0) :: %
        lazyEvaluate y
      explicitlyEmpty? y => coefficient(ups1,0) :: %
      yOrd : I :=
        explicitEntries? y => getExpon frst y
        retract(deref yRefer)@I
      compRefer := ref((-1) :: COM)
      makeSeries(compRefer,_
                 compose0(x,xRefer,y,yRefer,yOrd,ups2,1,0,compRefer,0))

    if Coef has Algebra Fraction Integer then

      integrate x == iMap2(1/(#2 + 1) * #1,#1 + 1,true,x)

--% Fixed point computations

      Ys ==> Y$ParadoxicalCombinatorsForStreams(Term)

      integ0: (ST,REF,REF,I) -> ST
      integ0(x,intRef,ansRef,n) == delay
        nLess1 := (n - 1) :: COM
        while (deref intRef) < nLess1 repeat lazyEvaluate x
        explicitlyEmpty? x => (setref(ansRef,plusInfinity()); empty())
        explicitEntries? x =>
          xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm
          setref(ansRef,(n1 := (nx + 1)) :: COM)
          concat(makeTerm(n1,inv(n1 :: RN) * xCoef),_
                 integ0(rst x,intRef,ansRef,n1))
        -- can't have deref(intRef) = infty unless all terms have been computed
        degr := retract(deref intRef)@I; setref(ansRef,(degr + 1) :: COM)
        integ0(x,intRef,ansRef,degr + 2)

      integ1: (ST,REF,REF) -> ST
      integ1(x,intRef,ansRef) == integ0(x,intRef,ansRef,1)

      lazyInteg: (Coef,() -> ST,REF,REF) -> ST
      lazyInteg(a,xf,intRef,ansRef) ==
        ansStr : ST := integ1(delay xf,intRef,ansRef)
        concat(makeTerm(0,a),ansStr)

      cPower(f,r) ==
        -- computes f^r.  f should have constant coefficient 1.
        fp := differentiate f
        fInv := iExquo(1,f,false) :: %; y := r * fp * fInv
        yRef := getRef y; yStr := getStream y
        intRef := ref((-1) :: COM); ansRef := ref(0 :: COM)
        ansStr := Ys lazyInteg(1,iTimes(#1,ansRef,yStr,yRef,intRef,0),_
                                 intRef,ansRef)
        makeSeries(ansRef,ansStr)

      iExp: (%,Coef) -> %
      iExp(f,cc) ==
        -- computes exp(f).  cc = exp coefficient(f,0)
        fp := differentiate f
        fpRef := getRef fp; fpStr := getStream fp
        intRef := ref((-1) :: COM); ansRef := ref(0 :: COM)
        ansStr := Ys lazyInteg(cc,iTimes(#1,ansRef,fpStr,fpRef,intRef,0),_
                                  intRef,ansRef)
        makeSeries(ansRef,ansStr)

      sincos0: (Coef,Coef,L ST,REF,REF,ST,REF,ST,REF) -> L ST
      sincos0(sinc,cosc,list,sinRef,cosRef,fpStr,fpRef,fpStr2,fpRef2) ==
        sinStr := first list; cosStr := second list
        prodRef1 := ref((-1) :: COM); prodRef2 := ref((-1) :: COM)
        prodStr1 := iTimes(cosStr,cosRef,fpStr,fpRef,prodRef1,0)
        prodStr2 := iTimes(sinStr,sinRef,fpStr2,fpRef2,prodRef2,0)
        [lazyInteg(sinc,prodStr1,prodRef1,sinRef),_
         lazyInteg(cosc,prodStr2,prodRef2,cosRef)]

      iSincos: (%,Coef,Coef,I) -> Record(%sin: %, %cos: %)
      iSincos(f,sinc,cosc,sign) ==
        fp := differentiate f
        fpRef := getRef fp; fpStr := getStream fp
        fp2 := (one? sign => fp; -fp)
        fpRef2 := getRef fp2; fpStr2 := getStream fp2
        sinRef := ref(0 :: COM); cosRef := ref(0 :: COM)
        sincos :=
          Ys(sincos0(sinc,cosc,#1,sinRef,cosRef,fpStr,fpRef,fpStr2,fpRef2),2)
        sinStr := (zero? sinc => rst first sincos; first sincos)
        cosStr := (zero? cosc => rst second sincos; second sincos)
        [makeSeries(sinRef,sinStr),makeSeries(cosRef,cosStr)]

      tan0: (Coef,ST,REF,ST,REF,I) -> ST
      tan0(cc,ansStr,ansRef,fpStr,fpRef,sign) ==
        sqRef := ref((-1) :: COM)
        sqStr := iTimes(ansStr,ansRef,ansStr,ansRef,sqRef,0)
        one : % := 1; oneStr := getStream one; oneRef := getRef one
        yRef := ref((-1) :: COM)
        yStr : ST :=
          one? sign => iPlus1(#1 + #2,oneStr,oneRef,sqStr,sqRef,yRef,0)
          iPlus1(#1 - #2,oneStr,oneRef,sqStr,sqRef,yRef,0)
        intRef := ref((-1) :: COM)
        lazyInteg(cc,iTimes(yStr,yRef,fpStr,fpRef,intRef,0),intRef,ansRef)

      iTan: (%,%,Coef,I) -> %
      iTan(f,fp,cc,sign) ==
        -- computes the tangent (and related functions) of f.
        fpRef := getRef fp; fpStr := getStream fp
        ansRef := ref(0 :: COM)
        ansStr := Ys tan0(cc,#1,ansRef,fpStr,fpRef,sign)
        zero? cc => makeSeries(ansRef,rst ansStr)
        makeSeries(ansRef,ansStr)

--% Error Reporting

      TRCONST : SG := "series expansion involves transcendental constants"
      NPOWERS : SG := "series expansion has terms of negative degree"
      FPOWERS : SG := "series expansion has terms of fractional degree"
      MAYFPOW : SG := "series expansion may have terms of fractional degree"
      LOGS : SG := "series expansion has logarithmic term"
      NPOWLOG : SG :=
         "series expansion has terms of negative degree or logarithmic term"
      NOTINV : SG := "leading coefficient not invertible"

--% Rational powers and transcendental functions

      orderOrFailed : % -> Union(I,"failed")
      orderOrFailed uts ==
      -- returns the order of x or "failed"
      -- if -1 is returned, the series is identically zero
        x := getStream uts
        for n in 0..1000 repeat
          explicitlyEmpty? x => return -1
          explicitEntries? x => return getExpon frst x
          lazyEvaluate x
        "failed"

      RATPOWERS : Boolean := Coef has "**": (Coef,RN) -> Coef
      TRANSFCN  : Boolean := Coef has TranscendentalFunctionCategory

      cRationalPower(uts,r) ==
        (ord0 := orderOrFailed uts) case "failed" =>
          error "**: series with many leading zero coefficients"
        order := ord0 :: I
        (n := order exquo denom(r)) case "failed" =>
          error "**: rational power does not exist"
        cc := coefficient(uts,order)
        (ccInv := recip cc) case "failed" => error concat("**: ",NOTINV)
        ccPow :=
          one? cc => cc
          one? denom r =>
            not negative?(num := numer r) => cc ** (num :: NNI)
            (ccInv :: Coef) ** ((-num) :: NNI)
          RATPOWERS => cc ** r
          error "** rational power of coefficient undefined"
        uts1 := (ccInv :: Coef) * uts
        uts2 := uts1 * monomial(1,-order)
        monomial(ccPow,(n :: I) * numer(r)) * cPower(uts2,r :: Coef)

      cExp uts ==
        zero?(cc := coefficient(uts,0)) => iExp(uts,1)
        TRANSFCN => iExp(uts,exp cc)
        error concat("exp: ",TRCONST)

      cLog uts ==
        zero?(cc := coefficient(uts,0)) =>
          error "log: constant coefficient should not be 0"
        one? cc => integrate(differentiate(uts) * (iExquo(1,uts,true) :: %))
        TRANSFCN =>
          y := iExquo(1,uts,true) :: %
          (log(cc) :: %) + integrate(y * differentiate(uts))
        error concat("log: ",TRCONST)

      sincos: % -> Record(%sin: %, %cos: %)
      sincos uts ==
        zero?(cc := coefficient(uts,0)) => iSincos(uts,0,1,-1)
        TRANSFCN => iSincos(uts,sin cc,cos cc,-1)
        error concat("sincos: ",TRCONST)

      cSin uts == sincos(uts).%sin
      cCos uts == sincos(uts).%cos

      cTan uts ==
        zero?(cc := coefficient(uts,0)) => iTan(uts,differentiate uts,0,1)
        TRANSFCN => iTan(uts,differentiate uts,tan cc,1)
        error concat("tan: ",TRCONST)

      cCot uts ==
        zero? uts => error "cot: cot(0) is undefined"
        zero?(cc := coefficient(uts,0)) => error error concat("cot: ",NPOWERS)
        TRANSFCN => iTan(uts,-differentiate uts,cot cc,1)
        error concat("cot: ",TRCONST)

      cSec uts ==
        zero?(cc := coefficient(uts,0)) => iExquo(1,cCos uts,true) :: %
        TRANSFCN =>
          cosUts := cCos uts
          zero? coefficient(cosUts,0) => error concat("sec: ",NPOWERS)
          iExquo(1,cosUts,true) :: %
        error concat("sec: ",TRCONST)

      cCsc uts ==
        zero? uts => error "csc: csc(0) is undefined"
        TRANSFCN =>
          sinUts := cSin uts
          zero? coefficient(sinUts,0) => error concat("csc: ",NPOWERS)
          iExquo(1,sinUts,true) :: %
        error concat("csc: ",TRCONST)

      cAsin uts ==
        zero?(cc := coefficient(uts,0)) =>
          integrate(cRationalPower(1 - uts*uts,-1/2) * differentiate(uts))
        TRANSFCN =>
          x := 1 - uts * uts
          cc = 1 or cc = -1 =>
            -- compute order of 'x'
            (ord := orderOrFailed x) case "failed" =>
              error concat("asin: ",MAYFPOW)
            (order := ord :: I) = -1 => return asin(cc) :: %
            odd? order => error concat("asin: ",FPOWERS)
            c0 := asin(cc) :: %
            c0 + integrate(cRationalPower(x,-1/2) * differentiate(uts))
          c0 := asin(cc) :: %
          c0 + integrate(cRationalPower(x,-1/2) * differentiate(uts))
        error concat("asin: ",TRCONST)

      cAcos uts ==
        zero? uts =>
          TRANSFCN => acos(0)$Coef :: %
          error concat("acos: ",TRCONST)
        TRANSFCN =>
          x := 1 - uts * uts
          cc := coefficient(uts,0)
          cc = 1 or cc = -1 =>
            -- compute order of 'x'
            (ord := orderOrFailed x) case "failed" =>
              error concat("acos: ",MAYFPOW)
            (order := ord :: I) = -1 => return acos(cc) :: %
            odd? order => error concat("acos: ",FPOWERS)
            c0 := acos(cc) :: %
            c0 + integrate(-cRationalPower(x,-1/2) * differentiate(uts))
          c0 := acos(cc) :: %
          c0 + integrate(-cRationalPower(x,-1/2) * differentiate(uts))
        error concat("acos: ",TRCONST)

      cAtan uts ==
        zero?(cc := coefficient(uts,0)) =>
          y := iExquo(1,(1 :: %) + uts*uts,true) :: %
          integrate(y * (differentiate uts))
        TRANSFCN =>
          (y := iExquo(1,(1 :: %) + uts*uts,true)) case "failed" =>
            error concat("atan: ",LOGS)
          (atan(cc) :: %) + integrate((y :: %) * (differentiate uts))
        error concat("atan: ",TRCONST)

      cAcot uts ==
        TRANSFCN =>
          (y := iExquo(1,(1 :: %) + uts*uts,true)) case "failed" =>
            error concat("acot: ",LOGS)
          cc := coefficient(uts,0)
          (acot(cc) :: %) + integrate(-(y :: %) * (differentiate uts))
        error concat("acot: ",TRCONST)

      cAsec uts ==
        zero?(cc := coefficient(uts,0)) =>
          error "asec: constant coefficient should not be 0"
        TRANSFCN =>
          x := uts * uts - 1
          y :=
            cc = 1 or cc = -1 =>
              -- compute order of 'x'
              (ord := orderOrFailed x) case "failed" =>
                error concat("asec: ",MAYFPOW)
              (order := ord :: I) = -1 => return asec(cc) :: %
              odd? order => error concat("asec: ",FPOWERS)
              cRationalPower(x,-1/2) * differentiate(uts)
            cRationalPower(x,-1/2) * differentiate(uts)
          (z := iExquo(y,uts,true)) case "failed" =>
            error concat("asec: ",NOTINV)
          (asec(cc) :: %) + integrate(z :: %)
        error concat("asec: ",TRCONST)

      cAcsc uts ==
        zero?(cc := coefficient(uts,0)) =>
          error "acsc: constant coefficient should not be 0"
        TRANSFCN =>
          x := uts * uts - 1
          y :=
            cc = 1 or cc = -1 =>
              -- compute order of 'x'
              (ord := orderOrFailed x) case "failed" =>
                error concat("acsc: ",MAYFPOW)
              (order := ord :: I) = -1 => return acsc(cc) :: %
              odd? order => error concat("acsc: ",FPOWERS)
              -cRationalPower(x,-1/2) * differentiate(uts)
            -cRationalPower(x,-1/2) * differentiate(uts)
          (z := iExquo(y,uts,true)) case "failed" =>
            error concat("asec: ",NOTINV)
          (acsc(cc) :: %) + integrate(z :: %)
        error concat("acsc: ",TRCONST)

      sinhcosh: % -> Record(%sinh: %, %cosh: %)
      sinhcosh uts ==
        zero?(cc := coefficient(uts,0)) =>
          tmp := iSincos(uts,0,1,1)
          [tmp.%sin,tmp.%cos]
        TRANSFCN =>
          tmp := iSincos(uts,sinh cc,cosh cc,1)
          [tmp.%sin,tmp.%cos]
        error concat("sinhcosh: ",TRCONST)

      cSinh uts == sinhcosh(uts).%sinh
      cCosh uts == sinhcosh(uts).%cosh

      cTanh uts ==
        zero?(cc := coefficient(uts,0)) => iTan(uts,differentiate uts,0,-1)
        TRANSFCN => iTan(uts,differentiate uts,tanh cc,-1)
        error concat("tanh: ",TRCONST)

      cCoth uts ==
        tanhUts := cTanh uts
        zero? tanhUts => error "coth: coth(0) is undefined"
        zero? coefficient(tanhUts,0) => error concat("coth: ",NPOWERS)
        iExquo(1,tanhUts,true) :: %

      cSech uts ==
        coshUts := cCosh uts
        zero? coefficient(coshUts,0) => error concat("sech: ",NPOWERS)
        iExquo(1,coshUts,true) :: %

      cCsch uts ==
        sinhUts := cSinh uts
        zero? coefficient(sinhUts,0) => error concat("csch: ",NPOWERS)
        iExquo(1,sinhUts,true) :: %

      cAsinh uts ==
        x := 1 + uts * uts
        zero?(cc := coefficient(uts,0)) => cLog(uts + cRationalPower(x,1/2))
        TRANSFCN =>
          (ord := orderOrFailed x) case "failed" =>
            error concat("asinh: ",MAYFPOW)
          (order := ord :: I) = -1 => return asinh(cc) :: %
          odd? order => error concat("asinh: ",FPOWERS)
          -- the argument to 'log' must have a non-zero constant term
          cLog(uts + cRationalPower(x,1/2))
        error concat("asinh: ",TRCONST)

      cAcosh uts ==
        zero? uts =>
          TRANSFCN => acosh(0)$Coef :: %
          error concat("acosh: ",TRCONST)
        TRANSFCN =>
          cc := coefficient(uts,0); x := uts*uts - 1
          cc = 1 or cc = -1 =>
            -- compute order of 'x'
            (ord := orderOrFailed x) case "failed" =>
              error concat("acosh: ",MAYFPOW)
            (order := ord :: I) = -1 => return acosh(cc) :: %
            odd? order => error concat("acosh: ",FPOWERS)
            -- the argument to 'log' must have a non-zero constant term
            cLog(uts + cRationalPower(x,1/2))
          cLog(uts + cRationalPower(x,1/2))
        error concat("acosh: ",TRCONST)

      cAtanh uts ==
        half := inv(2 :: RN) :: Coef
        zero?(cc := coefficient(uts,0)) =>
          half * (cLog(1 + uts) - cLog(1 - uts))
        TRANSFCN =>
          cc = 1 or cc = -1 => error concat("atanh: ",LOGS)
          half * (cLog(1 + uts) - cLog(1 - uts))
        error concat("atanh: ",TRCONST)

      cAcoth uts ==
        zero? uts =>
          TRANSFCN => acoth(0)$Coef :: %
          error concat("acoth: ",TRCONST)
        TRANSFCN =>
          cc := coefficient(uts,0); half := inv(2 :: RN) :: Coef
          cc = 1 or cc = -1 => error concat("acoth: ",LOGS)
          half * (cLog(uts + 1) - cLog(uts - 1))
        error concat("acoth: ",TRCONST)

      cAsech uts ==
        zero? uts => error "asech: asech(0) is undefined"
        TRANSFCN =>
          zero?(cc := coefficient(uts,0)) =>
            error concat("asech: ",NPOWLOG)
          x := 1 - uts * uts
          cc = 1 or cc = -1 =>
            -- compute order of 'x'
            (ord := orderOrFailed x) case "failed" =>
              error concat("asech: ",MAYFPOW)
            (order := ord :: I) = -1 => return asech(cc) :: %
            odd? order => error concat("asech: ",FPOWERS)
            (utsInv := iExquo(1,uts,true)) case "failed" =>
              error concat("asech: ",NOTINV)
            cLog((1 + cRationalPower(x,1/2)) * (utsInv :: %))
          (utsInv := iExquo(1,uts,true)) case "failed" =>
            error concat("asech: ",NOTINV)
          cLog((1 + cRationalPower(x,1/2)) * (utsInv :: %))
        error concat("asech: ",TRCONST)

      cAcsch uts ==
        zero? uts => error "acsch: acsch(0) is undefined"
        TRANSFCN =>
          zero?(cc := coefficient(uts,0)) => error concat("acsch: ",NPOWLOG)
          x := uts * uts + 1
          -- compute order of 'x'
          (ord := orderOrFailed x) case "failed" =>
            error concat("acsc: ",MAYFPOW)
          (order := ord :: I) = -1 => return acsch(cc) :: %
          odd? order => error concat("acsch: ",FPOWERS)
          (utsInv := iExquo(1,uts,true)) case "failed" =>
            error concat("acsch: ",NOTINV)
          cLog((1 + cRationalPower(x,1/2)) * (utsInv :: %))
        error concat("acsch: ",TRCONST)

--% Output forms

    -- check a global Lisp variable
    factorials?() == false

    termOutput(k,c,vv) ==
    -- creates a term c * vv ** k
      k = 0 => c :: OUT
      mon := (k = 1 => vv; vv ** (k :: OUT))
--       if factorials?() and k > 1 then
--         c := factorial(k)$IntegerCombinatoricFunctions * c
--         mon := mon / hconcat(k :: OUT,"!" :: OUT)
      c = 1 => mon
      c = -1 => -mon
      (c :: OUT) * mon

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

    seriesToOutputForm(st,refer,var,cen,r) ==
      vv :=
        zero? cen => var :: OUT
        paren(var :: OUT - cen :: OUT)
      l : L OUT := empty()
      while explicitEntries? st repeat
        term := frst st
        l := concat(termOutput(getExpon(term) * r,getCoef term,vv),l)
        st := rst st
      l :=
        explicitlyEmpty? st => l
        (deg := retractIfCan(deref refer)@Union(I,"failed")) case I =>
          concat(prefix("O" :: OUT,[vv ** ((((deg :: I) + 1) * r) :: OUT)]),l)
        l
      empty? l => (0$Coef) :: OUT
      reduce("+",reverse! l)

@
\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 ISUPS InnerSparseUnivariatePowerSeries>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}