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

)abbrev package INTPAF PureAlgebraicIntegration
++ Author: Manuel Bronstein
++ Date Created: 27 May 1988
++ Date Last Updated: 24 June 1994
++ Description:
++ Integration of pure algebraic functions;
++ This package provides functions for integration, limited integration,
++ extended integration and the risch differential equation for
++ pure algebraic integrands;

PureAlgebraicIntegration(R, F, L) : SIG == CODE where
  R : Join(GcdDomain,RetractableTo Integer,OrderedSet, CharacteristicZero,
           LinearlyExplicitRingOver Integer)
  F : Join(FunctionSpace R, AlgebraicallyClosedField,
           TranscendentalFunctionCategory)
  L : SetCategory

  SY  ==> Symbol
  N   ==> NonNegativeInteger
  K   ==> Kernel F
  P   ==> SparseMultivariatePolynomial(R, K)
  UP  ==> SparseUnivariatePolynomial F
  RF  ==> Fraction UP
  UPUP==> SparseUnivariatePolynomial RF
  IR  ==> IntegrationResult F
  IR2 ==> IntegrationResultFunctions2(curve, F)
  ALG ==> AlgebraicIntegrate(R, F, UP, UPUP, curve)
  LDALG ==> LinearOrdinaryDifferentialOperator1 curve
  RDALG ==> PureAlgebraicLODE(F, UP, UPUP, curve)
  LOG ==> Record(coeff:F, logand:F)
  REC ==> Record(particular:U1, basis:List F)
  CND ==> Record(left:UP, right:UP)
  CHV ==> Record(int:UPUP, left:UP, right:UP, den:RF, deg:N)
  U1  ==> Union(F, "failed")
  U2  ==> Union(Record(ratpart:F, coeff:F),"failed")
  U3  ==> Union(Record(mainpart:F, limitedlogs:List LOG), "failed")
  FAIL==> error "failed - cannot handle that integrand"

  SIG ==> with

    palgint : (F, K, K) -> IR
      ++ palgint(f, x, y) returns the integral of \spad{f(x,y)dx}
      ++ where y is an algebraic function of x.

    palgextint : (F, K, K, F) -> U2
      ++ palgextint(f, x, y, g) returns functions \spad{[h, c]} such that
      ++ \spad{dh/dx = f(x,y) - c g}, where y is an algebraic function of x;
      ++ returns "failed" if no such functions exist.

    palglimint : (F, K, K, List F) -> U3
      ++ palglimint(f, x, y, [u1,...,un]) returns functions
      ++ \spad{[h,[[ci, ui]]]} such that the ui's are among \spad{[u1,...,un]}
      ++ and \spad{d(h + sum(ci log(ui)))/dx = f(x,y)} if such functions exist,
      ++ "failed" otherwise;
      ++ y is an algebraic function of x.

    palgRDE : (F, F, F, K, K, (F, F, SY) -> U1) -> U1
      ++ palgRDE(nfp, f, g, x, y, foo) returns a function \spad{z(x,y)}
      ++ such that \spad{dz/dx + n * df/dx z(x,y) = g(x,y)} if such a z exists,
      ++ "failed" otherwise;
      ++ y is an algebraic function of x;
      ++ \spad{foo(a, b, x)} is a function that solves
      ++ \spad{du/dx + n * da/dx u(x) = u(x)}
      ++ for an unknown \spad{u(x)} not involving y.
      ++ \spad{nfp} is \spad{n * df/dx}.

    if L has LinearOrdinaryDifferentialOperatorCategory F then

      palgLODE : (L, F, K, K, SY) -> REC
        ++ palgLODE(op, g, kx, y, x) returns the solution of \spad{op f = g}.
        ++ y is an algebraic function of x.

  CODE ==> add

    import IntegrationTools(R, F)
    import RationalIntegration(F, UP)
    import GenusZeroIntegration(R, F, L)
    import ChangeOfVariable(F, UP, UPUP)
    import IntegrationResultFunctions2(F, F)
    import IntegrationResultFunctions2(RF, F)
    import SparseUnivariatePolynomialFunctions2(F, RF)
    import UnivariatePolynomialCommonDenominator(UP, RF, UPUP)
    import PolynomialCategoryQuotientFunctions(IndexedExponents K,
                                                        K, R, P, F)

    quadIfCan      : (K, K) -> Union(Record(coef:F, poly:UP), "failed")
    linearInXIfCan : (K, K) -> Union(Record(xsub:F, dxsub:RF), "failed")
    prootintegrate : (F, K, K) -> IR
    prootintegrate1: (UPUP, K, K, UPUP) -> IR
    prootextint    : (F, K, K, F) -> U2
    prootlimint    : (F, K, K, List F) -> U3
    prootRDE       : (F, F, F, K, K, (F, F, SY) -> U1) -> U1
    palgRDE1       : (F, F, K, K) -> U1
    palgLODE1      : (List F, F, K, K, SY) -> REC
    palgintegrate  : (F, K, K) -> IR
    palgext        : (F, K, K, F) -> U2
    palglim        : (F, K, K, List F) -> U3
    UPUP2F1        : (UPUP, RF, RF, K, K) -> F
    UPUP2F0        : (UPUP, K, K) -> F
    RF2UPUP        : (RF, UPUP) -> UPUP
    algaddx        : (IR, F) -> IR
    chvarIfCan     : (UPUP, RF, UP, RF) -> Union(UPUP, "failed")
    changeVarIfCan : (UPUP, RF, N) -> Union(CHV, "failed")
    rationalInt    : (UPUP, N, UP) -> IntegrationResult RF
    chv            : (UPUP, N, F, F) -> RF
    chv0           : (UPUP, N, F, F) -> F
    candidates     : UP -> List CND

    dummy := new()$SY
    dumk  := kernel(dummy)@K

    UPUP2F1(p, t, cf, kx, k) == UPUP2F0(eval(p, t, cf), kx, k)

    UPUP2F0(p, kx, k)        == multivariate(p, kx, k::F)

    chv(f, n, a, b)          == univariate(chv0(f, n, a, b), dumk)

    RF2UPUP(f, modulus) ==
      bc := extendedEuclidean(map((z1:F):RF+->z1::UP::RF, denom f), modulus,
                                      1)::Record(coef1:UPUP, coef2:UPUP)
      (map((x1:F):RF+->x1::UP::RF, numer f) * bc.coef1) rem modulus

    -- returns "failed", or (xx, c) such that f(x, y)dx = f(xx, y) c dy
    -- if p(x, y) = 0 is linear in x
    linearInXIfCan(x, y) ==
      a := b := 0$UP
      p := clearDenominator lift(minPoly y, x)
      while p ^= 0 repeat
        degree(q := numer leadingCoefficient p) > 1 => return "failed"
        a := a + monomial(coefficient(q, 1), d := degree p)
        b := b - monomial(coefficient(q, 0), d)
        p := reductum p
      xx:RF := b / a
      [xx(dumk::F), differentiate(xx, differentiate)]

    -- return Int(f(x,y)dx) where y is an n^th root of a rational function in x
    prootintegrate(f, x, y) ==
      modulus := lift(p := minPoly y, x)
      rf      := reductum(ff := univariate(f, x, y, p))
      ((r := retractIfCan(rf)@Union(RF,"failed")) case RF) and rf ^= 0 =>
            -- in this case, ff := lc(ff) y^i + r so we integrate both terms
            -- separately to gain time
            map(f1+->f1(x::F), integrate(r::RF)) +
                 prootintegrate1(leadingMonomial ff, x, y, modulus)
      prootintegrate1(ff, x, y, modulus)

    prootintegrate1(ff, x, y, modulus) ==
      chv:CHV
      r := radPoly(modulus)::Record(radicand:RF, deg:N)
      (uu := changeVarIfCan(ff, r.radicand, r.deg)) case CHV =>
        chv := uu::CHV
        newalg := nthRoot((chv.left)(dumk::F), chv.deg)
        kz := retract(numer newalg)@K
        newf := multivariate(chv.int, ku := dumk, newalg)
        vu := (chv.right)(x::F)
        vz := (chv.den)(x::F) * (y::F) * denom(newalg)::F
        map(x1+->eval(x1, [ku, kz], [vu, vz]), palgint(newf, ku, kz))
      cv     := chvar(ff, modulus)
      r      := radPoly(cv.poly)::Record(radicand:RF, deg:N)
      qprime := differentiate(q := retract(r.radicand)@UP)::RF
      not zero? qprime and
       ((u := chvarIfCan(cv.func, 1, q, inv qprime)) case UPUP) =>
         m := monomial(1, r.deg)$UPUP - q::RF::UPUP
         map(x1+->UPUP2F1(RF2UPUP(x1, m), cv.c1, cv.c2, x, y),
            rationalInt(u::UPUP, r.deg, monomial(1, 1)))
      curve  := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg)
      algaddx(map(x1+->UPUP2F1(lift x1, cv.c1, cv.c2, x, y),
        palgintegrate(reduce(cv.func), differentiate$UP)$ALG)$IR2, x::F)

    -- Do the rationalizing change of variable
    -- Int(f(x, y) dx) --> Int(n u^(n-1) f((u^n - b)/a, u) / a du) where
    -- u^n = y^n = g(x) = a x + b
    -- returns the integral as an integral of a rational function in u
    rationalInt(f, n, g) ==
      not ((degree g) = 1) => error "rationalInt: radicand must be linear"
      a := leadingCoefficient g
      integrate(n * monomial(inv a, (n-1)::N)$UP
                  * chv(f, n, a, leadingCoefficient reductum g))

    -- Do the rationalizing change of variable 
    -- f(x,y) --> f((u^n - b)/a, u) where
    -- u = y = (a x + b)^(1/n).
    -- Returns f((u^n - b)/a,u) as an element of F
    chv0(f, n, a, b) ==
      d := dumk::F
      (f (d::UP::RF)) ((d ** n - b) / a)

    -- candidates(p) returns a list of pairs [g, u] such that p(x) = g(u(x)),
    -- those u's are candidates for change of variables
    -- currently uses a dumb heuristic where the candidates u's are p itself
    -- and all the powers x^2, x^3, ..., x^{deg(p)},
    -- will use polynomial decomposition in smarter days   MB 8/93
    candidates p ==
      l:List(CND) := empty()
      ground? p => l
      for i in 2..degree p repeat
        if (u := composite(p, xi := monomial(1, i))) case UP then
          l := concat([u::UP, xi], l)
      concat([monomial(1, 1), p], l)

    -- checks whether Int(p(x, y) dx) can be rewritten as
    -- Int(r(u, z) du) where u is some polynomial of x,
    -- z = d y for some polynomial d, and z^m = g(u)
    -- returns either [r(u, z), g, u, d, m] or "failed"
    -- we have y^n = radi
    changeVarIfCan(p, radi, n) ==
      rec := rootPoly(radi, n)
      for cnd in candidates(rec.radicand) repeat
        (u := chvarIfCan(p, rec.coef, cnd.right,
              inv(differentiate(cnd.right)::RF))) case UPUP =>
                 return [u::UPUP, cnd.left, cnd.right, rec.coef, rec.exponent]
      "failed"

    -- checks whether Int(p(x, y) dx) can be rewritten as
    -- Int(r(u, z) du) where u is some polynomial of x and z = d y
    -- we have y^n = a(x)/d(x)
    -- returns either "failed" or r(u, z)
    chvarIfCan(p, d, u, u1) ==
      ans:UPUP := 0
      while p ^= 0 repeat
        (v := composite(u1 * leadingCoefficient(p) / d ** degree(p), u))
          case "failed" => return "failed"
        ans := ans + monomial(v::RF, degree p)
        p   := reductum p
      ans

    algaddx(i, xx) ==
      elem? i => i
      mkAnswer(ratpart i, logpart i,
                [[- ne.integrand / (xx**2), xx] for ne in notelem i])

    prootRDE(nfp, f, g, x, k, rde) ==
      modulus := lift(p := minPoly k, x)
      r       := radPoly(modulus)::Record(radicand:RF, deg:N)
      rec     := rootPoly(r.radicand, r.deg)
      dqdx    := inv(differentiate(q := rec.radicand)::RF)
      ((uf := chvarIfCan(ff := univariate(f,x,k,p),rec.coef,q,1)) case UPUP) _
       and _
        ((ug:=chvarIfCan(gg:=univariate(g,x,k,p),rec.coef,q,dqdx))_
             case UPUP) =>
          (u := rde(chv0(uf::UPUP, rec.exponent, 1, 0), rec.exponent *
                    (dumk::F) ** (rec.exponent * (rec.exponent - 1))
                      * chv0(ug::UPUP, rec.exponent, 1, 0),
                       symbolIfCan(dumk)::SY)) case "failed" => "failed"
          eval(u::F, dumk, k::F)
      ((rec.coef) = 1) =>
        curve  := RadicalFunctionField(F, UP, UPUP, q::RF, rec.exponent)
        rc := algDsolve(D()$LDALG + reduce(univariate(nfp, x, k, p))::LDALG,
                         reduce univariate(g, x, k, p))$RDALG
        rc.particular case "failed" => "failed"
        UPUP2F0(lift((rc.particular)::curve), x, k)
      palgRDE1(nfp, g, x, k)

    prootlimint(f, x, k, lu) ==
      modulus := lift(p := minPoly k, x)
      r       := radPoly(modulus)::Record(radicand:RF, deg:N)
      rec     := rootPoly(r.radicand, r.deg)
      dqdx    := inv(differentiate(q := rec.radicand)::RF)
      (uf:=chvarIfCan(ff := univariate(f,x,k,p),rec.coef,q,dqdx)) case UPUP =>
        l := empty()$List(RF)
        n := rec.exponent * monomial(1, (rec.exponent - 1)::N)$UP
        for u in lu repeat
         if ((v:=chvarIfCan(uu:=univariate(u,x,k,p),rec.coef,q,dqdx))case UPUP)
            then l := concat(n * chv(v::UPUP,rec.exponent, 1, 0), l) else FAIL
        m := monomial(1, rec.exponent)$UPUP - q::RF::UPUP
        map(x1+->UPUP2F0(RF2UPUP(x1,m), x, k),
            limitedint(n * chv(uf::UPUP, rec.exponent, 1, 0), reverse_! l))
      cv     := chvar(ff, modulus)
      r      := radPoly(cv.poly)::Record(radicand:RF, deg:N)
      dqdx   := inv(differentiate(q := retract(r.radicand)@UP)::RF)
      curve  := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg)
      (ui := palginfieldint(reduce(cv.func), differentiate$UP)$ALG)
        case "failed" => FAIL
      [UPUP2F1(lift(ui::curve), cv.c1, cv.c2, x, k), empty()]

    prootextint(f, x, k, g) ==
      modulus := lift(p := minPoly k, x)
      r       := radPoly(modulus)::Record(radicand:RF, deg:N)
      rec     := rootPoly(r.radicand, r.deg)
      dqdx    := inv(differentiate(q := rec.radicand)::RF)
      ((uf:=chvarIfCan(ff:=univariate(f,x,k,p),rec.coef,q,dqdx)) case UPUP) and
       ((ug:=chvarIfCan(gg:=univariate(g,x,k,p),rec.coef,q,dqdx)) case UPUP) =>
          m := monomial(1, rec.exponent)$UPUP - q::RF::UPUP
          n := rec.exponent * monomial(1, (rec.exponent - 1)::N)$UP
          map(x1+->UPUP2F0(RF2UPUP(x1,m), x, k),
              extendedint(n * chv(uf::UPUP, rec.exponent, 1, 0),
                          n * chv(ug::UPUP, rec.exponent, 1, 0)))
      cv     := chvar(ff, modulus)
      r      := radPoly(cv.poly)::Record(radicand:RF, deg:N)
      dqdx   := inv(differentiate(q := retract(r.radicand)@UP)::RF)
      curve  := RadicalFunctionField(F, UP, UPUP, q::RF, r.deg)
      (u := palginfieldint(reduce(cv.func), differentiate$UP)$ALG)
        case "failed" => FAIL
      [UPUP2F1(lift(u::curve), cv.c1, cv.c2, x, k), 0]

    palgRDE1(nfp, g, x, y) ==
      palgLODE1([nfp, 1], g, x, y, symbolIfCan(x)::SY).particular

    palgLODE1(eq, g, kx, y, x) ==
      modulus:= lift(p := minPoly y, kx)
      curve  := AlgebraicFunctionField(F, UP, UPUP, modulus)
      neq:LDALG := 0
      for f in eq for i in 0.. repeat
          neq := neq + monomial(reduce univariate(f, kx, y, p), i)
      empty? remove_!(y, remove_!(kx, varselect(kernels g, x))) =>
        rec := algDsolve(neq, reduce univariate(g, kx, y, p))$RDALG
        bas:List(F) := [UPUP2F0(lift h, kx, y) for h in rec.basis]
        rec.particular case "failed" => ["failed", bas]
        [UPUP2F0(lift((rec.particular)::curve), kx, y), bas]
      rec := algDsolve(neq, 0)
      ["failed", [UPUP2F0(lift h, kx, y) for h in rec.basis]]

    palgintegrate(f, x, k) ==
      modulus:= lift(p := minPoly k, x)
      cv     := chvar(univariate(f, x, k, p), modulus)
      curve  := AlgebraicFunctionField(F, UP, UPUP, cv.poly)
      knownInfBasis(cv.deg)
      algaddx(map(x1+->UPUP2F1(lift x1, cv.c1, cv.c2, x, k),
        palgintegrate(reduce(cv.func), differentiate$UP)$ALG)$IR2, x::F)

    palglim(f, x, k, lu) ==
      modulus:= lift(p := minPoly k, x)
      cv     := chvar(univariate(f, x, k, p), modulus)
      curve  := AlgebraicFunctionField(F, UP, UPUP, cv.poly)
      knownInfBasis(cv.deg)
      (u := palginfieldint(reduce(cv.func), differentiate$UP)$ALG)
        case "failed" => FAIL
      [UPUP2F1(lift(u::curve), cv.c1, cv.c2, x, k), empty()]

    palgext(f, x, k, g) ==
      modulus:= lift(p := minPoly k, x)
      cv     := chvar(univariate(f, x, k, p), modulus)
      curve  := AlgebraicFunctionField(F, UP, UPUP, cv.poly)
      knownInfBasis(cv.deg)
      (u := palginfieldint(reduce(cv.func), differentiate$UP)$ALG)
        case "failed" => FAIL
      [UPUP2F1(lift(u::curve), cv.c1, cv.c2, x, k), 0]

    palgint(f, x, y) ==
      (v := linearInXIfCan(x, y)) case "failed" =>
        (u := quadIfCan(x, y)) case "failed" =>
          is?(y, "nthRoot"::SY) => prootintegrate(f, x, y)
          is?(y,  "rootOf"::SY) => palgintegrate(f, x, y)
          FAIL
        palgint0(f, x, y, u.coef, u.poly)
      palgint0(f, x, y, dumk, v.xsub, v.dxsub)

    palgextint(f, x, y, g) ==
      (v := linearInXIfCan(x, y)) case "failed" =>
        (u := quadIfCan(x, y)) case "failed" =>
          is?(y, "nthRoot"::SY) => prootextint(f, x, y, g)
          is?(y,  "rootOf"::SY) => palgext(f, x, y, g)
          FAIL
        palgextint0(f, x, y, g, u.coef, u.poly)
      palgextint0(f, x, y, g, dumk, v.xsub, v.dxsub)

    palglimint(f, x, y, lu) ==
      (v := linearInXIfCan(x, y)) case "failed" =>
        (u := quadIfCan(x, y)) case "failed" =>
          is?(y, "nthRoot"::SY) => prootlimint(f, x, y, lu)
          is?(y,  "rootOf"::SY) => palglim(f, x, y, lu)
          FAIL
        palglimint0(f, x, y, lu, u.coef, u.poly)
      palglimint0(f, x, y, lu, dumk, v.xsub, v.dxsub)

    palgRDE(nfp, f, g, x, y, rde) ==
      (v := linearInXIfCan(x, y)) case "failed" =>
        (u := quadIfCan(x, y)) case "failed" =>
          is?(y, "nthRoot"::SY) => prootRDE(nfp, f, g, x, y, rde)
          palgRDE1(nfp, g, x, y)
        palgRDE0(f, g, x, y, rde, u.coef, u.poly)
      palgRDE0(f, g, x, y, rde, dumk, v.xsub, v.dxsub)

    -- returns "failed", or (d, P) such that (dy)**2 = P(x)
    -- and degree(P) = 2
    quadIfCan(x, y) ==
      (degree(p := minPoly y) = 2) and zero?(coefficient(p, 1)) =>
        d := denom(ff :=
                 univariate(- coefficient(p, 0) / coefficient(p, 2), x))
        degree(radi := d * numer ff) = 2 => [d(x::F), radi]
        "failed"
      "failed"

    if L has LinearOrdinaryDifferentialOperatorCategory F then

      palgLODE(eq, g, kx, y, x) ==
        (v := linearInXIfCan(kx, y)) case "failed" =>
          (u := quadIfCan(kx, y)) case "failed" =>
            palgLODE1([coefficient(eq, i) for i in 0..degree eq], g, kx, y, x)
          palgLODE0(eq, g, kx, y, u.coef, u.poly)
        palgLODE0(eq, g, kx, y, dumk, v.xsub, v.dxsub)