axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / LIMITPS.spad

)abbrev package LIMITPS PowerSeriesLimitPackage
++ Author: Clifton J. Williamson
++ Date Created: 21 March 1989
++ Date Last Updated: 30 March 1994
++ 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) : SIG == CODE where
  R  : Join(GcdDomain,OrderedSet,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)

  SIG ==> 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))}.

  CODE ==> 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(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
      is?(ker,"fresnelS" :: Symbol) => true
      is?(ker,"fresnelC" :: Symbol) => true
      error "finiteValueAtInfinity? true, but unknown value at infinity"

    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)]
      # expKers ^= 1 => "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))]
      #xkers ^= 1 => "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
            is?(ker,"fresnelS" :: Symbol) => -sqrt(pi()/(8::FE))
            is?(ker,"fresnelC" :: Symbol) => -sqrt(pi()/(8::FE))
            error "finiteValueAtInfinity? true, but unknown value at infinity"
            -- 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)