sign.spad

)abbrev package TOOLSIGN ToolsForSign
++ Tools for the sign finding utilities
++ Author: Manuel Bronstein
++ Date Created: 25 August 1989
++ 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
          r = 1 => 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
++ 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"
      (dir > 0) or (even? degree p) => u@Integer
      - (u@Integer)

)abbrev package SIGNRF RationalFunctionSign
--%% RationalFunctionSign
++ Author: Manuel Bronstein
++ Date Created: 23 August 1989
++ Description:
++ Find the sign of a rational function around a point or infinity.
RationalFunctionSign(R : GcdDomain) : Exports == Implementation where
  SY  ==> 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, SY, 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, SY, 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 from SGN
    import from InnerPolySign(RF, UP)
    import from PolynomialCategoryQuotientFunctions(IndexedExponents SY,
                                                      SY, 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 factorList 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
++ Description: Computation of limits for rational functions.
++ Keywords: limit, rational function.
-- Exposed, but unused internally
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
  SY      ==> 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 from ToolsForSign R
    import from InnerPolySign(RF, UP)
    import from RFSGN
    import from PolynomialCategoryQuotientFunctions(IndexedExponents SY,
                                                      SY, 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, SY, ORF) -> U
    locallimitcomplex : (RF, SY, OPF) -> OPF

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

    complexLimit(f : RF, eq : EQ RF) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx@SY; 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(SY,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx@SY
      a := map(y +-> y::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(SY,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx@SY
      a := map(y +-> y::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"
      si := n * (s@Z)
      si = 0 => error "fLimit: sign is 0"
      si > 0 => plusInfinity()
      minusInfinity()

    finiteComplexLimit(f, a) ==
      zero?(n := (numer f) a) => 0$RF::OPF
      zero?(d := (denom f) a) => infinity()
      (n / d)::OPF

    finiteLimit(f, a) ==
      zero?(n := (numer f) a) => 0$RF::ORF
      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"
        si := (sn@Z) * (sd@Z)
        si = 0 => error "locallimit: sign is 0"
        si > 0 => plusInfinity()
        minusInfinity()
      dn < dd => 0$RF::ORF
      ((leadingCoefficient numer g) / (leadingCoefficient denom g))::ORF

    limit(f, eq, st) ==
      (xx := retractIfCan(lhs eq)@Union(SY,"failed")) case "failed" =>
        error "limit: left hand side must be a variable"
      x := xx@SY; a := rhs eq
      zero?(n := (numer(g := univariate(f, x))) a) => 0$RF::ORF
      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$RF::OPF
      ((leadingCoefficient numer g) / (leadingCoefficient denom g))::OPF

--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.