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

)abbrev package IR2F IntegrationResultToFunction
++ Author: Manuel Bronstein
++ Date Created: 4 February 1988
++ Date Last Updated: 9 October 1991
++ Description:
++ Conversion of integration results to top-level expressions
++ This package allows a sum of logs over the roots of a polynomial
++ to be expressed as explicit logarithms and arc tangents, provided
++ that the indexing polynomial can be factored into quadratics.

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

  N   ==> NonNegativeInteger
  Z   ==> Integer
  Q   ==> Fraction Z
  K   ==> Kernel F
  P   ==> SparseMultivariatePolynomial(R, K)
  UP  ==> SparseUnivariatePolynomial F
  IR  ==> IntegrationResult F
  REC ==> Record(ans1:F, ans2:F)
  LOG ==> Record(scalar:Q, coeff:UP, logand:UP)

  SIG ==> with

    split : IR -> IR
       ++ split(u(x) + sum_{P(a)=0} Q(a,x)) returns
       ++ \spad{u(x) + sum_{P1(a)=0} Q(a,x) + ... + sum_{Pn(a)=0} Q(a,x)}
       ++ where P1,...,Pn are the factors of P.

    expand : IR -> List F
       ++ expand(i) returns the list of possible real functions
       ++ corresponding to i.

    complexExpand : IR -> F
       ++ complexExpand(i) returns the expanded complex function
       ++ corresponding to i.

  CODE ==> add

    import AlgebraicManipulations(R, F)
    import ElementaryFunctionSign(R, F)

    IR2F       : IR -> F
    insqrt     : F  -> Record(sqrt:REC, sgn:Z)
    pairsum    : (List F, List F) -> List F
    pairprod   : (F, List F) -> List F
    quadeval   : (UP, F, F, F) -> REC
    linear     : (UP, UP) -> F
    tantrick   : (F, F) -> F
    ilog       : (F, F, List K) -> F
    ilog0      : (F, F, UP, UP, F) -> F
    nlogs      : LOG -> List LOG
    lg2func    : LOG -> List F
    quadratic  : (UP, UP) -> List F
    mkRealFunc : List LOG -> List F
    lg2cfunc   : LOG -> F
    loglist    : (Q, UP, UP) -> List LOG
    cmplex     : (F, UP) -> F
    evenRoots  : F -> List F
    compatible?: (List F, List F) -> Boolean

    cmplex(alpha, p) == alpha * log p alpha

    IR2F i           == retract mkAnswer(ratpart i, empty(), notelem i)

    pairprod(x, l)   == [x * y for y in l]

    evenRoots x ==
      [first argument k for k in tower x |
       is?(k,"nthRoot"::Symbol) and even?(retract(second argument k)@Z)
         and (not empty? variables first argument k)]

    expand i ==
      j := split i
      pairsum([IR2F j], mkRealFunc logpart j)

    split i ==
      mkAnswer(ratpart i,concat [nlogs l for l in logpart i],notelem i)

    complexExpand i ==
      j := split i
      IR2F j + +/[lg.scalar::F * lg2cfunc lg for lg in logpart j]

    -- p = a t^2 + b t + c
    -- Expands sum_{p(t) = 0} t log(lg(t))
    quadratic(p, lg) ==
      zero?(delta := (b := coefficient(p, 1))**2 - 4 *
       (a := coefficient(p,2)) * (p0 := coefficient(p, 0))) =>
         [linear(monomial(1, 1) + (b / a)::UP, lg)]
      e := (q := quadeval(lg, c := - b * (d := inv(2*a)),d, delta)).ans1
      lgp  := c * log(nrm := (e**2 - delta * (f := q.ans2)**2))
      s    := (sqr := insqrt delta).sqrt
      pp := nn := 0$F
      if sqr.sgn >= 0 then
        sqrp := s.ans1 * rootSimp sqrt(s.ans2)
        pp := lgp + d * sqrp * log(((2 * e * f) / nrm) * sqrp
                                          + (e**2 + delta * f**2) / nrm)
      if sqr.sgn <= 0 then
        sqrn := s.ans1 * rootSimp sqrt(-s.ans2)
        nn := lgp + d * sqrn * ilog(e, f * sqrn,
                   setUnion(setUnion(kernels a, kernels b), kernels p0))
      sqr.sgn > 0 => [pp]
      sqr.sgn < 0 => [nn]
      [pp, nn]

    -- returns 2 atan(a/b) or 2 atan(-b/a) whichever looks better
    -- they differ by a constant so it's ok to do it from an IR
    tantrick(a, b) ==
      retractIfCan(a)@Union(Q, "failed") case Q => 2 * atan(-b/a)
      2 * atan(a/b)

    -- transforms i log((a + i b) / (a - i b)) into a sum of real
    -- arc-tangents using Rioboo's algorithm
    -- lk is a list of kernels which are parameters for the integral
    ilog(a, b, lk) ==
      l := setDifference(setUnion(variables numer a, variables numer b),
           setUnion(lk, setUnion(variables denom a, variables denom b)))
      empty? l => tantrick(a, b)
      k := "max"/l
      ilog0(a, b, numer univariate(a, k), numer univariate(b, k), k::F)

    -- transforms i log((a + i b) / (a - i b)) into a sum of real
    -- arc-tangents using Rioboo's algorithm
    -- the arc-tangents will not have k in the denominator
    -- we always keep upa(k) = a  and upb(k) = b
    ilog0(a, b, upa, upb, k) ==
      if degree(upa) < degree(upb) then
        (upa, upb) := (-upb, upa)
        (a, b) := (-b, a)
      zero? degree upb => tantrick(a, b)
      r := extendedEuclidean(upa, upb)
      (g:= retractIfCan(r.generator)@Union(F,"failed")) case "failed" =>
        tantrick(a, b)
      if degree(r.coef1) >= degree upb then
        qr := divide(r.coef1, upb)
        r.coef1 := qr.remainder
        r.coef2 := r.coef2 + qr.quotient * upa
      aa := (r.coef2) k
      bb := -(r.coef1) k
      tantrick(aa * a + bb * b, g::F) + ilog0(aa,bb,r.coef2,-r.coef1,k)

    lg2func lg ==
      zero?(d := degree(p := lg.coeff)) => error "poly has degree 0"
      (d = 1) => [linear(p, lg.logand)]
      d = 2  => quadratic(p, lg.logand)
      odd? d and
        ((r := retractIfCan(reductum p)@Union(F, "failed")) case F) =>
            pairsum([cmplex(alpha := rootSimp zeroOf p, lg.logand)],
                    lg2func [lg.scalar,
                     (p exquo (monomial(1, 1)$UP - alpha::UP))::UP,
                      lg.logand])
      [lg2cfunc lg]

    lg2cfunc lg ==
      +/[cmplex(alpha, lg.logand) for alpha in zerosOf(lg.coeff)]

    mkRealFunc l ==
      ans := empty()$List(F)
      for lg in l repeat
        ans := pairsum(ans, pairprod(lg.scalar::F, lg2func lg))
      ans

    -- returns a log(b)
    linear(p, lg) ==
      alpha := - coefficient(p, 0) / coefficient(p, 1)
      alpha * log lg alpha

    -- returns (c, d) s.t. p(a + b t) = c + d t, where t^2 = delta
    quadeval(p, a, b, delta) ==
      zero? p => [0, 0]
      bi := c := d := 0$F
      ai := 1$F
      v  := vectorise(p, 1 + degree p)
      for i in minIndex v .. maxIndex v repeat
        c    := c + qelt(v, i) * ai
        d    := d + qelt(v, i) * bi
        temp := a * ai + b * bi * delta
        bi   := a * bi + b * ai
        ai   := temp
      [c, d]

    compatible?(lx, ly) ==
      empty? ly => true
      for x in lx repeat
        for y in ly repeat
          ((s := sign(x*y)) case Z) and (s::Z < 0) => return false
      true

    pairsum(lx, ly) ==
      empty? lx => ly
      empty? ly => lx
      l := empty()$List(F)
      for x in lx repeat
        ls := evenRoots x
        if not empty?(ln :=
          [x + y for y in ly | compatible?(ls, evenRoots y)]) then
            l := removeDuplicates concat(l, ln)
      l

    -- returns [[a, b], s] where sqrt(y) = a sqrt(b) and
    -- s = 1 if b > 0, -1 if b < 0, 0 if the sign of b cannot be determined
    insqrt y ==
      rec := froot(y, 2)$PolynomialRoots(IndexedExponents K, K, R, P, F)
      ((rec.exponent) = 1) => [[rec.coef * rec.radicand, 1], 1]
      rec.exponent ^=2 => error "Should not happen"
      [[rec.coef, rec.radicand],
          ((s := sign(rec.radicand)) case "failed" => 0; s::Z)]

    nlogs lg ==
      [[f.exponent * lg.scalar, f.factor, lg.logand] for f in factors
         ffactor(primitivePart(lg.coeff)
                    )$FunctionSpaceUnivariatePolynomialFactor(R, F, UP)]