Initial population.

1 files changed, 1114 insertions, 0 deletions

diff --git a/src/algebra/sups.spad.pamphlet b/src/algebra/sups.spad.pamphlet
new file mode 100644
index 00000000..289a0f64
--- /dev/null
+++ b/src/algebra/sups.spad.pamphlet
@@ -0,0 +1,1114 @@
+\documentclass{article}
+\usepackage{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 => (setelt(refer,plusInfinity()); empty())
+      setelt(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) ==
+      (elt refer) < n =>
+        explicitlyEmpty? st => (setelt(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 => (setelt(refer,plusInfinity()); empty())
+      nn := n :: COM
+      while (elt xRefer) < nn repeat lazyEvaluate x
+      explicitEntries? x =>
+        (nx := getExpon(xTerm := frst x)) > maxExp =>
+          (setelt(refer,plusInfinity()); empty())
+        setelt(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 elt(xRefer) = infty unless all terms have been computed
+      degr := retract(elt xRefer)@I
+      setelt(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 elt(xRefer) = infty unless all terms have been computed
+      degr := retract(elt 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 elt(xRefer) = infty unless all terms have been computed
+      degr := retract(elt 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 := elt 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 => (setelt(refer,plusInfinity()); empty())
+      -- if terms in 'x' up to order n have not been computed,
+      -- apply lazy evaluation
+      nn := n :: COM
+      while (elt xRefer) < nn repeat lazyEvaluate x
+      -- 'x' may now be empty: retest
+      explicitlyEmpty? x => (setelt(refer,plusInfinity()); empty())
+      -- must have nx >= n
+      explicitEntries? x =>
+        xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm
+        newCoef := cFcn(xCoef,nx); m := eFcn nx
+        setelt(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 elt(xRefer) = infty unless all terms have been computed
+      degr := retract(elt xRefer)@I
+      setelt(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 elt(xRefer) = infty unless all terms have been computed
+      degr := retract(elt 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 (elt xRefer) < nn repeat lazyEvaluate x
+      while (elt 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 =>
+          setelt(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 =>
+          setelt(refer,nx :: COM)
+          concat(makeTerm(nx,op(xCoef,0)),_
+                 iPlus1(op,rst x,xRefer,y,yRefer,refer,nx + 1))
+        setelt(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 elt(yRefer) = infty unless all terms have been computed
+        (degr := retract(elt yRefer)@I) < nx =>
+          setelt(refer,elt yRefer)
+          iPlus1(op,x,xRefer,y,yRefer,refer,degr + 1)
+        setelt(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 elt(xRefer) = infty unless all terms have been computed
+        (degr := retract(elt xRefer)@I) < ny =>
+          setelt(refer,elt xRefer)
+          iPlus1(op,x,xRefer,y,yRefer,refer,degr + 1)
+        setelt(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
+      setelt(refer,xyRef := min(elt xRefer,elt 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 elt(xRefer) = infty unless all terms have been computed
+        retract(elt 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 elt(yRefer) = infty unless all terms have been computed
+        retract(elt 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 + elt(yRefer)) < nn repeat lazyEvaluate y
+        void()
+      explicitEntries? y =>
+        yDeg := (getExpon frst y) :: COM
+        while (yDeg + elt(xRefer)) < nn repeat lazyEvaluate x
+        void()
+      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) =>
+        (setelt(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
+        setelt(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 elt(yRefer) = infty unless all terms have been computed
+        degr := retract(elt yRefer)@I
+        setelt(refer,(xExpon + degr) :: COM)
+        iTimes(x,xRefer,y,yRefer,refer,xExpon + degr + 1)
+      explicitEntries? y =>
+        yExpon := getExpon frst y
+        -- can't have elt(xRefer) = infty unless all terms have been computed
+        degr := retract(elt xRefer)@I
+        setelt(refer,(yExpon + degr) :: COM)
+        iTimes(x,xRefer,y,yRefer,refer,yExpon + degr + 1)
+      -- can't have elt(xRefer) = infty unless all terms have been computed
+      xDegr := retract(elt xRefer)@I
+      yDegr := retract(elt yRefer)@I
+      setelt(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 elt(xRefer) = infty unless all terms have been computed
+        retract(elt xRefer)@I
+      yRefer := getRef ups2; y := getStream ups2
+      yDeg :=
+        explicitlyEmpty? y => return 0
+        explicitEntries? y => (getExpon frst y) - 1
+        -- can't have elt(yRefer) = infty unless all terms have been computed
+        retract(elt 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 => (setelt(refer,plusInfinity()); empty())
+      -- if terms up to order n - m have not been computed,
+      -- apply lazy evaluation
+      nm := (n + m) :: COM
+      while (elt xRefer) < nm repeat lazyEvaluate x
+      -- 'x' may now be empty: retest
+      explicitlyEmpty? x => (setelt(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)
+        setelt(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 elt(xRefer) = infty unless all terms have been computed
+      degr := retract(elt xRefer)@I
+      setelt(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 (elt(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 elt(xRefer) = infty unless all terms have been computed
+        retract(elt 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 => (setelt(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 * elt(xRefer)) < nn repeat lazyEvaluate x
+      explicitEntries? x =>
+        xCoef := getCoef(xTerm := frst x); n1 := getExpon xTerm
+        zero? n1 =>
+          setelt(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)
+        setelt(refer,(degr - 1) :: COM)
+        iPlus1(#1 + #2,prod,prodRefer,co,coRefer,refer,degr)
+      -- can't have elt(xRefer) = infty unless all terms have been computed
+      degr := yOrd * (retract(elt xRefer)@I + 1)
+      setelt(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 => coefficient(ups1,0) :: %
+        lazyEvaluate y
+      explicitlyEmpty? y => coefficient(ups1,0) :: %
+      yOrd : I :=
+        explicitEntries? y => getExpon frst y
+        retract(elt 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 (elt intRef) < nLess1 repeat lazyEvaluate x
+        explicitlyEmpty? x => (setelt(ansRef,plusInfinity()); empty())
+        explicitEntries? x =>
+          xCoef := getCoef(xTerm := frst x); nx := getExpon xTerm
+          setelt(ansRef,(n1 := (nx + 1)) :: COM)
+          concat(makeTerm(n1,inv(n1 :: RN) * xCoef),_
+                 integ0(rst x,intRef,ansRef,n1))
+        -- can't have elt(intRef) = infty unless all terms have been computed
+        degr := retract(elt intRef)@I; setelt(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)
+        fp2 := ((sign = 1) => 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)
+          (sign = 1) => 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
+          (cc = 1) => cc
+--          one? denom r =>
+          (denom r) = 1 =>
+            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) :: %))
+        (cc = 1) => 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(elt 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}