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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra limitps.spad}
\author{Clifton J. Williamson, Manuel Bronstein}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package LIMITPS PowerSeriesLimitPackage}
<<package LIMITPS PowerSeriesLimitPackage>>=
)abbrev package LIMITPS PowerSeriesLimitPackage
++ Author: Clifton J. Williamson
++ Date Created: 21 March 1989
++ Date Last Updated: 30 March 1994
++ Basic Operations:
++ Related Domains: UnivariateLaurentSeries(FE,x,a),
++   UnivariatePuiseuxSeries(FE,x,a),ExponentialExpansion(R,FE,x,a)
++ Also See:
++ AMS Classifications:
++ Keywords: limit, functional expression, power series
++ Examples:
++ References:
++ Description:
++   PowerSeriesLimitPackage implements limits of expressions
++   in one or more variables as one of the variables approaches a
++   limiting value.  Included are two-sided limits, left- and right-
++   hand limits, and limits at plus or minus infinity.
PowerSeriesLimitPackage(R,FE): Exports == Implementation where
  R  : Join(GcdDomain,RetractableTo Integer,_
            LinearlyExplicitRingOver Integer)
  FE : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
            FunctionSpace R)
  Z       ==> Integer
  RN      ==> Fraction Integer
  RF      ==> Fraction Polynomial R
  OFE     ==> OrderedCompletion FE
  OPF     ==> OnePointCompletion FE
  SY      ==> Symbol
  EQ      ==> Equation
  LF      ==> LiouvillianFunction
  UTS     ==> UnivariateTaylorSeries
  ULS     ==> UnivariateLaurentSeries
  UPXS    ==> UnivariatePuiseuxSeries
  EFULS   ==> ElementaryFunctionsUnivariateLaurentSeries
  EFUPXS  ==> ElementaryFunctionsUnivariatePuiseuxSeries
  FS2UPS  ==> FunctionSpaceToUnivariatePowerSeries
  FS2EXPXP ==> FunctionSpaceToExponentialExpansion
  Problem ==> Record(func:String,prob:String)
  RESULT  ==> Union(OFE,"failed")
  TwoSide ==> Record(leftHandLimit:RESULT,rightHandLimit:RESULT)
  U       ==> Union(OFE,TwoSide,"failed")
  SIGNEF  ==> ElementaryFunctionSign(R,FE)

  Exports ==> with

    limit: (FE,EQ OFE) -> U
      ++ limit(f(x),x = a) computes the real limit \spad{lim(x -> a,f(x))}.

    complexLimit: (FE,EQ OPF) -> Union(OPF, "failed")
      ++ complexLimit(f(x),x = a) computes the complex limit
      ++ \spad{lim(x -> a,f(x))}.

    limit: (FE,EQ FE,String) -> RESULT
      ++ limit(f(x),x=a,"left") computes the left hand real limit
      ++ \spad{lim(x -> a-,f(x))};
      ++ \spad{limit(f(x),x=a,"right")} computes the right hand real limit
      ++ \spad{lim(x -> a+,f(x))}.

  Implementation ==> add
    import ToolsForSign(R)
    import ElementaryFunctionStructurePackage(R,FE)

    zeroFE:FE := 0
    anyRootsOrAtrigs?   : FE -> Boolean
    complLimit  : (FE,SY) -> Union(OPF,"failed")
    okProblem?  : (String,String) -> Boolean
    realLimit   : (FE,SY) -> U
    xxpLimit    : (FE,SY) -> RESULT
    limitPlus   : (FE,SY) -> RESULT
    localsubst  : (FE,Kernel FE,Z,FE) -> FE
    locallimit  : (FE,SY,OFE) -> U
    locallimitcomplex : (FE,SY,OPF) -> Union(OPF,"failed")
    poleLimit:(RN,FE,SY) -> U
    poleLimitPlus:(RN,FE,SY) -> RESULT

    noX?: (FE,SY) -> Boolean
    noX?(fcn,x) == not member?(x,variables fcn)

    constant?: FE -> Boolean
    constant? fcn == empty? variables fcn

    firstNonLogPtr: (FE,SY) -> List Kernel FE
    firstNonLogPtr(fcn,x) ==
      -- returns a pointer to the first element of kernels(fcn) which
      -- has 'x' as a variable, which is not a logarithm, and which is
      -- not simply 'x'
      list := kernels fcn
      while not empty? list repeat
        ker := first list
        not is?(ker,"log" :: Symbol) and member?(x,variables(ker::FE)) _
               and not(x = name operator (ker)) =>
          return list
        list := rest list
      empty()

    finiteValueAtInfinity?: Kernel FE -> Boolean
    finiteValueAtInfinity? ker ==
      is?(ker,"erf" :: Symbol) => true
      is?(ker,"sech" :: Symbol) => true
      is?(ker,"csch" :: Symbol) => true
      is?(ker,"tanh" :: Symbol) => true
      is?(ker,"coth" :: Symbol) => true
      is?(ker,"atan" :: Symbol) => true
      is?(ker,"acot" :: Symbol) => true
      is?(ker,"asec" :: Symbol) => true
      is?(ker,"acsc" :: Symbol) => true
      is?(ker,"acsch" :: Symbol) => true
      is?(ker,"acoth" :: Symbol) => true
      false

    knownValueAtInfinity?: Kernel FE -> Boolean
    knownValueAtInfinity? ker ==
      is?(ker,"exp" :: Symbol) => true
      is?(ker,"sinh" :: Symbol) => true
      is?(ker,"cosh" :: Symbol) => true
      false

    leftOrRight: (FE,SY,FE) -> SingleInteger
    leftOrRight(fcn,x,limVal) ==
      -- function is called when limitPlus(fcn,x) = limVal
      -- determines whether the limiting value is approached
      -- from the left or from the right
      (value := limitPlus(inv(fcn - limVal),x)) case "failed" => 0
      (inf := whatInfinity(val := value :: OFE)) = 0 =>
         error "limit package: internal error"
      inf

    specialLimit1: (FE,SY) -> RESULT
    specialLimitKernel: (Kernel FE,SY) -> RESULT
    specialLimitNormalize: (FE,SY) -> RESULT
    specialLimit: (FE, SY) -> RESULT

    specialLimit(fcn, x) ==
      xkers := [k for k in kernels fcn | member?(x,variables(k::FE))]
      #xkers = 1 => specialLimit1(fcn,x)
      num := numerator fcn
      den := denominator fcn
      for k in xkers repeat
        (fval := limitPlus(k::FE,x)) case "failed" =>
            return specialLimitNormalize(fcn,x)
        whatInfinity(val := fval::OFE) ~= 0 =>
            return specialLimitNormalize(fcn,x)
        (valu := retractIfCan(val)@Union(FE,"failed")) case "failed" =>
            return specialLimitNormalize(fcn,x)
        finVal := valu :: FE
        num := eval(num, k, finVal)
        den := eval(den, k, finVal)
        den = 0 => return specialLimitNormalize(fcn,x)
      (num/den) :: OFE :: RESULT

    specialLimitNormalize(fcn,x) == -- tries to normalize result first
      nfcn := normalize(fcn)
      fcn ~= nfcn => limitPlus(nfcn,x)
      xkers := [k for k in tower fcn | member?(x,variables(k::FE))]
      # xkers ~= 2 => "failed"
      expKers := [k for k in xkers | is?(k, "exp" :: Symbol)]
      not one?(# expKers) => "failed"
    -- fcn is a rational function of x and exp(g(x)) for some rational function g
      expKer := first expKers
      (fval := limitPlus(expKer::FE,x)) case "failed" => "failed"
      vv := new()$SY; eq : EQ FE := equation(expKer :: FE,vv :: FE)
      cc := eval(fcn,eq)
      expKerLim := fval :: OFE
        -- following test for "failed" is needed due to compiler bug
        -- limVal case OFE generates EQCAR(limVal, 1) which fails on atom "failed"
      (limVal := locallimit(cc,vv,expKerLim)) case "failed" => "failed"
      limVal case OFE =>
         limm := limVal :: OFE
         (lim := retractIfCan(limm)@Union(FE,"failed")) case "failed" =>
               "failed" -- need special handling for directions at infinity
         limitPlus(lim, x)
      "failed"

    -- limit of expression having only 1 kernel involving x
    specialLimit1(fcn,x) ==
      -- find the first interesting kernel in tower(fcn)
      xkers := [k for k in kernels fcn | member?(x,variables(k::FE))]
      not one?(#xkers) => "failed"
      ker := first xkers
      vv := new()$SY; eq : EQ FE := equation(ker :: FE,vv :: FE)
      cc := eval(fcn,eq)
      member?(x,variables cc) => "failed"
      (lim := specialLimitKernel(ker, x)) case "failed" => lim
      argLim : OFE := lim :: OFE
      (limVal := locallimit(cc,vv,argLim)) case "failed" => "failed"
      limVal case OFE => limVal :: OFE
      "failed"

    -- limit of single kernel involving x
    specialLimitKernel(ker,x) ==
      is?(ker,"log" :: Symbol) =>
          args := argument ker
          empty? args => "failed" -- error "No argument"
          not empty? rest args => "failed" -- error "Too many arugments"
          arg := first args
          -- compute limit(x -> 0+,arg)
          (limm := limitPlus(arg,x)) case "failed" => "failed"
          lim := limm :: OFE
          (inf := whatInfinity lim) = -1 => "failed"
          argLim : OFE :=
            -- log(+infinity) = +infinity
            inf = 1 => lim
            -- now 'lim' must be finite
            (li := retractIfCan(lim)@Union(FE,"failed") :: FE) = 0 =>
              -- log(0) = -infinity
              leftOrRight(arg,x,0) = 1 => minusInfinity()
              return "failed"
            log(li) :: OFE
      -- kernel should be a function of one argument f(arg)
      args := argument(ker)
      empty? args => "failed"  -- error "No argument"
      not empty? rest args => "failed" -- error "Too many arugments"
      arg := first args
      -- compute limit(x -> 0+,arg)
      (limm := limitPlus(arg,x)) case "failed" => "failed"
      lim := limm :: OFE
      f := elt(operator ker,(var := new()$SY) :: FE)
      -- compute limit(x -> 0+,f(arg))
      -- case where 'lim' is finite
      (inf := whatInfinity lim) = 0 =>
         is?(ker,"erf" :: Symbol) => erf(retract(lim)@FE)$LF(R,FE) :: OFE
         (kerValue := locallimit(f,var,lim)) case "failed" => "failed"
         kerValue case OFE => kerValue :: OFE
         "failed"
      -- case where 'lim' is plus infinity
      inf = 1 =>
        finiteValueAtInfinity? ker =>
          val : FE :=
            is?(ker,"erf" :: Symbol) => 1
            is?(ker,"sech" :: Symbol) => 0
            is?(ker,"csch" :: Symbol) => 0
            is?(ker,"tanh" :: Symbol) => 0
            is?(ker,"coth" :: Symbol) => 0
            is?(ker,"atan" :: Symbol) => pi()/(2 :: FE)
            is?(ker,"acot" :: Symbol) => 0
            is?(ker,"asec" :: Symbol) => pi()/(2 :: FE)
            is?(ker,"acsc" :: Symbol) => 0
            is?(ker,"acsch" :: Symbol) => 0
            -- ker must be acoth
            0
          val :: OFE
        knownValueAtInfinity? ker =>
          lim -- limit(exp, cosh, sinh ,x=inf) = inf
        "failed"
      -- case where 'lim' is minus infinity
      finiteValueAtInfinity? ker =>
        val : FE :=
          is?(ker,"erf" :: Symbol) => -1
          is?(ker,"sech" :: Symbol) => 0
          is?(ker,"csch" :: Symbol) => 0
          is?(ker,"tanh" :: Symbol) => 0
          is?(ker,"coth" :: Symbol) => 0
          is?(ker,"atan" :: Symbol) => -pi()/(2 :: FE)
          is?(ker,"acot" :: Symbol) => pi()
          is?(ker,"asec" :: Symbol) => -pi()/(2 :: FE)
          is?(ker,"acsc" :: Symbol) => -pi()
          is?(ker,"acsch" :: Symbol) => 0
          -- ker must be acoth
          0
        val :: OFE
      knownValueAtInfinity? ker =>
        is?(ker,"exp" :: Symbol) => (0@FE) :: OFE
        is?(ker,"sinh" :: Symbol) => lim
        is?(ker,"cosh" :: Symbol) => plusInfinity()
        "failed"
      "failed"

    logOnlyLimit: (FE,SY) -> RESULT
    logOnlyLimit(coef,x) ==
      -- this function is called when the 'constant' coefficient involves
      -- the variable 'x'. Its purpose is to compute a right hand limit
      -- of an expression involving log x. Here log x is replaced by -1/v,
      -- where v is a new variable. If the new expression no longer involves
      -- x, then take the right hand limit as v -> 0+
      vv := new()$SY
      eq : EQ FE := equation(log(x :: FE),-inv(vv :: FE))
      member?(x,variables(cc := eval(coef,eq))) => "failed"
      limitPlus(cc,vv)

    locallimit(fcn,x,a) ==
      -- Here 'fcn' is a function f(x) = f(x,...) in 'x' and possibly
      -- other variables, and 'a' is a limiting value.  The function
      -- computes lim(x -> a,f(x)).
      xK := retract(x::FE)@Kernel(FE)
      (n := whatInfinity a) = 0 =>
        realLimit(localsubst(fcn,xK,1,retract(a)@FE),x)
      (u := limitPlus(eval(fcn,xK,n * inv(xK::FE)),x))
                                                case "failed" => "failed"
      u::OFE

    localsubst(fcn, k, n, a) ==
      a = 0 and n = 1 => fcn
      eval(fcn,k,n * (k::FE) + a)

    locallimitcomplex(fcn,x,a) ==
      xK := retract(x::FE)@Kernel(FE)
      (g := retractIfCan(a)@Union(FE,"failed")) case FE =>
        complLimit(localsubst(fcn,xK,1,g::FE),x)
      complLimit(eval(fcn,xK,inv(xK::FE)),x)

    limit(fcn,eq,str) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "limit:left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      xK := retract(x::FE)@Kernel(FE)
      limitPlus(localsubst(fcn,xK,direction str,a),x)

    anyRootsOrAtrigs? fcn ==
      -- determines if 'fcn' has any kernels which are roots
      -- or if 'fcn' has any kernels which are inverse trig functions
      -- which could produce series expansions with fractional exponents
      for kernel in tower fcn repeat
        is?(kernel,"nthRoot" :: Symbol) => return true
        is?(kernel,"asin" :: Symbol) => return true
        is?(kernel,"acos" :: Symbol) => return true
        is?(kernel,"asec" :: Symbol) => return true
        is?(kernel,"acsc" :: Symbol) => return true
      false

    complLimit(fcn,x) ==
      -- computes lim(x -> 0,fcn) using a Puiseux expansion of fcn,
      -- if fcn is an expression involving roots, and using a Laurent
      -- expansion of fcn otherwise
      lim : FE :=
        anyRootsOrAtrigs? fcn =>
          ppack := FS2UPS(R,FE,RN,_
              UPXS(FE,x,zeroFE),EFUPXS(FE,ULS(FE,x,zeroFE),UPXS(FE,x,zeroFE),_
              EFULS(FE,UTS(FE,x,zeroFE),ULS(FE,x,zeroFE))),x)
          pseries := exprToUPS(fcn,false,"complex")$ppack
          pseries case %problem => return "failed"
          if pole?(upxs := pseries.%series) then upxs := map(normalize,upxs)
          pole? upxs => return infinity()
          coefficient(upxs,0)
        lpack := FS2UPS(R,FE,Z,ULS(FE,x,zeroFE),_
                 EFULS(FE,UTS(FE,x,zeroFE),ULS(FE,x,zeroFE)),x)
        lseries := exprToUPS(fcn,false,"complex")$lpack
        lseries case %problem => return "failed"
        if pole?(uls := lseries.%series) then uls := map(normalize,uls)
        pole? uls => return infinity()
        coefficient(uls,0)
      -- can the following happen?
      member?(x,variables lim) =>
        member?(x,variables(answer := normalize lim)) =>
          error "limit: can't evaluate limit"
        answer :: OPF
      lim :: FE :: OPF

    okProblem?(function,problem) ==
      (function = "log") or (function = "nth root") =>
        (problem = "series of non-zero order") or _
               (problem = "negative leading coefficient")
      (function = "atan") => problem = "branch problem"
      (function = "erf") => problem = "unknown kernel"
      problem = "essential singularity"

    poleLimit(order,coef,x) ==
      -- compute limit for function with pole
      not member?(x,variables coef) =>
        (s := sign(coef)$SIGNEF) case Integer =>
          rtLim := (s :: Integer) * plusInfinity()
          even? numer order => rtLim
          even? denom order => ["failed",rtLim]$TwoSide
          [-rtLim,rtLim]$TwoSide
        -- infinite limit, but cannot determine sign
        "failed"
      error "limit: can't evaluate limit"

    poleLimitPlus(order,coef,x) ==
      -- compute right hand limit for function with pole
      not member?(x,variables coef) =>
        (s := sign(coef)$SIGNEF) case Integer =>
          (s :: Integer) * plusInfinity()
        -- infinite limit, but cannot determine sign
        "failed"
      (clim := specialLimit(coef,x)) case "failed" => "failed"
      zero? (lim := clim :: OFE) =>
        -- in this event, we need to determine if the limit of
        -- the coef is 0+ or 0-
        (cclim := specialLimit(inv coef,x)) case "failed" => "failed"
        ss := whatInfinity(cclim :: OFE) :: Z
        zero? ss =>
          error "limit: internal error"
        ss * plusInfinity()
      t := whatInfinity(lim :: OFE) :: Z
      zero? t =>
        (tt := sign(coef)$SIGNEF) case Integer =>
          (tt :: Integer) * plusInfinity()
        -- infinite limit, but cannot determine sign
        "failed"
      t * plusInfinity()

    realLimit(fcn,x) ==
      -- computes lim(x -> 0,fcn) using a Puiseux expansion of fcn,
      -- if fcn is an expression involving roots, and using a Laurent
      -- expansion of fcn otherwise
      lim : Union(FE,"failed") :=
        anyRootsOrAtrigs? fcn =>
          ppack := FS2UPS(R,FE,RN,_
               UPXS(FE,x,zeroFE),EFUPXS(FE,ULS(FE,x,zeroFE),UPXS(FE,x,zeroFE),_
               EFULS(FE,UTS(FE,x,zeroFE),ULS(FE,x,zeroFE))),x)
          pseries := exprToUPS(fcn,true,"real: two sides")$ppack
          pseries case %problem =>
            trouble := pseries.%problem
            function := trouble.func; problem := trouble.prob
            okProblem?(function,problem) =>
              left :=
                xK : Kernel FE := kernel x
                fcn0 := eval(fcn,xK,-(xK :: FE))
                limitPlus(fcn0,x)
              right := limitPlus(fcn,x)
              (left case "failed") and (right case "failed") =>
                return "failed"
              if (left case OFE) and (right case OFE) then
                (left :: OFE) = (right :: OFE) => return (left :: OFE)
              return([left,right]$TwoSide)
            return "failed"
          if pole?(upxs := pseries.%series) then upxs := map(normalize,upxs)
          pole? upxs =>
            cp := coefficient(upxs,ordp := order upxs)
            return poleLimit(ordp,cp,x)
          coefficient(upxs,0)
        lpack := FS2UPS(R,FE,Z,ULS(FE,x,zeroFE),_
                 EFULS(FE,UTS(FE,x,zeroFE),ULS(FE,x,zeroFE)),x)
        lseries := exprToUPS(fcn,true,"real: two sides")$lpack
        lseries case %problem =>
          trouble := lseries.%problem
          function := trouble.func; problem := trouble.prob
          okProblem?(function,problem) =>
            left :=
              xK : Kernel FE := kernel x
              fcn0 := eval(fcn,xK,-(xK :: FE))
              limitPlus(fcn0,x)
            right := limitPlus(fcn,x)
            (left case "failed") and (right case "failed") =>
              return "failed"
            if (left case OFE) and (right case OFE) then
              (left :: OFE) = (right :: OFE) => return (left :: OFE)
            return([left,right]$TwoSide)
          return "failed"
        if pole?(uls := lseries.%series) then uls := map(normalize,uls)
        pole? uls =>
          cl := coefficient(uls,ordl := order uls)
          return poleLimit(ordl :: RN,cl,x)
        coefficient(uls,0)
      lim case "failed" => "failed"
      member?(x,variables(lim :: FE)) =>
        member?(x,variables(answer := normalize(lim :: FE))) =>
          error "limit: can't evaluate limit"
        answer :: OFE
      lim :: FE :: OFE

    xxpLimit(fcn,x) ==
      -- computes lim(x -> 0+,fcn) using an exponential expansion of fcn
      xpack := FS2EXPXP(R,FE,x,zeroFE)
      xxp := exprToXXP(fcn,true)$xpack
      xxp case %problem => "failed"
      limitPlus(xxp.%expansion)

    limitPlus(fcn,x) ==
      -- computes lim(x -> 0+,fcn) using a generalized Puiseux expansion
      -- of fcn, if fcn is an expression involving roots, and using a
      -- generalized Laurent expansion of fcn otherwise
      lim : Union(FE,"failed") :=
        anyRootsOrAtrigs? fcn =>
          ppack := FS2UPS(R,FE,RN,_
               UPXS(FE,x,zeroFE),EFUPXS(FE,ULS(FE,x,zeroFE),UPXS(FE,x,zeroFE),_
               EFULS(FE,UTS(FE,x,zeroFE),ULS(FE,x,zeroFE))),x)
          pseries := exprToGenUPS(fcn,true,"real: right side")$ppack
          pseries case %problem =>
            trouble := pseries.%problem
            ff := trouble.func; pp := trouble.prob
            (pp = "negative leading coefficient") => return "failed"
            "failed"
          -- pseries case %problem => return "failed"
          if pole?(upxs := pseries.%series) then upxs := map(normalize,upxs)
          pole? upxs =>
            cp := coefficient(upxs,ordp := order upxs)
            return poleLimitPlus(ordp,cp,x)
          coefficient(upxs,0)
        lpack := FS2UPS(R,FE,Z,ULS(FE,x,zeroFE),_
                 EFULS(FE,UTS(FE,x,zeroFE),ULS(FE,x,zeroFE)),x)
        lseries := exprToGenUPS(fcn,true,"real: right side")$lpack
        lseries case %problem =>
          trouble := lseries.%problem
          ff := trouble.func; pp := trouble.prob
          (pp = "negative leading coefficient") => return "failed"
          "failed"
        -- lseries case %problem => return "failed"
        if pole?(uls := lseries.%series) then uls := map(normalize,uls)
        pole? uls =>
          cl := coefficient(uls,ordl := order uls)
          return poleLimitPlus(ordl :: RN,cl,x)
        coefficient(uls,0)
      lim case "failed" =>
        (xLim := xxpLimit(fcn,x)) case "failed" => specialLimit(fcn,x)
        xLim
      member?(x,variables(lim :: FE)) =>
        member?(x,variables(answer := normalize(lim :: FE))) =>
          (xLim := xxpLimit(answer,x)) case "failed" => specialLimit(answer,x)
          xLim
        answer :: OFE
      lim :: FE :: OFE

    limit(fcn:FE,eq:EQ OFE) ==
      (f := retractIfCan(lhs eq)@Union(FE,"failed")) case "failed" =>
        error "limit:left hand side must be a variable"
      (xx := retractIfCan(f)@Union(SY,"failed")) case "failed" =>
        error "limit:left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      locallimit(fcn,x,a)

    complexLimit(fcn:FE,eq:EQ OPF) ==
      (f := retractIfCan(lhs eq)@Union(FE,"failed")) case "failed" =>
        error "limit:left hand side must be a variable"
      (xx := retractIfCan(f)@Union(SY,"failed")) case "failed" =>
        error "limit:left hand side must be a variable"
      x := xx :: SY; a := rhs eq
      locallimitcomplex(fcn,x,a)

@
\section{package SIGNEF ElementaryFunctionSign}
<<package SIGNEF ElementaryFunctionSign>>=
)abbrev package SIGNEF ElementaryFunctionSign
++ Author: Manuel Bronstein
++ Date Created: 25 Aug 1989
++ Date Last Updated: 4 May 1992
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: elementary function, sign
++ Examples:
++ References:
++ Description:
++   This package provides functions to determine the sign of an
++   elementary function around a point or infinity.
ElementaryFunctionSign(R,F): Exports == Implementation where
  R : Join(IntegralDomain,RetractableTo Integer,_
           LinearlyExplicitRingOver Integer,GcdDomain)
  F : Join(AlgebraicallyClosedField,TranscendentalFunctionCategory,_
            FunctionSpace R)

  N  ==> NonNegativeInteger
  Z  ==> Integer
  SY ==> Symbol
  RF ==> Fraction Polynomial R
  ORF ==> OrderedCompletion RF
  OFE ==> OrderedCompletion F
  K  ==> Kernel F
  P  ==> SparseMultivariatePolynomial(R, K)
  U  ==> Union(Z, "failed")
  FS2 ==> FunctionSpaceFunctions2

  Exports ==> with
    sign: F -> U
      ++ sign(f) returns the sign of f if it is constant everywhere.
    sign: (F, SY, OFE) -> U
      ++ sign(f, x, a) returns the sign of f as x nears \spad{a}, from both
      ++ sides if \spad{a} is finite.
    sign: (F, SY, F, String) -> U
      ++ sign(f, x, a, s) returns the sign of f as x nears \spad{a} from below
      ++ if s is "left", or above if s is "right".

  Implementation ==> add
    macro POSIT == 'positive
    macro NEGAT == 'negative
    import ToolsForSign R
    import RationalFunctionSign(R)
    import PowerSeriesLimitPackage(R, F)
    import TrigonometricManipulations(R, F)

    smpsign : P -> U
    sqfrSign: P -> U
    termSign: P -> U
    kerSign : K -> U
    listSign: (List P,Z) -> U
    insign  : (F,SY,OFE, N) -> U
    psign   : (F,SY,F,String, N) -> U
    ofesign : OFE -> U
    overRF  : OFE -> Union(ORF, "failed")

    sign(f, x, a) ==
      not real? f => "failed"
      insign(f, x, a, 0)

    sign(f, x, a, st) ==
      not real? f => "failed"
      psign(f, x, a, st, 0)

    sign f ==
      not real? f => "failed"
      (u := retractIfCan(f)@Union(RF,"failed")) case RF => sign(u::RF)
      (un := smpsign numer f) case Z and (ud := smpsign denom f) case Z =>
        un::Z * ud::Z
      --abort if there are any variables
      not empty? variables f => "failed"
      -- abort in the presence of algebraic numbers
      member?('rootOf,map(name,operators f)$ListFunctions2(BasicOperator,Symbol)) => "failed"
      -- In the last resort try interval evaluation where feasible.
      if R has ConvertibleTo Float then
        import Interval(Float)
        import Expression(Interval Float)
        mapfun : (R -> Interval(Float)) := interval(convert(#1)$R)
        f2 : Expression(Interval Float) := map(mapfun,f)$FS2(R,F,Interval(Float),Expression(Interval Float))
        r : Union(Interval(Float),"failed") := retractIfCan f2
        if r case "failed" then  return "failed"
        negative? r => return(-1)
        positive? r => return 1
        zero? r => return 0
        "failed"
      "failed"

    overRF a ==
      (n := whatInfinity a) = 0 =>
        (u := retractIfCan(retract(a)@F)@Union(RF,"failed")) _
               case "failed" => "failed"
        u::RF::ORF
      n * plusInfinity()$ORF

    ofesign a ==
      (n := whatInfinity a) ~= 0 => convert(n)@Z
      sign(retract(a)@F)

    insign(f, x, a, m) ==
      m > 10 => "failed"                 -- avoid infinite loops for now
      (uf := retractIfCan(f)@Union(RF,"failed")) case RF and
                   (ua := overRF a) case ORF => sign(uf::RF, x, ua::ORF)
      eq : Equation OFE := equation(x :: F :: OFE,a)
      (u := limit(f,eq)) case "failed" => "failed"
      u case OFE =>
        (n := whatInfinity(u::OFE)) ~= 0 => convert(n)@Z
        (v := retract(u::OFE)@F) = 0 =>
          (s := insign(differentiate(f, x), x, a, m + 1)) case "failed"
                                                             => "failed"
          - s::Z * n
        sign v
      (u.leftHandLimit case "failed") or
         (u.rightHandLimit case "failed") => "failed"
      (ul := ofesign(u.leftHandLimit::OFE))  case "failed" => "failed"
      (ur := ofesign(u.rightHandLimit::OFE)) case "failed" => "failed"
      (ul::Z) = (ur::Z) => ul
      "failed"

    psign(f, x, a, st, m) ==
      m > 10 => "failed"                 -- avoid infinite loops for now
      f = 0 => 0
      (uf := retractIfCan(f)@Union(RF,"failed")) case RF and
           (ua := retractIfCan(a)@Union(RF,"failed")) case RF =>
            sign(uf::RF, x, ua::RF, st)
      eq : Equation F := equation(x :: F,a)
      u := limit(f,eq,st)
      u  case "failed" => "failed"
      (n := whatInfinity(u@OFE)) ~= 0 => convert(n)@Z
      v := retract(u@OFE)@F
      v = 0 =>
        s := psign(differentiate(f,x),x,a,st,m + 1)
        s case "failed"=> "failed"
        direction(st) * s@Z
      sign v

    smpsign p ==
      (r := retractIfCan(p)@Union(R,"failed")) case R => sign(r::R)
      (u := sign(retract(unit(s := squareFree p))@R)) case "failed" =>
        "failed"
      ans := u::Z
      for term in factorList s | odd?(term.xpnt) repeat
        (u := sqfrSign(term.fctr)) case "failed" => return "failed"
        ans := ans * u::Z
      ans

    sqfrSign p ==
      (u := termSign first(l := monomials p)) case "failed" => "failed"
      listSign(rest l, u::Z)

    listSign(l, s) ==
      for term in l repeat
        (u := termSign term) case "failed" => return "failed"
        not(s = u::Z) => return "failed"
      s

    termSign term ==
      (us := sign leadingCoefficient term) case "failed" => "failed"
      for var in (lv := variables term) repeat
        odd? degree(term, var) =>
          empty? rest lv and (vs := kerSign first lv) case Z =>
                                                   return(us::Z * vs::Z)
          return "failed"
      us::Z

    kerSign k ==
      has?(op := operator k, NEGAT) => -1
      has?(op, POSIT) or is?(op,  'pi) or is?(op,'exp) or
                           is?(op,'cosh) or is?(op,'sech) => 1
      empty?(arg := argument k) => "failed"
      (s := sign first arg) case "failed" =>
        is?(op,"nthRoot" :: SY) =>
          even?(retract(second arg)@Z) => 1
          "failed"
        "failed"
      is?(op,"log" :: SY) =>
        negative?(s::Z) => "failed"
        sign(first arg - 1)
      is?(op,"tanh" :: SY) or is?(op,"sinh" :: SY) or
                     is?(op,"csch" :: SY) or is?(op,"coth" :: SY) => s
      is?(op,"nthRoot" :: SY) =>
        even?(retract(second arg)@Z) =>
          negative?(s::Z) => "failed"
          s
        s
      "failed"

@
\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 LIMITPS PowerSeriesLimitPackage>>
<<package SIGNEF ElementaryFunctionSign>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}