Initial population.

1 files changed, 812 insertions, 0 deletions

diff --git a/src/algebra/fs2ups.spad.pamphlet b/src/algebra/fs2ups.spad.pamphlet
new file mode 100644
index 00000000..9751dd40
--- /dev/null
+++ b/src/algebra/fs2ups.spad.pamphlet
@@ -0,0 +1,812 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra fs2ups.spad}
+\author{Clifton J. Williamson}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package FS2UPS FunctionSpaceToUnivariatePowerSeries}
+<<package FS2UPS FunctionSpaceToUnivariatePowerSeries>>=
+)abbrev package FS2UPS FunctionSpaceToUnivariatePowerSeries
+++ Author: Clifton J. Williamson
+++ Date Created: 21 March 1989
+++ Date Last Updated: 2 December 1994
+++ Basic Operations:
+++ Related Domains:
+++ Also See:
+++ AMS Classifications:
+++ Keywords: elementary function, power series
+++ Examples:
+++ References:
+++ Description:
+++   This package converts expressions in some function space to power
+++   series in a variable x with coefficients in that function space.
+++   The function \spadfun{exprToUPS} converts expressions to power series
+++   whose coefficients do not contain the variable x. The function
+++   \spadfun{exprToGenUPS} converts functional expressions to power series
+++   whose coefficients may involve functions of \spad{log(x)}.
+FunctionSpaceToUnivariatePowerSeries(R,FE,Expon,UPS,TRAN,x):_
+ Exports == Implementation where
+  R     : Join(GcdDomain,OrderedSet,RetractableTo Integer,_
+               LinearlyExplicitRingOver Integer)
+  FE    : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
+               FunctionSpace R)
+            with
+              coerce: Expon -> %
+                ++ coerce(e) converts an 'exponent' e to an 'expression'
+  Expon : OrderedRing
+  UPS   : Join(UnivariatePowerSeriesCategory(FE,Expon),Field,_
+               TranscendentalFunctionCategory)
+            with
+              differentiate: % -> %
+                ++ differentiate(x) returns the derivative of x since we 
+                ++ need to be able to differentiate a power series
+              integrate: % -> %
+                ++ integrate(x) returns the integral of x since
+                ++ we need to be able to integrate a power series
+  TRAN  : PartialTranscendentalFunctions UPS
+  x     : Symbol
+  B       ==> Boolean
+  BOP     ==> BasicOperator
+  I       ==> Integer
+  NNI     ==> NonNegativeInteger
+  K       ==> Kernel FE
+  L       ==> List
+  RN      ==> Fraction Integer
+  S       ==> String
+  SY      ==> Symbol
+  PCL     ==> PolynomialCategoryLifting(IndexedExponents K,K,R,SMP,FE)
+  POL     ==> Polynomial R
+  SMP     ==> SparseMultivariatePolynomial(R,K)
+  SUP     ==> SparseUnivariatePolynomial Polynomial R
+  Problem ==> Record(func:String,prob:String)
+  Result  ==> Union(%series:UPS,%problem:Problem)
+  SIGNEF  ==> ElementaryFunctionSign(R,FE)
+
+  Exports ==> with
+    exprToUPS : (FE,B,S) -> Result
+      ++ exprToUPS(fcn,posCheck?,atanFlag) converts the expression
+      ++ \spad{fcn} to a power series.  If \spad{posCheck?} is true,
+      ++ log's of negative numbers are not allowed nor are nth roots of
+      ++ negative numbers with n even.  If \spad{posCheck?} is false,
+      ++ these are allowed.  \spad{atanFlag} determines how the case
+      ++ \spad{atan(f(x))}, where \spad{f(x)} has a pole, will be treated.
+      ++ The possible values of \spad{atanFlag} are \spad{"complex"},
+      ++ \spad{"real: two sides"}, \spad{"real: left side"},
+      ++ \spad{"real: right side"}, and \spad{"just do it"}.
+      ++ If \spad{atanFlag} is \spad{"complex"}, then no series expansion
+      ++ will be computed because, viewed as a function of a complex
+      ++ variable, \spad{atan(f(x))} has an essential singularity.
+      ++ Otherwise, the sign of the leading coefficient of the series
+      ++ expansion of \spad{f(x)} determines the constant coefficient
+      ++ in the series expansion of \spad{atan(f(x))}.  If this sign cannot
+      ++ be determined, a series expansion is computed only when
+      ++ \spad{atanFlag} is \spad{"just do it"}.  When the leading term
+      ++ in the series expansion of \spad{f(x)} is of odd degree (or is a
+      ++ rational degree with odd numerator), then the constant coefficient
+      ++ in the series expansion of \spad{atan(f(x))} for values to the
+      ++ left differs from that for values to the right.  If \spad{atanFlag}
+      ++ is \spad{"real: two sides"}, no series expansion will be computed.
+      ++ If \spad{atanFlag} is \spad{"real: left side"} the constant
+      ++ coefficient for values to the left will be used and if \spad{atanFlag}
+      ++ \spad{"real: right side"} the constant coefficient for values to the
+      ++ right will be used.
+      ++ If there is a problem in converting the function to a power series,
+      ++ a record containing the name of the function that caused the problem
+      ++ and a brief description of the problem is returned.
+      ++ When expanding the expression into a series it is assumed that
+      ++ the series is centered at 0.  For a series centered at a, the
+      ++ user should perform the substitution \spad{x -> x + a} before calling
+      ++ this function.
+
+    exprToGenUPS : (FE,B,S) -> Result
+      ++ exprToGenUPS(fcn,posCheck?,atanFlag) converts the expression
+      ++ \spad{fcn} to a generalized power series.  If \spad{posCheck?}
+      ++ is true, log's of negative numbers are not allowed nor are nth roots
+      ++ of negative numbers with n even. If \spad{posCheck?} is false,
+      ++ these are allowed.  \spad{atanFlag} determines how the case
+      ++ \spad{atan(f(x))}, where \spad{f(x)} has a pole, will be treated.
+      ++ The possible values of \spad{atanFlag} are \spad{"complex"},
+      ++ \spad{"real: two sides"}, \spad{"real: left side"},
+      ++ \spad{"real: right side"}, and \spad{"just do it"}.
+      ++ If \spad{atanFlag} is \spad{"complex"}, then no series expansion
+      ++ will be computed because, viewed as a function of a complex
+      ++ variable, \spad{atan(f(x))} has an essential singularity.
+      ++ Otherwise, the sign of the leading coefficient of the series
+      ++ expansion of \spad{f(x)} determines the constant coefficient
+      ++ in the series expansion of \spad{atan(f(x))}.  If this sign cannot
+      ++ be determined, a series expansion is computed only when
+      ++ \spad{atanFlag} is \spad{"just do it"}.  When the leading term
+      ++ in the series expansion of \spad{f(x)} is of odd degree (or is a
+      ++ rational degree with odd numerator), then the constant coefficient
+      ++ in the series expansion of \spad{atan(f(x))} for values to the
+      ++ left differs from that for values to the right.  If \spad{atanFlag}
+      ++ is \spad{"real: two sides"}, no series expansion will be computed.
+      ++ If \spad{atanFlag} is \spad{"real: left side"} the constant
+      ++ coefficient for values to the left will be used and if \spad{atanFlag}
+      ++ \spad{"real: right side"} the constant coefficient for values to the
+      ++ right will be used.
+      ++ If there is a problem in converting the function to a power
+      ++ series, we return a record containing the name of the function
+      ++ that caused the problem and a brief description of the problem.
+      ++ When expanding the expression into a series it is assumed that
+      ++ the series is centered at 0.  For a series centered at a, the
+      ++ user should perform the substitution \spad{x -> x + a} before calling
+      ++ this function.
+    localAbs: FE -> FE
+      ++ localAbs(fcn) = \spad{abs(fcn)} or \spad{sqrt(fcn**2)} depending
+      ++ on whether or not FE has a function \spad{abs}.  This should be
+      ++ a local function, but the compiler won't allow it.
+
+  Implementation ==> add
+
+    ratIfCan            : FE -> Union(RN,"failed")
+    carefulNthRootIfCan : (UPS,NNI,B,B) -> Result
+    stateProblem        : (S,S) -> Result
+    polyToUPS           : SUP -> UPS
+    listToUPS           : (L FE,(FE,B,S) -> Result,B,S,UPS,(UPS,UPS) -> UPS)_
+                                            -> Result
+    isNonTrivPower      : FE -> Union(Record(val:FE,exponent:I),"failed")
+    powerToUPS          : (FE,I,B,S) -> Result
+    kernelToUPS         : (K,B,S) -> Result
+    nthRootToUPS        : (FE,NNI,B,S) -> Result
+    logToUPS            : (FE,B,S) -> Result
+    atancotToUPS        : (FE,B,S,I) -> Result
+    applyIfCan          : (UPS -> Union(UPS,"failed"),FE,S,B,S) -> Result
+    tranToUPS           : (K,FE,B,S) -> Result
+    powToUPS            : (L FE,B,S) -> Result
+    newElem             : FE -> FE
+    smpElem             : SMP -> FE
+    k2Elem              : K -> FE
+    contOnReals?        : S -> B
+    bddOnReals?         : S -> B
+    iExprToGenUPS       : (FE,B,S) -> Result
+    opsInvolvingX       : FE -> L BOP
+    opInOpList?         : (SY,L BOP) -> B
+    exponential?        : FE -> B
+    productOfNonZeroes? : FE -> B
+    powerToGenUPS       : (FE,I,B,S) -> Result
+    kernelToGenUPS      : (K,B,S) -> Result
+    nthRootToGenUPS     : (FE,NNI,B,S) -> Result
+    logToGenUPS         : (FE,B,S) -> Result
+    expToGenUPS         : (FE,B,S) -> Result
+    expGenUPS           : (UPS,B,S) -> Result
+    atancotToGenUPS     : (FE,FE,B,S,I) -> Result
+    genUPSApplyIfCan    : (UPS -> Union(UPS,"failed"),FE,S,B,S) -> Result
+    applyBddIfCan       : (FE,UPS -> Union(UPS,"failed"),FE,S,B,S) -> Result
+    tranToGenUPS        : (K,FE,B,S) -> Result
+    powToGenUPS         : (L FE,B,S) -> Result
+
+    ZEROCOUNT : I := 1000
+    -- number of zeroes to be removed when taking logs or nth roots
+
+    ratIfCan fcn == retractIfCan(fcn)@Union(RN,"failed")
+
+    carefulNthRootIfCan(ups,n,posCheck?,rightOnly?) ==
+      -- similar to 'nthRootIfCan', but it is fussy about the series
+      -- it takes as an argument.  If 'n' is EVEN and 'posCheck?'
+      -- is truem then the leading coefficient of the series must
+      -- be POSITIVE.  In this case, if 'rightOnly?' is false, the
+      -- order of the series must be zero.  The idea is that the
+      -- series represents a real function of a real variable, and
+      -- we want a unique real nth root defined on a neighborhood
+      -- of zero.
+      n < 1 => error "nthRoot: n must be positive"
+      deg := degree ups
+      if (coef := coefficient(ups,deg)) = 0 then
+        deg := order(ups,deg + ZEROCOUNT :: Expon)
+        (coef := coefficient(ups,deg)) = 0 =>
+          error "log of series with many leading zero coefficients"
+      -- if 'posCheck?' is true, we do not allow nth roots of negative
+      -- numbers when n in even
+      if even?(n :: I) then
+        if posCheck? and ((signum := sign(coef)$SIGNEF) case I) then
+          (signum :: I) = -1 =>
+            return stateProblem("nth root","negative leading coefficient")
+          not rightOnly? and not zero? deg => -- nth root not unique
+            return stateProblem("nth root","series of non-zero order")
+      (ans := nthRootIfCan(ups,n)) case "failed" =>
+        stateProblem("nth root","no nth root")
+      [ans :: UPS]
+
+    stateProblem(function,problem) ==
+      -- records the problem which occured in converting an expression
+      -- to a power series
+      [[function,problem]]
+
+    exprToUPS(fcn,posCheck?,atanFlag) ==
+      -- converts a functional expression to a power series
+      --!! The following line is commented out so that expressions of
+      --!! the form a**b will be normalized to exp(b * log(a)) even if
+      --!! 'a' and 'b' do not involve the limiting variable 'x'.
+      --!!                         - cjw 1 Dec 94
+      --not member?(x,variables fcn) => [monomial(fcn,0)]
+      (poly := retractIfCan(fcn)@Union(POL,"failed")) case POL =>
+        [polyToUPS univariate(poly :: POL,x)]
+      (sum := isPlus fcn) case L(FE) =>
+        listToUPS(sum :: L(FE),exprToUPS,posCheck?,atanFlag,0,#1 + #2)
+      (prod := isTimes fcn) case L(FE) =>
+        listToUPS(prod :: L(FE),exprToUPS,posCheck?,atanFlag,1,#1 * #2)
+      (expt := isNonTrivPower fcn) case Record(val:FE,exponent:I) =>
+        power := expt :: Record(val:FE,exponent:I)
+        powerToUPS(power.val,power.exponent,posCheck?,atanFlag)
+      (ker := retractIfCan(fcn)@Union(K,"failed")) case K =>
+        kernelToUPS(ker :: K,posCheck?,atanFlag)
+      error "exprToUPS: neither a sum, product, power, nor kernel"
+
+    polyToUPS poly ==
+      -- converts a polynomial to a power series
+      zero? poly => 0
+      -- we don't start with 'ans := 0' as this may lead to an
+      -- enormous number of leading zeroes in the power series
+      deg  := degree poly
+      coef := leadingCoefficient(poly) :: FE
+      ans  := monomial(coef,deg :: Expon)$UPS
+      poly := reductum poly
+      while not zero? poly repeat
+        deg  := degree poly
+        coef := leadingCoefficient(poly) :: FE
+        ans  := ans + monomial(coef,deg :: Expon)$UPS
+        poly := reductum poly
+      ans
+
+    listToUPS(list,feToUPS,posCheck?,atanFlag,ans,op) ==
+      -- converts each element of a list of expressions to a power
+      -- series and returns the sum of these series, when 'op' is +
+      -- and 'ans' is 0, or the product of these series, when 'op' is *
+      -- and 'ans' is 1
+      while not null list repeat
+        (term := feToUPS(first list,posCheck?,atanFlag)) case %problem =>
+          return term
+        ans := op(ans,term.%series)
+        list := rest list
+      [ans]
+
+    isNonTrivPower fcn ==
+      -- is the function a power with exponent other than 0 or 1?
+      (expt := isPower fcn) case "failed" => "failed"
+      power := expt :: Record(val:FE,exponent:I)
+--      one? power.exponent => "failed"
+      (power.exponent = 1) => "failed"
+      power
+
+    powerToUPS(fcn,n,posCheck?,atanFlag) ==
+      -- converts an integral power to a power series
+      (b := exprToUPS(fcn,posCheck?,atanFlag)) case %problem => b
+      n > 0 => [(b.%series) ** n]
+      -- check lowest order coefficient when n < 0
+      ups := b.%series; deg := degree ups
+      if (coef := coefficient(ups,deg)) = 0 then
+        deg := order(ups,deg + ZEROCOUNT :: Expon)
+        (coef := coefficient(ups,deg)) = 0 =>
+          error "inverse of series with many leading zero coefficients"
+      [ups ** n]
+
+    kernelToUPS(ker,posCheck?,atanFlag) ==
+      -- converts a kernel to a power series
+      (sym := symbolIfCan(ker)) case Symbol =>
+        (sym :: Symbol) = x => [monomial(1,1)]
+        [monomial(ker :: FE,0)]
+      empty?(args := argument ker) => [monomial(ker :: FE,0)]
+      not member?(x, variables(ker :: FE)) => [monomial(ker :: FE,0)]
+      empty? rest args =>
+        arg := first args
+        is?(ker,"abs" :: Symbol) =>
+          nthRootToUPS(arg*arg,2,posCheck?,atanFlag)
+        is?(ker,"%paren" :: Symbol) => exprToUPS(arg,posCheck?,atanFlag)
+        is?(ker,"log" :: Symbol) => logToUPS(arg,posCheck?,atanFlag)
+        is?(ker,"exp" :: Symbol) =>
+          applyIfCan(expIfCan,arg,"exp",posCheck?,atanFlag)
+        tranToUPS(ker,arg,posCheck?,atanFlag)
+      is?(ker,"%power" :: Symbol) => powToUPS(args,posCheck?,atanFlag)
+      is?(ker,"nthRoot" :: Symbol) =>
+        n := retract(second args)@I
+        nthRootToUPS(first args,n :: NNI,posCheck?,atanFlag)
+      stateProblem(string name ker,"unknown kernel")
+
+    nthRootToUPS(arg,n,posCheck?,atanFlag) ==
+      -- converts an nth root to a power series
+      -- this is not used in the limit package, so the series may
+      -- have non-zero order, in which case nth roots may not be unique
+      (result := exprToUPS(arg,posCheck?,atanFlag)) case %problem => result
+      ans := carefulNthRootIfCan(result.%series,n,posCheck?,false)
+      ans case %problem => ans
+      [ans.%series]
+
+    logToUPS(arg,posCheck?,atanFlag) ==
+      -- converts a logarithm log(f(x)) to a power series
+      -- f(x) must have order 0 and if 'posCheck?' is true,
+      -- then f(x) must have a non-negative leading coefficient
+      (result := exprToUPS(arg,posCheck?,atanFlag)) case %problem => result
+      ups := result.%series
+      not zero? order(ups,1) =>
+        stateProblem("log","series of non-zero order")
+      coef := coefficient(ups,0)
+      -- if 'posCheck?' is true, we do not allow logs of negative numbers
+      if posCheck? then
+        if ((signum := sign(coef)$SIGNEF) case I) then
+          (signum :: I) = -1 =>
+            return stateProblem("log","negative leading coefficient")
+      [logIfCan(ups) :: UPS]
+
+    if FE has abs: FE -> FE then
+      localAbs fcn == abs fcn
+    else
+      localAbs fcn == sqrt(fcn * fcn)
+
+    signOfExpression: FE -> FE
+    signOfExpression arg == localAbs(arg)/arg
+
+    atancotToUPS(arg,posCheck?,atanFlag,plusMinus) ==
+      -- converts atan(f(x)) to a power series
+      (result := exprToUPS(arg,posCheck?,atanFlag)) case %problem => result
+      ups := result.%series; coef := coefficient(ups,0)
+      (ord := order(ups,0)) = 0 and coef * coef = -1 =>
+        -- series involves complex numbers
+        return stateProblem("atan","logarithmic singularity")
+      cc : FE :=
+        ord < 0 =>
+          atanFlag = "complex" =>
+            return stateProblem("atan","essential singularity")
+          (rn := ratIfCan(ord :: FE)) case "failed" =>
+            -- this condition usually won't occur because exponents will
+            -- be integers or rational numbers
+            return stateProblem("atan","branch problem")
+          if (atanFlag = "real: two sides") and (odd? numer(rn :: RN)) then
+            -- expansions to the left and right of zero have different
+            -- constant coefficients
+            return stateProblem("atan","branch problem")
+          lc := coefficient(ups,ord)
+          (signum := sign(lc)$SIGNEF) case "failed" =>
+            -- can't determine sign
+            atanFlag = "just do it" =>
+              plusMinus = 1 => pi()/(2 :: FE)
+              0
+            posNegPi2 := signOfExpression(lc) * pi()/(2 :: FE)
+            plusMinus = 1 => posNegPi2
+            pi()/(2 :: FE) - posNegPi2
+            --return stateProblem("atan","branch problem")
+          left? : B := atanFlag = "real: left side"; n := signum :: Integer
+          (left? and n = 1) or (not left? and n = -1) =>
+            plusMinus = 1 => -pi()/(2 :: FE)
+            pi()
+          plusMinus = 1 => pi()/(2 :: FE)
+          0
+        atan coef
+      [(cc :: UPS) + integrate(plusMinus * differentiate(ups)/(1 + ups*ups))]
+
+    applyIfCan(fcn,arg,fcnName,posCheck?,atanFlag) ==
+      -- converts fcn(arg) to a power series
+      (ups := exprToUPS(arg,posCheck?,atanFlag)) case %problem => ups
+      ans := fcn(ups.%series)
+      ans case "failed" => stateProblem(fcnName,"essential singularity")
+      [ans :: UPS]
+
+    tranToUPS(ker,arg,posCheck?,atanFlag) ==
+      -- converts ker to a power series for certain functions
+      -- in trig or hyperbolic trig categories
+      is?(ker,"sin" :: SY) =>
+        applyIfCan(sinIfCan,arg,"sin",posCheck?,atanFlag)
+      is?(ker,"cos" :: SY) =>
+        applyIfCan(cosIfCan,arg,"cos",posCheck?,atanFlag)
+      is?(ker,"tan" :: SY) =>
+        applyIfCan(tanIfCan,arg,"tan",posCheck?,atanFlag)
+      is?(ker,"cot" :: SY) =>
+        applyIfCan(cotIfCan,arg,"cot",posCheck?,atanFlag)
+      is?(ker,"sec" :: SY) =>
+        applyIfCan(secIfCan,arg,"sec",posCheck?,atanFlag)
+      is?(ker,"csc" :: SY) =>
+        applyIfCan(cscIfCan,arg,"csc",posCheck?,atanFlag)
+      is?(ker,"asin" :: SY) =>
+        applyIfCan(asinIfCan,arg,"asin",posCheck?,atanFlag)
+      is?(ker,"acos" :: SY) =>
+        applyIfCan(acosIfCan,arg,"acos",posCheck?,atanFlag)
+      is?(ker,"atan" :: SY) => atancotToUPS(arg,posCheck?,atanFlag,1)
+      is?(ker,"acot" :: SY) => atancotToUPS(arg,posCheck?,atanFlag,-1)
+      is?(ker,"asec" :: SY) =>
+        applyIfCan(asecIfCan,arg,"asec",posCheck?,atanFlag)
+      is?(ker,"acsc" :: SY) =>
+        applyIfCan(acscIfCan,arg,"acsc",posCheck?,atanFlag)
+      is?(ker,"sinh" :: SY) =>
+        applyIfCan(sinhIfCan,arg,"sinh",posCheck?,atanFlag)
+      is?(ker,"cosh" :: SY) =>
+        applyIfCan(coshIfCan,arg,"cosh",posCheck?,atanFlag)
+      is?(ker,"tanh" :: SY) =>
+        applyIfCan(tanhIfCan,arg,"tanh",posCheck?,atanFlag)
+      is?(ker,"coth" :: SY) =>
+        applyIfCan(cothIfCan,arg,"coth",posCheck?,atanFlag)
+      is?(ker,"sech" :: SY) =>
+        applyIfCan(sechIfCan,arg,"sech",posCheck?,atanFlag)
+      is?(ker,"csch" :: SY) =>
+        applyIfCan(cschIfCan,arg,"csch",posCheck?,atanFlag)
+      is?(ker,"asinh" :: SY) =>
+        applyIfCan(asinhIfCan,arg,"asinh",posCheck?,atanFlag)
+      is?(ker,"acosh" :: SY) =>
+        applyIfCan(acoshIfCan,arg,"acosh",posCheck?,atanFlag)
+      is?(ker,"atanh" :: SY) =>
+        applyIfCan(atanhIfCan,arg,"atanh",posCheck?,atanFlag)
+      is?(ker,"acoth" :: SY) =>
+        applyIfCan(acothIfCan,arg,"acoth",posCheck?,atanFlag)
+      is?(ker,"asech" :: SY) =>
+        applyIfCan(asechIfCan,arg,"asech",posCheck?,atanFlag)
+      is?(ker,"acsch" :: SY) =>
+        applyIfCan(acschIfCan,arg,"acsch",posCheck?,atanFlag)
+      stateProblem(string name ker,"unknown kernel")
+
+    powToUPS(args,posCheck?,atanFlag) ==
+      -- converts a power f(x) ** g(x) to a power series
+      (logBase := logToUPS(first args,posCheck?,atanFlag)) case %problem =>
+        logBase
+      (expon := exprToUPS(second args,posCheck?,atanFlag)) case %problem =>
+        expon
+      ans := expIfCan((expon.%series) * (logBase.%series))
+      ans case "failed" => stateProblem("exp","essential singularity")
+      [ans :: UPS]
+
+-- Generalized power series: power series in x, where log(x) and
+-- bounded functions of x are allowed to appear in the coefficients
+-- of the series.  Used for evaluating REAL limits at x = 0.
+
+    newElem f ==
+    -- rewrites a functional expression; all trig functions are
+    -- expressed in terms of sin and cos; all hyperbolic trig
+    -- functions are expressed in terms of exp
+      smpElem(numer f) / smpElem(denom f)
+
+    smpElem p == map(k2Elem,#1::FE,p)$PCL
+
+    k2Elem k ==
+    -- rewrites a kernel; all trig functions are
+    -- expressed in terms of sin and cos; all hyperbolic trig
+    -- functions are expressed in terms of exp
+      null(args := [newElem a for a in argument k]) => k::FE
+      iez  := inv(ez  := exp(z := first args))
+      sinz := sin z; cosz := cos z
+      is?(k,"tan" :: Symbol)  => sinz / cosz
+      is?(k,"cot" :: Symbol)  => cosz / sinz
+      is?(k,"sec" :: Symbol)  => inv cosz
+      is?(k,"csc" :: Symbol)  => inv sinz
+      is?(k,"sinh" :: Symbol) => (ez - iez) / (2 :: FE)
+      is?(k,"cosh" :: Symbol) => (ez + iez) / (2 :: FE)
+      is?(k,"tanh" :: Symbol) => (ez - iez) / (ez + iez)
+      is?(k,"coth" :: Symbol) => (ez + iez) / (ez - iez)
+      is?(k,"sech" :: Symbol) => 2 * inv(ez + iez)
+      is?(k,"csch" :: Symbol) => 2 * inv(ez - iez)
+      (operator k) args
+
+    CONTFCNS : L S := ["sin","cos","atan","acot","exp","asinh"]
+    -- functions which are defined and continuous at all real numbers
+
+    BDDFCNS : L S := ["sin","cos","atan","acot"]
+    -- functions which are bounded on the reals
+
+    contOnReals? fcn == member?(fcn,CONTFCNS)
+    bddOnReals? fcn  == member?(fcn,BDDFCNS)
+
+    exprToGenUPS(fcn,posCheck?,atanFlag) ==
+      -- converts a functional expression to a generalized power
+      -- series; "generalized" means that log(x) and bounded functions
+      -- of x are allowed to appear in the coefficients of the series
+      iExprToGenUPS(newElem fcn,posCheck?,atanFlag)
+
+    iExprToGenUPS(fcn,posCheck?,atanFlag) ==
+      -- converts a functional expression to a generalized power
+      -- series without first normalizing the expression
+      --!! The following line is commented out so that expressions of
+      --!! the form a**b will be normalized to exp(b * log(a)) even if
+      --!! 'a' and 'b' do not involve the limiting variable 'x'.
+      --!!                         - cjw 1 Dec 94
+      --not member?(x,variables fcn) => [monomial(fcn,0)]
+      (poly := retractIfCan(fcn)@Union(POL,"failed")) case POL =>
+        [polyToUPS univariate(poly :: POL,x)]
+      (sum := isPlus fcn) case L(FE) =>
+        listToUPS(sum :: L(FE),iExprToGenUPS,posCheck?,atanFlag,0,#1 + #2)
+      (prod := isTimes fcn) case L(FE) =>
+        listToUPS(prod :: L(FE),iExprToGenUPS,posCheck?,atanFlag,1,#1 * #2)
+      (expt := isNonTrivPower fcn) case Record(val:FE,exponent:I) =>
+        power := expt :: Record(val:FE,exponent:I)
+        powerToGenUPS(power.val,power.exponent,posCheck?,atanFlag)
+      (ker := retractIfCan(fcn)@Union(K,"failed")) case K =>
+        kernelToGenUPS(ker :: K,posCheck?,atanFlag)
+      error "exprToGenUPS: neither a sum, product, power, nor kernel"
+
+    opsInvolvingX fcn ==
+      opList := [op for k in tower fcn | unary?(op := operator k) _
+                 and member?(x,variables first argument k)]
+      removeDuplicates opList
+
+    opInOpList?(name,opList) ==
+      for op in opList repeat
+        is?(op,name) => return true
+      false
+
+    exponential? fcn ==
+      -- is 'fcn' of the form exp(f)?
+      (ker := retractIfCan(fcn)@Union(K,"failed")) case K =>
+        is?(ker :: K,"exp" :: Symbol)
+      false
+
+    productOfNonZeroes? fcn ==
+      -- is 'fcn' a product of non-zero terms, where 'non-zero'
+      -- means an exponential or a function not involving x
+      exponential? fcn => true
+      (prod := isTimes fcn) case "failed" => false
+      for term in (prod :: L(FE)) repeat
+        (not exponential? term) and member?(x,variables term) =>
+          return false
+      true
+
+    powerToGenUPS(fcn,n,posCheck?,atanFlag) ==
+      -- converts an integral power to a generalized power series
+      -- if n < 0 and the lowest order coefficient of the series
+      -- involves x, we are careful about inverting this coefficient
+      -- the coefficient is inverted only if
+      -- (a) the only function involving x is 'log', or
+      -- (b) the lowest order coefficient is a product of exponentials
+      --     and functions not involving x
+      (b := exprToGenUPS(fcn,posCheck?,atanFlag)) case %problem => b
+      n > 0 => [(b.%series) ** n]
+      -- check lowest order coefficient when n < 0
+      ups := b.%series; deg := degree ups
+      if (coef := coefficient(ups,deg)) = 0 then
+        deg := order(ups,deg + ZEROCOUNT :: Expon)
+        (coef := coefficient(ups,deg)) = 0 =>
+          error "inverse of series with many leading zero coefficients"
+      xOpList := opsInvolvingX coef
+      -- only function involving x is 'log'
+      (null xOpList) => [ups ** n]
+      (null rest xOpList and is?(first xOpList,"log" :: SY)) =>
+        [ups ** n]
+      -- lowest order coefficient is a product of exponentials and
+      -- functions not involving x
+      productOfNonZeroes? coef => [ups ** n]
+      stateProblem("inv","lowest order coefficient involves x")
+
+    kernelToGenUPS(ker,posCheck?,atanFlag) ==
+      -- converts a kernel to a generalized power series
+      (sym := symbolIfCan(ker)) case Symbol =>
+        (sym :: Symbol) = x => [monomial(1,1)]
+        [monomial(ker :: FE,0)]
+      empty?(args := argument ker) => [monomial(ker :: FE,0)]
+      empty? rest args =>
+        arg := first args
+        is?(ker,"abs" :: Symbol) =>
+          nthRootToGenUPS(arg*arg,2,posCheck?,atanFlag)
+        is?(ker,"%paren" :: Symbol) => iExprToGenUPS(arg,posCheck?,atanFlag)
+        is?(ker,"log" :: Symbol) => logToGenUPS(arg,posCheck?,atanFlag)
+        is?(ker,"exp" :: Symbol) => expToGenUPS(arg,posCheck?,atanFlag)
+        tranToGenUPS(ker,arg,posCheck?,atanFlag)
+      is?(ker,"%power" :: Symbol) => powToGenUPS(args,posCheck?,atanFlag)
+      is?(ker,"nthRoot" :: Symbol) =>
+        n := retract(second args)@I
+        nthRootToGenUPS(first args,n :: NNI,posCheck?,atanFlag)
+      stateProblem(string name ker,"unknown kernel")
+
+    nthRootToGenUPS(arg,n,posCheck?,atanFlag) ==
+      -- convert an nth root to a power series
+      -- used for computing right hand limits, so the series may have
+      -- non-zero order, but may not have a negative leading coefficient
+      -- when n is even
+      (result := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem =>
+        result
+      ans := carefulNthRootIfCan(result.%series,n,posCheck?,true)
+      ans case %problem => ans
+      [ans.%series]
+
+    logToGenUPS(arg,posCheck?,atanFlag) ==
+      -- converts a logarithm log(f(x)) to a generalized power series
+      (result := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem =>
+        result
+      ups := result.%series; deg := degree ups
+      if (coef := coefficient(ups,deg)) = 0 then
+        deg := order(ups,deg + ZEROCOUNT :: Expon)
+        (coef := coefficient(ups,deg)) = 0 =>
+          error "log of series with many leading zero coefficients"
+      -- if 'posCheck?' is true, we do not allow logs of negative numbers
+      if posCheck? then
+        if ((signum := sign(coef)$SIGNEF) case I) then
+          (signum :: I) = -1 =>
+            return stateProblem("log","negative leading coefficient")
+      -- create logarithmic term, avoiding log's of negative rationals
+      lt := monomial(coef,deg)$UPS; cen := center lt
+      -- check to see if lowest order coefficient is a negative rational
+      negRat? : Boolean :=
+        ((rat := ratIfCan coef) case RN) =>
+          (rat :: RN) < 0 => true
+          false
+        false
+      logTerm : FE :=
+        mon : FE := (x :: FE) - (cen :: FE)
+        pow : FE := mon ** (deg :: FE)
+        negRat? => log(coef * pow)
+        term1 : FE := (deg :: FE) * log(mon)
+        log(coef) + term1
+      [monomial(logTerm,0) + log(ups/lt)]
+
+    expToGenUPS(arg,posCheck?,atanFlag) ==
+      -- converts an exponential exp(f(x)) to a generalized
+      -- power series
+      (ups := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem => ups
+      expGenUPS(ups.%series,posCheck?,atanFlag)
+
+    expGenUPS(ups,posCheck?,atanFlag) ==
+      -- computes the exponential of a generalized power series.
+      -- If the series has order zero and the constant term a0 of the
+      -- series involves x, the function tries to expand exp(a0) as
+      -- a power series.
+      (deg := order(ups,1)) < 0 =>
+        stateProblem("exp","essential singularity")
+      deg > 0 => [exp ups]
+      lc := coefficient(ups,0); xOpList := opsInvolvingX lc
+      not opInOpList?("log" :: SY,xOpList) => [exp ups]
+      -- try to fix exp(lc) if necessary
+      expCoef :=
+        normalize(exp lc,x)$ElementaryFunctionStructurePackage(R,FE)
+      opInOpList?("log" :: SY,opsInvolvingX expCoef) =>
+        stateProblem("exp","logs in constant coefficient")
+      result := exprToGenUPS(expCoef,posCheck?,atanFlag)
+      result case %problem => result
+      [(result.%series) * exp(ups - monomial(lc,0))]
+
+    atancotToGenUPS(fe,arg,posCheck?,atanFlag,plusMinus) ==
+      -- converts atan(f(x)) to a generalized power series
+      (result := exprToGenUPS(arg,posCheck?,atanFlag)) case %problem =>
+        trouble := result.%problem
+        trouble.prob = "essential singularity" => [monomial(fe,0)$UPS]
+        result
+      ups := result.%series; coef := coefficient(ups,0)
+      -- series involves complex numbers
+      (ord := order(ups,0)) = 0 and coef * coef = -1 =>
+        y := differentiate(ups)/(1 + ups*ups)
+        yCoef := coefficient(y,-1)
+        [monomial(log yCoef,0) + integrate(y - monomial(yCoef,-1)$UPS)]
+      cc : FE :=
+        ord < 0 =>
+          atanFlag = "complex" =>
+            return stateProblem("atan","essential singularity")
+          (rn := ratIfCan(ord :: FE)) case "failed" =>
+            -- this condition usually won't occur because exponents will
+            -- be integers or rational numbers
+            return stateProblem("atan","branch problem")
+          if (atanFlag = "real: two sides") and (odd? numer(rn :: RN)) then
+            -- expansions to the left and right of zero have different
+            -- constant coefficients
+            return stateProblem("atan","branch problem")
+          lc := coefficient(ups,ord)
+          (signum := sign(lc)$SIGNEF) case "failed" =>
+            -- can't determine sign
+            atanFlag = "just do it" =>
+              plusMinus = 1 => pi()/(2 :: FE)
+              0
+            posNegPi2 := signOfExpression(lc) * pi()/(2 :: FE)
+            plusMinus = 1 => posNegPi2
+            pi()/(2 :: FE) - posNegPi2
+            --return stateProblem("atan","branch problem")
+          left? : B := atanFlag = "real: left side"; n := signum :: Integer
+          (left? and n = 1) or (not left? and n = -1) =>
+            plusMinus = 1 => -pi()/(2 :: FE)
+            pi()
+          plusMinus = 1 => pi()/(2 :: FE)
+          0
+        atan coef
+      [(cc :: UPS) + integrate(differentiate(ups)/(1 + ups*ups))]
+
+    genUPSApplyIfCan(fcn,arg,fcnName,posCheck?,atanFlag) ==
+      -- converts fcn(arg) to a generalized power series
+      (series := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem =>
+        series
+      ups := series.%series
+      (deg := order(ups,1)) < 0 =>
+        stateProblem(fcnName,"essential singularity")
+      deg > 0 => [fcn(ups) :: UPS]
+      lc := coefficient(ups,0); xOpList := opsInvolvingX lc
+      null xOpList => [fcn(ups) :: UPS]
+      opInOpList?("log" :: SY,xOpList) =>
+        stateProblem(fcnName,"logs in constant coefficient")
+      contOnReals? fcnName => [fcn(ups) :: UPS]
+      stateProblem(fcnName,"x in constant coefficient")
+
+    applyBddIfCan(fe,fcn,arg,fcnName,posCheck?,atanFlag) ==
+      -- converts fcn(arg) to a generalized power series, where the
+      -- function fcn is bounded for real values
+      -- if fcn(arg) has an essential singularity as a complex
+      -- function, we return fcn(arg) as a monomial of degree 0
+      (ups := iExprToGenUPS(arg,posCheck?,atanFlag)) case %problem =>
+        trouble := ups.%problem
+        trouble.prob = "essential singularity" => [monomial(fe,0)$UPS]
+        ups
+      (ans := fcn(ups.%series)) case "failed" => [monomial(fe,0)$UPS]
+      [ans :: UPS]
+
+    tranToGenUPS(ker,arg,posCheck?,atanFlag) ==
+      -- converts op(arg) to a power series for certain functions
+      -- op in trig or hyperbolic trig categories
+      -- N.B. when this function is called, 'k2elem' will have been
+      -- applied, so the following functions cannot appear:
+      -- tan, cot, sec, csc, sinh, cosh, tanh, coth, sech, csch
+      is?(ker,"sin" :: SY) =>
+        applyBddIfCan(ker :: FE,sinIfCan,arg,"sin",posCheck?,atanFlag)
+      is?(ker,"cos" :: SY) =>
+        applyBddIfCan(ker :: FE,cosIfCan,arg,"cos",posCheck?,atanFlag)
+      is?(ker,"asin" :: SY) =>
+        genUPSApplyIfCan(asinIfCan,arg,"asin",posCheck?,atanFlag)
+      is?(ker,"acos" :: SY) =>
+        genUPSApplyIfCan(acosIfCan,arg,"acos",posCheck?,atanFlag)
+      is?(ker,"atan" :: SY) =>
+        atancotToGenUPS(ker :: FE,arg,posCheck?,atanFlag,1)
+      is?(ker,"acot" :: SY) =>
+        atancotToGenUPS(ker :: FE,arg,posCheck?,atanFlag,-1)
+      is?(ker,"asec" :: SY) =>
+        genUPSApplyIfCan(asecIfCan,arg,"asec",posCheck?,atanFlag)
+      is?(ker,"acsc" :: SY) =>
+        genUPSApplyIfCan(acscIfCan,arg,"acsc",posCheck?,atanFlag)
+      is?(ker,"asinh" :: SY) =>
+        genUPSApplyIfCan(asinhIfCan,arg,"asinh",posCheck?,atanFlag)
+      is?(ker,"acosh" :: SY) =>
+        genUPSApplyIfCan(acoshIfCan,arg,"acosh",posCheck?,atanFlag)
+      is?(ker,"atanh" :: SY) =>
+        genUPSApplyIfCan(atanhIfCan,arg,"atanh",posCheck?,atanFlag)
+      is?(ker,"acoth" :: SY) =>
+        genUPSApplyIfCan(acothIfCan,arg,"acoth",posCheck?,atanFlag)
+      is?(ker,"asech" :: SY) =>
+        genUPSApplyIfCan(asechIfCan,arg,"asech",posCheck?,atanFlag)
+      is?(ker,"acsch" :: SY) =>
+        genUPSApplyIfCan(acschIfCan,arg,"acsch",posCheck?,atanFlag)
+      stateProblem(string name ker,"unknown kernel")
+
+    powToGenUPS(args,posCheck?,atanFlag) ==
+      -- converts a power f(x) ** g(x) to a generalized power series
+      (logBase := logToGenUPS(first args,posCheck?,atanFlag)) case %problem =>
+        logBase
+      expon := iExprToGenUPS(second args,posCheck?,atanFlag)
+      expon case %problem => expon
+      expGenUPS((expon.%series) * (logBase.%series),posCheck?,atanFlag)
+
+@
+\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>>
+
+<<package FS2UPS FunctionSpaceToUnivariatePowerSeries>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}