defintef.spad

)abbrev package DEFINTEF ElementaryFunctionDefiniteIntegration
++ Definite integration of elementary functions.
++ Author: Manuel Bronstein
++ Date Created: 14 April 1992
++ Description:
++   \spadtype{ElementaryFunctionDefiniteIntegration}
++   provides functions to compute definite
++   integrals of elementary functions.
ElementaryFunctionDefiniteIntegration(R, F) : Exports == Implementation where
  R : Join(PolynomialFactorizationExplicit, Comparable, CharacteristicZero,
           RetractableTo Integer, LinearlyExplicitOver Integer)
  F : Join(TranscendentalFunctionCategory, PrimitiveFunctionCategory,
           AlgebraicallyClosedFunctionSpace R)

  B   ==> Boolean
  SY  ==> Symbol
  Z   ==> Integer
  P   ==> SparseMultivariatePolynomial(R, K)
  K   ==> Kernel F
  UP  ==> SparseUnivariatePolynomial F
  OFE ==> OrderedCompletion F
  U   ==> Union(f1:OFE, f2:List OFE, fail:"failed", pole:"potentialPole")

  Exports ==> with
    integrate : (F, SegmentBinding OFE) -> U
      ++ integrate(f, x = a..b) returns the integral of
      ++ \spad{f(x)dx} from a to b.
      ++ Error: if f has a pole for x between a and b.
    integrate : (F, SegmentBinding OFE, String) -> U
      ++ integrate(f, x = a..b, "noPole") returns the
      ++ integral of \spad{f(x)dx} from a to b.
      ++ If it is not possible to check whether f has a pole for x
      ++ between a and b (because of parameters), then this function
      ++ will assume that f has no such pole.
      ++ Error: if f has a pole for x between a and b or
      ++ if the last argument is not "noPole".
    innerint : (F, SY, OFE, OFE, B) -> U
      ++ innerint(f, x, a, b, ignore?) should be local but conditional

  Implementation ==> add
    import from ElementaryFunctionSign(R, F)
    import from DefiniteIntegrationTools(R, F)
    import from FunctionSpaceIntegration(R, F)

    polyIfCan   : (P, K) -> Union(UP, "failed")
    int         : (F, SY, OFE, OFE, B) -> U
    nopole      : (F, SY, K, OFE, OFE) -> U
    checkFor0   : (P, K, OFE, OFE) -> Union(B, "failed")
    checkSMP    : (P, SY, K, OFE, OFE) -> Union(B, "failed")
    checkForPole : (F, SY, K, OFE, OFE) -> Union(B, "failed")
    posit       : (F, SY, K, OFE, OFE) -> Union(B, "failed")
    negat       : (F, SY, K, OFE, OFE) -> Union(B, "failed")
    moreThan    : (OFE, Fraction Z) -> Union(B, "failed")

    if R has Join(ConvertibleTo Pattern Integer, PatternMatchable Integer)
      and F has SpecialFunctionCategory then
        import from PatternMatchIntegration(R, F)

        innerint(f, x, a, b, ignor?) ==
          ((u := int(f, x, a, b, ignor?)) case f1) or (u case f2)
            or ((v := pmintegrate(f, x, a, b)) case "failed") => u
          [v@F::OFE]

    else
      innerint(f, x, a, b, ignor?) == int(f, x, a, b, ignor?)

    integrate(f : F, s : SegmentBinding OFE) ==
      innerint(f, variable s, low(segment(s)), high(segment(s)), false)

    integrate(f : F, s : SegmentBinding OFE, str : String) ==
      innerint(f, variable s, low(segment(s)), high(segment(s)), ignore? str)

    int(f, x, a, b, ignor?) ==
      a = b => [0$F::OFE]
      k := kernel(x)@Kernel(F)
      (z := checkForPole(f, x, k, a, b)) case "failed" =>
        ignor? => nopole(f, x, k, a, b)
        ["potentialPole"]
      z@B => error "integrate: pole in path of integration"
      nopole(f, x, k, a, b)

    checkForPole(f, x, k, a, b) ==
      ((u := checkFor0(d := denom f, k, a, b)) case "failed") or (u@B) => u
      ((u := checkSMP(d, x, k, a, b)) case "failed") or (u@B) => u
      checkSMP(numer f, x, k, a, b)

-- true if p has a zero between a and b exclusive
    checkFor0(p, x, a, b) ==
      (u := polyIfCan(p, x)) case UP => checkForZero(u@UP, a, b, false)
      (v := isTimes p) case List(P) =>
         for t in v@List(P) repeat
           ((w := checkFor0(t, x, a, b)) case "failed") or (w@B) => return w
         false
      (z := isExpt p) case "failed" => "failed"
      k := z.var
-- functions with no real zeros
      -- FIXME: maybe also is?(k, 'nthRoot)
      is?(k, 'exp) or is?(k, 'acot) or is?(k, 'cosh) => false
-- special case for log
      is?(k, 'log) =>
        (w := moreThan(b, 1)) case "failed" or not(w@B) => w
        moreThan(-a, -1)
      "failed"

-- returns true if a > b, false if a < b, "failed" if can't decide
    moreThan(a, b) ==
      (r := retractIfCan(a)@Union(F, "failed")) case "failed" =>  -- infinite
        whatInfinity(a) > 0
      (u := retractIfCan(r@F)@Union(Fraction Z, "failed")) case "failed" =>
        "failed"
      u@Fraction(Z) > b

-- true if p has a pole between a and b
    checkSMP(p, x, k, a, b) ==
      (u := polyIfCan(p, k)) case UP => false
      (v := isTimes p) case List(P) =>
         for t in v@List(P) repeat
           ((w := checkSMP(t, x, k, a, b)) case "failed") or (w@B) => return w
         false
      (v := isPlus p) case List(P) =>
         n := 0              -- number of summand having a pole
         for t in v@List(P) repeat
           (w := checkSMP(t, x, k, a, b)) case "failed" => return w
           if w@B then n := n + 1
         zero? n => false    -- no summand has a pole
         (n = 1) => true      -- only one summand has a pole
         "failed"            -- at least 2 summands have a pole
      (z := isExpt(p)) case "failed" => "failed"
      kk := z.var
      -- nullary operators have no poles
      nullary? operator kk => false
      ak := argument kk
      f := first ak
      -- functions which are defined over all the reals:
      is?(kk, 'exp) or is?(kk, 'sin) or is?(kk, 'cos)
        or is?(kk, 'sinh) or is?(kk, 'cosh) or is?(kk, 'tanh)
          or is?(kk, 'sech) or is?(kk, 'atan) or is?(kk, 'acot)
            or is?(kk, 'asinh) or is?(kk, 'erf)  or is?(kk, 'erfi)
              or is?(kk, 'fresnelC)  or is?(kk, 'fresnelS)
                or is?(kk, 'Si) or is?(kk, 'Shi) => checkForPole(f, x, k, a, b)
      -- functions which are defined on (-1, +1):
      is?(kk, 'asin) or is?(kk, 'acos) or is?(kk, 'atanh) =>
        ((w := checkForPole(f, x, k, a, b)) case "failed") or (w@B) => w
        ((w := posit(f - 1, x, k, a, b)) case "failed") or (w@B) => w
        negat(f + 1, x, k, a, b)
      -- functions which are defined on (+1, +infty):
      is?(kk, 'acosh) =>
        ((w := checkForPole(f, x, k, a, b)) case "failed") or (w@B) => w
        negat(f - 1, x, k, a, b)
      -- functions which are defined on (0, +infty):
      is?(kk, 'log) or is?(kk, 'nthRoot) or is?(kk, 'Ei) or is?(kk, 'Ci)
        or is?(kk, 'Chi) or is?(kk, 'dilog) =>
          ((w := checkForPole(f, x, k, a, b)) case "failed") or (w@B) => w
          negat(f, x, k, a, b)
      is?(kk, 'Gamma) and #ak = 2 and D(f, x) = 0 =>
          f2 := ak(2)
          ((w := checkForPole(f2, x, k, a, b)) case "failed") or (w@B) => w
          negat(f2, x, k, a, b)
      -- defined on (-infty, 1)
      is?(kk, 'polylog) and D(f, x) = 0 =>
          f2 := ak(2)
          ((w := checkForPole(f2, x, k, a, b)) case "failed") or (w@B) => w
          posit(f - 1, x, k, a, b)
      "failed"

-- returns true if it is certain that f takes at least one strictly positive
-- value for x in (a, b), false if it is certain that f takes no strictly
-- positive value in (a,b), "failed" otherwise
-- f must be known to have no poles in (a, b)
    posit(f, x, k, a, b) ==
      z :=
          (r := retractIfCan(a)@Union(F, "failed")) case "failed" =>
            sign(f, x, a)
          sign(f, x, r@F, "right")
      (b1 := z case Z) and z::Z > 0 => true
      z :=
          (r := retractIfCan(b)@Union(F, "failed")) case "failed" =>
              sign(f, x, b)
          sign(f, x, r@F, "left")
      (b2 := z case Z) and z::Z > 0 => true
      b1 and b2 =>
        ((w := checkFor0(numer f, k, a, b)) case "failed") or (w@B) => "failed"
        false
      "failed"

-- returns true if it is certain that f takes at least one strictly negative
-- value for x in (a, b), false if it is certain that f takes no strictly
-- negative value in (a,b), "failed" otherwise
-- f must be known to have no poles in (a, b)
    negat(f, x, k, a, b) ==
      z :=
        (r := retractIfCan(a)@Union(F, "failed")) case "failed" => sign(f, x, a)
        sign(f, x, r@F, "right")
      (b1 := z case Z) and z::Z < 0 => true
      z :=
        (r := retractIfCan(b)@Union(F, "failed")) case "failed" => sign(f, x, b)
        sign(f, x, r@F, "left")
      (b2 := z case Z) and z::Z < 0 => true
      b1 and b2 =>
        ((w := checkFor0(numer f, k, a, b)) case "failed") or (w@B) => "failed"
        false
      "failed"

-- returns a UP if p is only a poly w.r.t. the kernel x
    polyIfCan(p, x) ==
      q := univariate(p, x)
      ans : UP := 0
      while q ~= 0 repeat
        member?(x, tower(c := leadingCoefficient(q)::F)) => return "failed"
        ans := ans + monomial(c, degree q)
        q := reductum q
      ans

-- integrate f for x between a and b assuming that f has no pole in between
    nopole(f, x, k, a, b) ==
      (u := integrate(f, x)) case F =>
        (v := computeInt(k, u::F, a, b, false)) case "failed" => ["failed"]
        [v@OFE]
      ans := empty()$List(OFE)
      for g in u::List(F) repeat
        (v := computeInt(k, g, a, b, false)) case "failed" => return ["failed"]
        ans := concat!(ans, [v@OFE])
      [ans]

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