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

\documentclass{article}
\usepackage{open-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,RetractableTo Integer,_
               LinearlyExplicitRingOver Integer)
  FE    : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
               FunctionSpace R) with coerce: Expon -> %
  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

    powerToUPS(fcn,n,posCheck?,atanFlag) ==
      -- converts an integral power to a power series
      (b := exprToUPS(fcn,posCheck?,atanFlag)) case %problem => b
      positive? n => [(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 operator 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 :=
        negative? ord =>
          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 operator 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
      positive? n => [(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 operator 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) =>
          negative?(rat :: RN) => 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.
      negative?(deg := order(ups,1)) =>
        stateProblem("exp","essential singularity")
      positive? deg => [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 :=
        negative? ord =>
          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
      negative?(deg := order(ups,1)) =>
        stateProblem(fcnName,"essential singularity")
      positive? deg => [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 operator 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.
--Copyright (C) 2007-2009, Gabriel Dos Reis.
--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}