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

)abbrev package INTPM PatternMatchIntegration
++ Author: Manuel Bronstein
++ Date Created: 5 May 1992
++ Date Last Updated: 27 September 1995
++ Description:
++ \spadtype{PatternMatchIntegration} provides functions that use
++ the pattern matcher to find some indefinite and definite integrals
++ involving special functions and found in the litterature.

PatternMatchIntegration(R, F) : SIG == CODE where
  R : Join(OrderedSet, RetractableTo Integer, GcdDomain,
           LinearlyExplicitRingOver Integer)
  F : Join(AlgebraicallyClosedField, TranscendentalFunctionCategory,
           FunctionSpace R)

  N   ==> NonNegativeInteger
  Z   ==> Integer
  SY  ==> Symbol
  K   ==> Kernel F
  P   ==> SparseMultivariatePolynomial(R, K)
  SUP ==> SparseUnivariatePolynomial F
  PAT ==> Pattern Z
  RES ==> PatternMatchResult(Z, F)
  OFE ==> OrderedCompletion F
  REC ==> Record(which: Z, exponent: F, coeff: F)
  ANS ==> Record(special:F, integrand:F)
  NONE ==> 0
  EI   ==> 1
  ERF  ==> 2
  SI   ==> 3
  CI   ==> 4
  GAM2 ==> 5
  CI0  ==> 6

  SIG ==> with

    splitConstant: (F, SY) -> Record(const:F, nconst:F)
      ++ splitConstant(f, x) returns \spad{[c, g]} such that
      ++ \spad{f = c * g} and \spad{c} does not involve \spad{t}.

    if R has ConvertibleTo Pattern Integer and
       R has PatternMatchable Integer then

         if F has LiouvillianFunctionCategory then

           pmComplexintegrate : (F, SY) -> Union(ANS, "failed")
             ++ pmComplexintegrate(f, x) returns either "failed" or
             ++ \spad{[g,h]} such that
             ++ \spad{integrate(f,x) = g + integrate(h,x)}.
             ++ It only looks for special complex integrals that pmintegrate
             ++ does not return.

           pmintegrate : (F, SY) -> Union(ANS, "failed")
             ++ pmintegrate(f, x) returns either "failed" or \spad{[g,h]} such
             ++ that \spad{integrate(f,x) = g + integrate(h,x)}.

         if F has SpecialFunctionCategory then

           pmintegrate : (F, SY, OFE, OFE) -> Union(F, "failed")
             ++ pmintegrate(f, x = a..b) returns the integral of
             ++ \spad{f(x)dx} from a to b
             ++ if it can be found by the built-in pattern matching rules.

  CODE ==> add

    import PatternMatch(Z, F, F)
    import ElementaryFunctionSign(R, F)
    import FunctionSpaceAssertions(R, F)
    import TrigonometricManipulations(R, F)
    import FunctionSpaceAttachPredicates(R, F, F)

    mkalist   : RES -> AssociationList(SY, F)

    pm := new()$SY
    pmw := new pm
    pmm := new pm
    pms := new pm
    pmc := new pm
    pma := new pm
    pmb := new pm

    c := optional(pmc::F)

    w := suchThat(optional(pmw::F), 
      (x1:F):Boolean +-> empty? variables x1)

    s := suchThat(optional(pms::F), 
      (x1:F):Boolean +-> empty? variables x1 and real? x1)

    m := suchThat(optional(pmm::F),
      (x1:F):Boolean+->(retractIfCan(x1)@Union(Z,"failed") case Z) and x1 >= 0)

    spi := sqrt(pi()$F)

    half := 1::F / 2::F

    mkalist res == construct destruct res

    splitConstant(f, x) ==
      not member?(x, variables f) => [f, 1]
      (retractIfCan(f)@Union(K, "failed")) case K => [1, f]
      (u := isTimes f) case List(F) =>
        cc := nc := 1$F
        for g in u::List(F) repeat
          rec := splitConstant(g, x)
          cc  := cc * rec.const
          nc  := nc * rec.nconst
        [cc, nc]
      (u := isPlus f) case List(F) =>
        rec := splitConstant(first(u::List(F)), x)
        cc  := rec.const
        nc  := rec.nconst
        for g in rest(u::List(F)) repeat
          rec := splitConstant(g, x)
          if rec.nconst = nc then cc := cc + rec.const
          else if rec.nconst = -nc then cc := cc - rec.const
          else return [1, f]
        [cc, nc]
      if (v := isPower f) case Record(val:F, exponent:Z) then
        vv := v::Record(val:F, exponent:Z)
        (vv.exponent ^= 1) =>
          rec := splitConstant(vv.val, x)
          return [rec.const ** vv.exponent, rec.nconst ** vv.exponent]
      error "splitConstant: should not happen"

    if R has ConvertibleTo Pattern Integer and
       R has PatternMatchable Integer then
         if F has LiouvillianFunctionCategory then
           import ElementaryFunctionSign(R, F)

           insqrt     : F -> F
           matchei    : (F, SY) -> REC
           matcherfei : (F, SY, Boolean) -> REC
           matchsici  : (F, SY) -> REC
           matchli    : (F, SY) -> List F
           matchli0   : (F, K, SY) -> List F
           matchdilog : (F, SY) -> List F
           matchdilog0: (F, K, SY, P, F) -> List F
           goodlilog? : (K, P) -> Boolean
           gooddilog? : (K, P, P) -> Boolean

           goodlilog?(k, p) == is?(k, "log"::SY) and (minimumDegree(p, k) = 1)

           gooddilog?(k, p, q) ==
             is?(k, "log"::SY) and (degree(p, k) = 1) and zero? degree(q, k)

           -- matches the integral to a result of the form d*erf(u) or d*ei(u)
           -- returns [case, u, d]
           matcherfei(f, x, comp?) ==
             res0 := new()$RES
             pat := c * exp(pma::F)
             failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               comp? => [NONE, 0,0]
               matchei(f,x)
             l := mkalist res
             da := differentiate(a := l.pma, x)
             d := a * (cc := l.pmc) / da
             zero? differentiate(d, x) => [EI, a, d]
             comp? or (((u := sign a) case Z) and (u::Z) < 0) =>
               d := cc * (sa := insqrt(- a)) / da
               zero? differentiate(d, x) => [ERF, sa, - d * spi]
               [NONE, 0, 0]
             [NONE, 0, 0]

           -- matches the integral to a result of the form d * ei(k * log u)
           -- returns [case, k * log u, d]
           matchei(f, x) ==
             res0 := new()$RES
             a := pma::F
             pat := c * a**w / log a
             failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               [NONE, 0, 0]
             l := mkalist res
             da := differentiate(a := l.pma, x)
             d := (cc := l.pmc) / da
             zero? differentiate(d, x) => [EI, (1 + l.pmw) * log a, d]
             [NONE, 0, 0]

           -- matches the integral to a result of the form d*dilog(u) + int(v),
           -- returns [u,d,v] or []
           matchdilog(f, x) ==
             n := numer f
             df := (d := denom f)::F
             for k in select_!(
              (x1:K):Boolean +-> gooddilog?(x1,n,d),variables n)$List(K) repeat
                 not empty?(l := matchdilog0(f, k, x, n, df)) => return l
             empty()

           -- matches the integral to a result of the form d*dilog(a) + int(v)
           -- where k = log(a)
           -- returns [a,d,v] or []
           matchdilog0(f, k, x, p, q) ==
             zero?(da := differentiate(a := first argument k, x)) => empty()
             a1 := 1 - a
             d := coefficient(univariate(p, k), 1)::F * a1 / (q * da)
             zero? differentiate(d, x) => [a, d, f - d * da * (k::F) / a1]
             empty()

           -- matches the integral to a result of the form d * li(u) + int(v),
           -- returns [u,d,v] or []
           matchli(f, x) ==
             d := denom f
             for k in select_!(
               (x1:K):Boolean+->goodlilog?(x1,d), variables d)$List(K) repeat
                 not empty?(l := matchli0(f, k, x)) => return l
             empty()

           -- matches the integral to a result of the form d * li(a) + int(v)
           -- where k = log(a)
           -- returns [a,d,v] or []
           matchli0(f, k, x) ==
             g := (lg := k::F) * f
             zero?(da := differentiate(a := first argument k, x)) => empty()
             zero? differentiate(d := g / da, x) => [a, d, 0]
             ug := univariate(g, k)
             (u:=retractIfCan(ug)@Union(SUP,"failed")) case "failed" => empty()
             degree(p := u::SUP) > 1 => empty()
             zero? differentiate(d := coefficient(p, 0) / da, x) =>
               [a, d, leadingCoefficient p]
             empty()
           -- matches the integral to a result of the form 
           -- d * Si(u) or d * Ci(u) returns [case, u, d]
           matchsici(f, x) ==
             res0 := new()$RES
             b := pmb::F
             t := tan(a := pma::F)
             patsi := c * t / (patden := b + b * t**2)
             patci := (c - c * t**2) / patden
             patci0 := c / patden
             ci0?:Boolean
             (ci? := failed?(res := patternMatch(f, convert(patsi)@PAT, res0)))
              and (ci0?:=failed?(res:=patternMatch(f,convert(patci)@PAT,res0)))
                and failed?(res := patternMatch(f,convert(patci0)@PAT,res0)) =>
                   [NONE, 0, 0]
             l := mkalist res
             (b := l.pmb) ^= 2 * (a := l.pma) => [NONE, 0, 0]
             db := differentiate(b, x)
             d := (cc := l.pmc) / db
             zero? differentiate(d, x) =>
               ci? =>
                  ci0? => [CI0, b, d / (2::F)]
                  [CI, b, d]
               [SI, b, d / (2::F)]
             [NONE, 0, 0]

           -- returns a simplified sqrt(y)
           insqrt y ==
             rec := froot(y, 2)$PolynomialRoots(IndexedExponents K, K, R, P, F)
             ((rec.exponent) = 1) => rec.coef * rec.radicand
             rec.exponent ^=2 => error "insqrt: hould not happen"
             rec.coef * sqrt(rec.radicand)

           pmintegrate(f, x) ==
             (rc := splitConstant(f, x)).const ^= 1 =>
               (u := pmintegrate(rc.nconst, x)) case "failed" => "failed"
               rec := u::ANS
               [rc.const * rec.special, rc.const * rec.integrand]
             not empty?(l := matchli(f, x)) => [second l * li first l, third l]
             not empty?(l := matchdilog(f, x)) =>
                                            [second l * dilog first l, third l]
             cse := (rec := matcherfei(f, x, false)).which
             cse = EI   => [rec.coeff * Ei(rec.exponent), 0]
             cse = ERF  => [rec.coeff * erf(rec.exponent), 0]
             cse := (rec := matchsici(f, x)).which
             cse = SI => [rec.coeff * Si(rec.exponent), 0]
             cse = CI => [rec.coeff * Ci(rec.exponent), 0]
             cse = CI0 => [rec.coeff * Ci(rec.exponent)
                           + rec.coeff * log(rec.exponent), 0]
             "failed"

           pmComplexintegrate(f, x) ==
             (rc := splitConstant(f, x)).const ^= 1 =>
               (u := pmintegrate(rc.nconst, x)) case "failed" => "failed"
               rec := u::ANS
               [rc.const * rec.special, rc.const * rec.integrand]
             cse := (rec := matcherfei(f, x, true)).which
             cse = ERF  => [rec.coeff * erf(rec.exponent), 0]
             "failed"

         if F has SpecialFunctionCategory then

           match1    : (F, SY, F, F) -> List F
           formula1  : (F, SY, F, F) -> Union(F, "failed")

           -- tries only formula (1) of the Geddes & al, AAECC 1 (1990) paper
           formula1(f, x, t, cc) ==
             empty?(l := match1(f, x, t, cc)) => "failed"
             mw := first l
             zero?(ms := third l) or ((sgs := sign ms) case "failed")=>_
                "failed"
             ((sgz := sign(z := (mw + 1) / ms)) case "failed") or (sgz::Z < 0)
                => "failed"
             mmi := retract(mm := second l)@Z
             sgs * (last l) * ms**(- mmi - 1) *
                eval(differentiate(Gamma(x::F), x, mmi::N), [kernel(x)@K], [z])

           -- returns [w, m, s, c] or []
           -- matches only formula (1) of the Geddes & al, AAECC 1 (1990) paper
           match1(f, x, t, cc) ==
             res0 := new()$RES
             pat := cc * log(t)**m * exp(-t**s)
             not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               l := mkalist res
               [0, l.pmm, l.pms, l.pmc]
             pat := cc * t**w * exp(-t**s)
             not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               l := mkalist res
               [l.pmw, 0, l.pms, l.pmc]
             pat := cc / t**w * exp(-t**s)
             not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               l := mkalist res
               [- l.pmw, 0, l.pms, l.pmc]
             pat := cc * t**w * log(t)**m * exp(-t**s)
             not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               l := mkalist res
               [l.pmw, l.pmm, l.pms, l.pmc]
             pat := cc / t**w * log(t)**m * exp(-t**s)
             not failed?(res := patternMatch(f, convert(pat)@PAT, res0)) =>
               l := mkalist res
               [- l.pmw, l.pmm, l.pms, l.pmc]
             empty()

           pmintegrate(f, x, a, b) ==
             zero? a and ((whatInfinity b) = 1) =>
               formula1(f, x, constant(x::F), 
                 suchThat(c, (x1:F):Boolean +-> freeOf?(x1, x)))
             "failed"