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

)abbrev domain IR IntegrationResult
++ Author: Barry Trager, Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 12 August 1992
++ Description:
++ The result of a transcendental integration.
++ If a function f has an elementary integral g, then g can be written
++ in the form \spad{g = h + c1 log(u1) + c2 log(u2) + ... + cn log(un)}
++ where h, which is in the same field than f, is called the rational
++ part of the integral, and \spad{c1 log(u1) + ... cn log(un)} is called the
++ logarithmic part of the integral. This domain manipulates integrals
++ represented in that form, by keeping both parts separately. The logs
++ are not explicitly computed.

IntegrationResult(F) : SIG == CODE where
  F : Field

  O   ==> OutputForm
  B   ==> Boolean
  Z   ==> Integer
  Q   ==> Fraction Integer
  SE  ==> Symbol
  UP  ==> SparseUnivariatePolynomial F
  LOG ==> Record(scalar:Q, coeff:UP, logand:UP)
  NE  ==> Record(integrand:F, intvar:F)

  SIG ==> (Module Q, RetractableTo F) with

    mkAnswer : (F, List LOG, List NE) -> %
      ++ mkAnswer(r,l,ne) creates an integration result from
      ++ a rational part r, a logarithmic part l, and a non-elementary part ne.

    ratpart : % -> F
      ++ ratpart(ir) returns the rational part of an integration result

    logpart : % -> List LOG
      ++ logpart(ir) returns the logarithmic part of an integration result

    notelem : % -> List NE
      ++ notelem(ir) returns the non-elementary part of an integration result

    elem? : % -> B
      ++ elem?(ir) tests if an integration result is elementary over F?

    integral : (F, F) -> %
      ++ integral(f,x) returns the formal integral of f with respect to x

    differentiate : (%, F -> F) -> F
      ++ differentiate(ir,D) differentiates ir with respect to the derivation D

    if F has PartialDifferentialRing(SE) then

      differentiate : (%, Symbol) -> F
        ++ differentiate(ir,x) differentiates ir with respect to x

    if F has RetractableTo Symbol then

      integral : (F, Symbol) -> %
        ++ integral(f,x) returns the formal integral of f with respect to x

  CODE ==> add

    Rep := Record(ratp: F, logp: List LOG, nelem: List NE)

    timelog : (Q, LOG) -> LOG
    timene  : (Q, NE)  -> NE
    LOG2O   : LOG      -> O
    NE2O    : NE       -> O
    Q2F     : Q        -> F
    nesimp  : List NE  -> List NE
    neselect: (List NE, F) -> F
    pLogDeriv: (LOG, F -> F) -> F
    pNeDeriv : (NE,  F -> F) -> F

    alpha:O := new()$Symbol :: O

    - u               == (-1$Z) * u

    0                 == mkAnswer(0, empty(), empty())

    coerce(x:F):%     == mkAnswer(x, empty(), empty())

    ratpart u         == u.ratp

    logpart u         == u.logp

    notelem u         == u.nelem

    elem? u           == empty? notelem u

    mkAnswer(x, l, n) == [x, l, nesimp n]

    timelog(r, lg)    == [r * lg.scalar, lg.coeff, lg.logand]

    integral(f:F,x:F) == (zero? f => 0; mkAnswer(0, empty(), [[f, x]]))

    timene(r, ne)     == [Q2F(r) * ne.integrand, ne.intvar]

    n:Z * u:%         == (n::Q) * u

    Q2F r             == numer(r)::F / denom(r)::F

    neselect(l, x)    == _+/[ne.integrand for ne in l | ne.intvar = x]

    if F has RetractableTo Symbol then
      integral(f:F, x:Symbol):% == integral(f, x::F)

    LOG2O rec ==
      (degree rec.coeff) = 1 =>
        -- deg 1 minimal poly doesn't get sigma
        lastc := - coefficient(rec.coeff, 0) / coefficient(rec.coeff, 1)
        lg    := (rec.logand) lastc
        logandp := prefix("log"::Symbol::O, [lg::O])
        (cc := Q2F(rec.scalar) * lastc) = 1 => logandp
        cc = -1 => - logandp
        cc::O * logandp
      coeffp:O := (outputForm(rec.coeff, alpha) = 0::Z::O)@O
      logandp :=
           alpha * prefix("log"::Symbol::O, [outputForm(rec.logand, alpha)])
      if (cc := Q2F(rec.scalar)) ^= 1 then
        logandp := cc::O * logandp
      sum(logandp, coeffp)

    nesimp l ==
      [[u,x] for x in removeDuplicates_!([ne.intvar for ne in l]$List(F))
                                           | (u := neselect(l, x)) ^= 0]

    if (F has LiouvillianFunctionCategory) _
     and (F has RetractableTo Symbol) then

      retractIfCan u ==
        empty? logpart u =>
          ratpart u +
             _+/[integral(ne.integrand, retract(ne.intvar)@Symbol)$F
                for ne in notelem u]
        "failed"

    else

      retractIfCan u ==
        elem? u and empty? logpart u => ratpart u
        "failed"

    r:Q * u:% ==
      r = 0 => 0
      mkAnswer(Q2F(r) * ratpart u, map(x1+->timelog(r, x1), logpart u),
                                       map(x2+->timene(r, x2), notelem u))

    -- Initial attempt, quick and dirty, no simplification
    u + v ==
      mkAnswer(ratpart u + ratpart v, concat(logpart u, logpart v),
                                    nesimp concat(notelem u, notelem v))

    if F has PartialDifferentialRing(Symbol) then
      differentiate(u:%, x:Symbol):F == 
        differentiate(u, x1+->differentiate(x1, x))

    differentiate(u:%, derivation:F -> F):F ==
      derivation ratpart u +
          _+/[pLogDeriv(log, derivation) for log in logpart u]
               + _+/[pNeDeriv(ne, derivation) for ne in notelem u]

    pNeDeriv(ne, derivation) ==
      (derivation(ne.intvar) = 1) => ne.integrand
      zero? derivation(ne.integrand) => 0
      error "pNeDeriv: cannot differentiate not elementary part into F"

    pLogDeriv(log, derivation) ==
      map(derivation, log.coeff) ^= 0 =>
        error "pLogDeriv: can only handle logs with constant coefficients"
      ((n := degree(log.coeff)) = 1) =>
        c   := - (leadingCoefficient reductum log.coeff)
                                        / (leadingCoefficient log.coeff)
        ans := (log.logand) c
        Q2F(log.scalar) * c * derivation(ans) / ans
      numlog := map(derivation, log.logand)
      diflog := extendedEuclidean(log.logand, log.coeff,
                                    numlog)::Record(coef1:UP, coef2:UP)
      algans := diflog.coef1
      ans:F := 0
      for i in 0..(n-1) repeat
        algans := algans * monomial(1, 1) rem log.coeff
        ans := ans + coefficient(algans, i)
      Q2F(log.scalar) * ans

    coerce(u:%):O ==
      (r := retractIfCan u) case F => r::F::O
      l := reverse_! [LOG2O f for f in logpart u]$List(O)
      if ratpart u ^= 0 then l := concat(ratpart(u)::O, l)
      if not elem? u then l := concat([NE2O f for f in notelem u], l)
      null l => 0::O
      reduce("+", l)

    NE2O ne ==
      int((ne.integrand)::O * hconcat ["d"::Symbol::O, (ne.intvar)::O])