--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.
 
)abbrev package TOOLSIGN ToolsForSign
++ Tools for the sign finding utilities
++ Author: Manuel Bronstein
++ Date Created: 25 August 1989
++ Date Last Updated: 26 November 1991
++ Description: Tools for the sign finding utilities.
ToolsForSign(R:Ring): with
    sign     : R      -> Union(Integer, "failed")
	++ sign(r) \undocumented
    nonQsign : R      -> Union(Integer, "failed")
	++ nonQsign(r) \undocumented
    direction: String -> Integer
	++ direction(s) \undocumented
 == add
 
    if R is AlgebraicNumber then
      nonQsign r ==
        sign((r pretend AlgebraicNumber)::Expression(
                  Integer))$ElementaryFunctionSign(Integer, Expression Integer)
    else
      nonQsign r == "failed"
 
    if R has RetractableTo Fraction Integer then
      sign r ==
        (u := retractIfCan(r)@Union(Fraction Integer, "failed"))
          case Fraction(Integer) => sign(u::Fraction Integer)
        nonQsign r
 
    else
      if R has RetractableTo Integer then
        sign r ==
          (u := retractIfCan(r)@Union(Integer, "failed"))
            case "failed" => "failed"
          sign(u::Integer)
 
      else
        sign r ==
          zero? r => 0
          one? r => 1
          r = -1 => -1
          "failed"
 
    direction st ==
      st = "right" => 1
      st = "left" => -1
      error "Unknown option"

)abbrev package INPSIGN InnerPolySign
--%% InnerPolySign
++ Author: Manuel Bronstein
++ Date Created: 23 Aug 1989
++ Date Last Updated: 19 Feb 1990
++ Description:
++ Find the sign of a polynomial around a point or infinity.
InnerPolySign(R, UP): Exports == Implementation where
  R : Ring
  UP: UnivariatePolynomialCategory R
 
  U ==> Union(Integer, "failed")
 
  Exports ==> with
    signAround: (UP,    Integer, R -> U) -> U
	++ signAround(u,i,f) \undocumented
    signAround: (UP, R, Integer, R -> U) -> U
	++ signAround(u,r,i,f) \undocumented
    signAround: (UP, R,          R -> U) -> U
	++ signAround(u,r,f) \undocumented
 
  Implementation ==> add
    signAround(p:UP, x:R, rsign:R -> U) ==
      (ur := signAround(p, x,  1, rsign)) case "failed" => "failed"
      (ul := signAround(p, x, -1, rsign)) case "failed" => "failed"
      (ur::Integer) = (ul::Integer) => ur
      "failed"
 
    signAround(p, x, dir, rsign) ==
      zero? p => 0
      zero?(r := p x) =>
        (u := signAround(differentiate p, x, dir, rsign)) case "failed"
          => "failed"
        dir * u::Integer
      rsign r
 
    signAround(p:UP, dir:Integer, rsign:R -> U) ==
      zero? p => 0
      (u := rsign leadingCoefficient p) case "failed" => "failed"
      positive? dir or (even? degree p) => u::Integer
      - (u::Integer)

)abbrev package SIGNRF RationalFunctionSign
--%% RationalFunctionSign
++ Author: Manuel Bronstein
++ Date Created: 23 August 1989
++ Date Last Updated: 26 November 1991
++ Description:
++ Find the sign of a rational function around a point or infinity.
RationalFunctionSign(R:GcdDomain): Exports == Implementation where
  SE  ==> Symbol
  P   ==> Polynomial R
  RF  ==> Fraction P
  ORF ==> OrderedCompletion RF
  UP  ==> SparseUnivariatePolynomial RF
  U   ==> Union(Integer, "failed")
  SGN ==> ToolsForSign(R)
 
  Exports ==> with
    sign: RF -> U
      ++ sign f returns the sign of f if it is constant everywhere.
    sign: (RF, SE, ORF) -> U
      ++ sign(f, x, a) returns the sign of f as x approaches \spad{a},
      ++ from both sides if \spad{a} is finite.
    sign: (RF, SE, RF, String) -> U
      ++ sign(f, x, a, s) returns the sign of f as x nears \spad{a} from
      ++ the left (below) if s is the string \spad{"left"},
      ++ or from the right (above) if s is the string \spad{"right"}.
 
  Implementation ==> add
    import SGN
    import InnerPolySign(RF, UP)
    import PolynomialCategoryQuotientFunctions(IndexedExponents SE,
                                                      SE, R, P, RF)
 
    psign     : P -> U
    sqfrSign  : P -> U
    termSign  : P -> U
    listSign  : (List P, Integer) -> U
    finiteSign: (Fraction UP, RF) -> U
 
    sign f ==
      (un := psign numer f) case "failed" => "failed"
      (ud := psign denom f) case "failed" => "failed"
      (un::Integer) * (ud::Integer)
 
    finiteSign(g, a) ==
      (ud := signAround(denom g, a, sign$%)) case "failed" => "failed"
      (un := signAround(numer g, a, sign$%)) case "failed" => "failed"
      (un::Integer) * (ud::Integer)
 
    sign(f, x, a) ==
      g := univariate(f, x)
      zero?(n := whatInfinity a) => finiteSign(g, retract a)
      (ud := signAround(denom g, n, sign$%)) case "failed" => "failed"
      (un := signAround(numer g, n, sign$%)) case "failed" => "failed"
      (un::Integer) * (ud::Integer)
 
    sign(f, x, a, st) ==
      (ud := signAround(denom(g := univariate(f, x)), a,
                    d := direction st, sign$%)) case "failed" => "failed"
      (un := signAround(numer g, a, d, sign$%)) case "failed" => "failed"
      (un::Integer) * (ud::Integer)
 
    psign p ==
      (r := retractIfCan(p)@Union(R, "failed")) case R => sign(r::R)$SGN
      (u := sign(retract(unit(s := squareFree p))@R)$SGN) case "failed" =>
        "failed"
      ans := u::Integer
      for term in factors s | odd?(term.exponent) repeat
        (u := sqfrSign(term.factor)) case "failed" => return "failed"
        ans := ans * (u::Integer)
      ans
 
    sqfrSign p ==
      (u := termSign first(l := monomials p)) case "failed" => "failed"
      listSign(rest l, u::Integer)
 
    listSign(l, s) ==
      for term in l repeat
        (u := termSign term) case "failed" => return "failed"
        u::Integer ~= s => return "failed"
      s
 
    termSign term ==
      for var in variables term repeat
        odd? degree(term, var) => return "failed"
      sign(leadingCoefficient term)$SGN

)abbrev package LIMITRF RationalFunctionLimitPackage
++ Computation of limits for rational functions
++ Author: Manuel Bronstein
++ Date Created: 4 October 1989
++ Date Last Updated: 26 November 1991
++ Description: Computation of limits for rational functions.
++ Keywords: limit, rational function.
RationalFunctionLimitPackage(R:GcdDomain):Exports==Implementation where
  Z       ==> Integer
  P       ==> Polynomial R
  RF      ==> Fraction P
  EQ      ==> Equation
  ORF     ==> OrderedCompletion RF
  OPF     ==> OnePointCompletion RF
  UP      ==> SparseUnivariatePolynomial RF
  SE      ==> Symbol
  QF      ==> Fraction SparseUnivariatePolynomial RF
  Result  ==> Union(ORF, "failed")
  TwoSide ==> Record(leftHandLimit:Result, rightHandLimit:Result)
  U       ==> Union(ORF, TwoSide, "failed")
  RFSGN   ==> RationalFunctionSign(R)
 
  Exports ==> with
-- The following are the one we really want, but the interpreter cannot
-- handle them...
--  limit: (RF,EQ ORF) -> U
--  ++ limit(f(x),x,a) computes the real two-sided limit lim(x -> a,f(x))
 
--  complexLimit: (RF,EQ OPF) -> OPF
--  ++ complexLimit(f(x),x,a) computes the complex limit lim(x -> a,f(x))
 
-- ... so we replace them by the following 4:
    limit: (RF,EQ OrderedCompletion P) -> U
      ++ limit(f(x),x = a) computes the real two-sided limit
      ++ of f as its argument x approaches \spad{a}.
    limit: (RF,EQ RF) -> U
      ++ limit(f(x),x = a) computes the real two-sided limit
      ++ of f as its argument x approaches \spad{a}.
    complexLimit: (RF,EQ OnePointCompletion P) -> OPF
      ++ \spad{complexLimit(f(x),x = a)} computes the complex limit
      ++ of \spad{f} as its argument x approaches \spad{a}.
    complexLimit: (RF,EQ RF) -> OPF
      ++ complexLimit(f(x),x = a) computes the complex limit
      ++ of f as its argument x approaches \spad{a}.
    limit: (RF,EQ RF,String) -> Result
      ++ limit(f(x),x,a,"left") computes the real limit
      ++ of f as its argument x approaches \spad{a} from the left;
      ++ limit(f(x),x,a,"right") computes the corresponding limit as x
      ++ approaches \spad{a} from the right.
 
  Implementation ==> add
    import ToolsForSign R
    import InnerPolySign(RF, UP)
    import RFSGN
    import PolynomialCategoryQuotientFunctions(IndexedExponents SE,
                                                      SE, R, P, RF)
 
    finiteComplexLimit: (QF, RF) -> OPF
    finiteLimit       : (QF, RF) -> U
    fLimit            : (Z, UP, RF, Z) -> Result
 
-- These 2 should be exported, see comment above
    locallimit       : (RF, SE, ORF) -> U
    locallimitcomplex: (RF, SE, OPF) -> OPF
 
    limit(f:RF,eq:EQ RF) ==
      (xx := retractIfCan(lhs eq)@Union(SE,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx :: SE; a := rhs eq
      locallimit(f,x,a::ORF)
 
    complexLimit(f:RF,eq:EQ RF) ==
      (xx := retractIfCan(lhs eq)@Union(SE,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx :: SE; a := rhs eq
      locallimitcomplex(f,x,a::OPF)
 
    limit(f:RF,eq:EQ OrderedCompletion P) ==
      (p := retractIfCan(lhs eq)@Union(P,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      (xx := retractIfCan(p)@Union(SE,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx :: SE
      a := map(#1::RF,rhs eq)$OrderedCompletionFunctions2(P,RF)
      locallimit(f,x,a)
 
    complexLimit(f:RF,eq:EQ OnePointCompletion P) ==
      (p := retractIfCan(lhs eq)@Union(P,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      (xx := retractIfCan(p)@Union(SE,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx :: SE
      a := map(#1::RF,rhs eq)$OnePointCompletionFunctions2(P,RF)
      locallimitcomplex(f,x,a)
 
    fLimit(n, d, a, dir) ==
      (s := signAround(d, a, dir, sign$RFSGN)) case "failed" => "failed"
      n * (s::Z) * plusInfinity()
 
    finiteComplexLimit(f, a) ==
      zero?(n := (numer f) a) => 0
      zero?(d := (denom f) a) => infinity()
      (n / d)::OPF
 
    finiteLimit(f, a) ==
      zero?(n := (numer f) a) => 0
      zero?(d := (denom f) a) =>
        (s := sign(n)$RFSGN) case "failed" => "failed"
        rhsl := fLimit(s::Z, denom f, a, 1)
        lhsl := fLimit(s::Z, denom f, a, -1)
        rhsl case "failed" =>
          lhsl case "failed" => "failed"
          [lhsl, rhsl]
        lhsl case "failed" => [lhsl, rhsl]
        rhsl::ORF = lhsl::ORF => lhsl::ORF
        [lhsl, rhsl]
      (n / d)::ORF
 
    locallimit(f,x,a) ==
      g := univariate(f, x)
      zero?(n := whatInfinity a) => finiteLimit(g, retract a)
      (dn := degree numer g) > (dd := degree denom g) =>
        (sn := signAround(numer g, n, sign$RFSGN)) case "failed" => "failed"
        (sd := signAround(denom g, n, sign$RFSGN)) case "failed" => "failed"
        (sn::Z) * (sd::Z) * plusInfinity()
      dn < dd => 0
      ((leadingCoefficient numer g) / (leadingCoefficient denom g))::ORF
 
    limit(f,eq,st) ==
      (xx := retractIfCan(lhs eq)@Union(SE,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx :: SE; a := rhs eq
      zero?(n := (numer(g := univariate(f, x))) a) => 0
      zero?(d := (denom g) a) =>
        (s := sign(n)$RFSGN) case "failed" => "failed"
        fLimit(s::Z, denom g, a, direction st)
      (n / d)::ORF
 
    locallimitcomplex(f,x,a) ==
      g := univariate(f, x)
      (r := retractIfCan(a)@Union(RF, "failed")) case RF =>
        finiteComplexLimit(g, r::RF)
      (dn := degree numer g) > (dd := degree denom g) => infinity()
      dn < dd => 0
      ((leadingCoefficient numer g) / (leadingCoefficient denom g))::OPF

