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

)abbrev package DFINTTLS DefiniteIntegrationTools
++ Author: Manuel Bronstein
++ Date Created: 15 April 1992
++ Date Last Updated: 24 February 1993
++ Description:
++ \spadtype{DefiniteIntegrationTools} provides common tools used
++ by the definite integration of both rational and elementary functions.

DefiniteIntegrationTools(R, F) : SIG == CODE where
  R : Join(GcdDomain, OrderedSet, RetractableTo Integer,
           LinearlyExplicitRingOver Integer)
  F : Join(TranscendentalFunctionCategory,
           AlgebraicallyClosedFunctionSpace R)

  B   ==> Boolean
  Z   ==> Integer
  Q   ==> Fraction Z
  SE  ==> Symbol
  P   ==> Polynomial R
  RF  ==> Fraction P
  UP  ==> SparseUnivariatePolynomial F
  K   ==> Kernel F
  OFE ==> OrderedCompletion F
  UPZ ==> SparseUnivariatePolynomial Z
  UPQ ==> SparseUnivariatePolynomial Q
  REC ==> Record(left:Q, right:Q)
  REC2==> Record(endpoint:Q, dir:Z)
  U   ==> Union(fin:REC, halfinf:REC2, all:"all", failed:"failed")
  IGNOR ==> "noPole"

  SIG ==> with

    ignore? : String -> B
      ++ ignore?(s) is true if s is the string that tells the integrator
      ++ to assume that the function has no pole in the integration interval.

    computeInt : (K, F, OFE, OFE, B) -> Union(OFE, "failed")
      ++ computeInt(x, g, a, b, eval?) returns the integral of \spad{f} for x
      ++ between a and b, assuming that g is an indefinite integral of
      ++ \spad{f} and \spad{f} has no pole between a and b.
      ++ If \spad{eval?} is true, then \spad{g} can be evaluated safely
      ++ at \spad{a} and \spad{b}, provided that they are finite values.
      ++ Otherwise, limits must be computed.

    checkForZero : (P,  SE, OFE, OFE, B) -> Union(B, "failed")
      ++ checkForZero(p, x, a, b, incl?) is true if p has a zero for x between
      ++ a and b, false otherwise, "failed" if this cannot be determined.
      ++ Check for a and b inclusive if incl? is true, exclusive otherwise.

    checkForZero : (UP, OFE, OFE, B) -> Union(B, "failed")
      ++ checkForZero(p, a, b, incl?) is true if p has a zero between
      ++ a and b, false otherwise, "failed" if this cannot be determined.
      ++ Check for a and b inclusive if incl? is true, exclusive otherwise.

  CODE ==> add

    import RealZeroPackage UPZ
    import InnerPolySign(F, UP)
    import ElementaryFunctionSign(R, F)
    import PowerSeriesLimitPackage(R, F)
    import UnivariatePolynomialCommonDenominator(Z, Q, UPQ)

    mkLogPos    : F -> F
    keeprec?    : (Q, REC) -> B
    negative    : F -> Union(B, "failed")
    mkKerPos    : K -> Union(F, "positive")
    posRoot     : (UP, B) -> Union(B, "failed")
    realRoot    : UP -> Union(B, "failed")
    var         : UP -> Union(Z, "failed")
    maprat      : UP -> Union(UPZ, "failed")
    variation   : (UP, F) -> Union(Z, "failed")
    infeval     : (UP, OFE) -> Union(F, "failed")
    checkHalfAx : (UP, F, Z, B) -> Union(B, "failed")
    findLimit   : (F, K, OFE, String, B) -> Union(OFE, "failed")
    checkBudan  : (UP, OFE, OFE, B) -> Union(B, "failed")
    checkDeriv  : (UP, OFE, OFE) -> Union(B, "failed")
    sameSign    : (UP, OFE, OFE) -> Union(B, "failed")
    intrat      : (OFE, OFE) -> U
    findRealZero: (UPZ, U, B) -> List REC

    variation(p, a)      == var p(monomial(1, 1)$UP - a::UP)

    keeprec?(a, rec)     == (a > rec.right) or (a < rec.left)

    checkHalfAx(p, a, d, incl?) ==
      posRoot(p(d * (monomial(1, 1)$UP - a::UP)), incl?)

    ignore? str ==
      str = IGNOR => true
      error "integrate: last argument must be 'noPole'"

    computeInt(k, f, a, b, eval?) ==
      is?(f, "integral"::SE) => "failed"
      if not eval? then f := mkLogPos f
      ((ib := findLimit(f, k, b, "left", eval?)) case "failed") or
         ((ia := findLimit(f, k, a, "right", eval?)) case "failed") => "failed"
      infinite?(ia::OFE) and (ia::OFE = ib::OFE) => "failed"
      ib::OFE - ia::OFE

    findLimit(f, k, a, dir, eval?) ==
      r := retractIfCan(a)@Union(F, "failed")
      r case F =>
        eval? => mkLogPos(eval(f, k, r::F))::OFE
        (u := limit(f, equation(k::F, r::F), dir)) case OFE => u::OFE
        "failed"
      (u := limit(f, equation(k::F::OFE, a))) case OFE => u::OFE
      "failed"

    mkLogPos f ==
      lk := empty()$List(K)
      lv := empty()$List(F)
      for k in kernels f | is?(k, "log"::SE) repeat
        if (v := mkKerPos k) case F then
          lk := concat(k, lk)
          lv := concat(v::F, lv)
      eval(f, lk, lv)

    mkKerPos k ==
      (u := negative(f := first argument k)) case "failed" =>
                                                     log(f**2) / (2::F)
      u::B => log(-f)
      "positive"

    negative f ==
      (u := sign f) case "failed" => "failed"
      u::Z < 0

    checkForZero(p, x, a, b, incl?) ==
      checkForZero(
        map(s+->s::F, univariate(p, x))_
         $SparseUnivariatePolynomialFunctions2(P, F),
            a, b, incl?)

    checkForZero(q, a, b, incl?) ==
      ground? q => false
      (d := maprat q) case UPZ and not((i := intrat(a, b)) case failed) =>
          not empty? findRealZero(d::UPZ, i, incl?)
      (u := checkBudan(q, a, b, incl?)) case "failed" =>
         incl? => checkDeriv(q, a, b)
         "failed"
      u::B

    maprat p ==
      ans:UPQ := 0
      while p ^= 0 repeat
        (r := retractIfCan(c := leadingCoefficient p)@Union(Q,"failed"))
          case "failed"  => return "failed"
        ans := ans + monomial(r::Q, degree p)
        p   := reductum p
      map(numer,(splitDenominator ans).num
         )$SparseUnivariatePolynomialFunctions2(Q, Z)

    intrat(a, b) ==
      (n := whatInfinity a) ^= 0 =>
        (r := retractIfCan(b)@Union(F,"failed")) case "failed" => ["all"]
        (q := retractIfCan(r::F)@Union(Q, "failed")) case "failed" =>
          ["failed"]
        [[q::Q, n]]
      (q := retractIfCan(retract(a)@F)@Union(Q,"failed")) case "failed"
        => ["failed"]
      (n := whatInfinity b) ^= 0 => [[q::Q, n]]
      (t := retractIfCan(retract(b)@F)@Union(Q,"failed")) case "failed"
        => ["failed"]
      [[q::Q, t::Q]]

    findRealZero(p, i, incl?) ==
      i case fin =>
        l := realZeros(p, r := i.fin)
        incl? => l
        select_!(s+->keeprec?(r.left, s) and keeprec?(r.right, s), l)
      i case all => realZeros p
      i case halfinf =>
        empty?(l := realZeros p) => empty()
        bounds:REC :=
          i.halfinf.dir > 0 => [i.halfinf.endpoint, "max"/[t.right for t in l]]
          ["min"/[t.left for t in l], i.halfinf.endpoint]
        l := [u::REC for t in l | (u := refine(p, t, bounds)) case REC]
        incl? => l
        ep := i.halfinf.endpoint
        select_!(s+->keeprec?(ep, s), l)
      error "findRealZero: should not happpen"

    checkBudan(p, a, b, incl?) ==
      r := retractIfCan(b)@Union(F, "failed")
      (n := whatInfinity a) ^= 0 =>
        r case "failed" => realRoot p
        checkHalfAx(p, r::F, n, incl?)
      (za? := zero? p(aa := retract(a)@F)) and incl? => true
      (n := whatInfinity b) ^= 0 => checkHalfAx(p, aa, n, incl?)
      (zb? := zero? p(bb := r::F)) and incl? => true
      (va := variation(p, aa)) case "failed" or
                   (vb := variation(p, bb)) case "failed" => "failed"
      m:Z := 0
      if za? then m := inc m
      if zb? then m := inc m
      odd?(v := va::Z - vb::Z) =>          -- p has an odd number of roots
        incl? or even? m => true
        (v = 1) => false
        "failed"
      zero? v => false                     -- p has no roots
      (m = 1) => true                     -- p has an even number > 0 of roots
      "failed"

    checkDeriv(p, a, b) ==
      (r := retractIfCan(p)@Union(F, "failed")) case F => zero?(r::F)
      (s := sameSign(p, a, b)) case "failed" => "failed"
      s::B =>                  -- p has the same nonzero sign at a and b
        (u := checkDeriv(differentiate p,a,b)) case "failed" => "failed"
        u::B => "failed"
        false
      true

    realRoot p ==
      (b := posRoot(p, true)) case "failed" => "failed"
      b::B => true
      posRoot(p(p - monomial(1, 1)$UP), true)

    sameSign(p, a, b) ==
      (ea := infeval(p, a)) case "failed" => "failed"
      (eb := infeval(p, b)) case "failed" => "failed"
      (s := sign(ea::F * eb::F)) case "failed" => "failed"
      s::Z > 0

    -- returns true if p has a positive root. Include 0 is incl0? is true
    posRoot(p, incl0?) ==
      (z0? := zero?(coefficient(p, 0))) and incl0? => true
      (v := var p) case "failed" => "failed"
      odd?(v::Z) =>            -- p has an odd number of positive roots
        incl0? or not(z0?) => true
        (v::Z) = 1 => false
        "failed"
      zero?(v::Z) => false     -- p has no positive roots
      z0? => true              -- p has an even number > 0 of positive roots
      "failed"

    infeval(p, a) ==
      zero?(n := whatInfinity a) => p(retract(a)@F)
      (u := signAround(p, n, sign)) case "failed" => "failed"
      u::Z::F

    var q ==
      i:Z := 0
      (lastCoef := negative leadingCoefficient q) case "failed" =>
        "failed"
      while ((q := reductum q) ^= 0) repeat
        (next := negative leadingCoefficient q) case "failed" =>
          return "failed"
        if ((not(lastCoef::B)) and next::B) or
                        ((not(next::B)) and lastCoef::B) then i := i + 1
        lastCoef := next
      i