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

)abbrev package RDETR TranscendentalRischDE
++ Author: Manuel Bronstein
++ Date Created: Jan 1988
++ Date Last Updated: 2 November 1995
++ References:
++ Bron90 The Transcendental Risch Differential Equation
++ Description:
++ Risch differential equation, transcendental case.

TranscendentalRischDE(F, UP) : SIG == CODE where
  F  : Join(Field, CharacteristicZero, RetractableTo Integer)
  UP : UnivariatePolynomialCategory F

  N   ==> NonNegativeInteger
  Z   ==> Integer
  RF  ==> Fraction UP
  REC ==> Record(a:UP, b:UP, c:UP, t:UP)
  SPE ==> Record(b:UP, c:UP, m:Z, alpha:UP, beta:UP)
  PSOL==> Record(ans:UP, nosol:Boolean)
  ANS ==> Union(ans:PSOL, eq:SPE)
  PSQ ==> Record(ans:RF, nosol:Boolean)

  SIG ==> with

    monomRDE : (RF,RF,UP->UP) -> Union(Record(a:UP,b:RF,c:RF,t:UP), "failed")
      ++ monomRDE(f,g,D) returns \spad{[A, B, C, T]} such that
      ++ \spad{y' + f y = g} has a solution if and only if \spad{y = Q / T},
      ++ where Q satisfies \spad{A Q' + B Q = C} and has no normal pole.
      ++ A and T are polynomials and B and C have no normal poles.
      ++ D is the derivation to use.

    baseRDE : (RF, RF) -> PSQ
      ++ baseRDE(f, g) returns a \spad{[y, b]} such that \spad{y' + fy = g}
      ++ if \spad{b = true}, y is a partial solution otherwise (no solution
      ++ in that case).
      ++ D is the derivation to use.

    polyRDE : (UP, UP, UP, Z, UP -> UP) -> ANS
      ++ polyRDE(a, B, C, n, D) returns either:
      ++ 1. \spad{[Q, b]} such that \spad{degree(Q) <= n} and
      ++    \spad{a Q'+ B Q = C} if \spad{b = true}, Q is a partial solution
      ++    otherwise.
      ++ 2. \spad{[B1, C1, m, \alpha, \beta]} such that any polynomial solution
      ++    of degree at most n of \spad{A Q' + BQ = C} must be of the form
      ++    \spad{Q = \alpha H + \beta} where \spad{degree(H) <= m} and
      ++    H satisfies \spad{H' + B1 H = C1}.
      ++ D is the derivation to use.

  CODE ==> add

    import MonomialExtensionTools(F, UP)

    getBound     : (UP, UP, Z) -> Z
    SPDEnocancel1: (UP, UP, Z, UP -> UP) -> PSOL
    SPDEnocancel2: (UP, UP, Z, Z, F, UP -> UP) -> ANS
    SPDE         : (UP, UP, UP, Z, UP -> UP) -> Union(SPE, "failed")

-- cancellation at infinity is possible, A is assumed nonzero
-- needs tagged union because of branch choice problem
-- always returns a PSOL in the base case (never a SPE)
    polyRDE(aa, bb, cc, d, derivation) ==
      n:Z
      (u := SPDE(aa, bb, cc, d, derivation)) case "failed" => [[0, true]]
      zero?(u.c) => [[u.beta, false]]
      baseCase? := ((dt := derivation monomial(1, 1)) = 1)
      n := degree(dt)::Z - 1
      b0? := zero?(u.b)
      (~b0?) and (baseCase? or degree(u.b) > max(0, n)) =>
          answ := SPDEnocancel1(u.b, u.c, u.m, derivation)
          [[u.alpha * answ.ans + u.beta, answ.nosol]]
      (n > 0) and (b0? or degree(u.b) < n) =>
          uansw:= SPDEnocancel2(u.b,u.c,u.m,n,leadingCoefficient dt,derivation)
          uansw case ans=> [[u.alpha * uansw.ans.ans + u.beta,uansw.ans.nosol]]
          [[uansw.eq.b, uansw.eq.c, uansw.eq.m,
            u.alpha * uansw.eq.alpha, u.alpha * uansw.eq.beta + u.beta]]
      b0? and baseCase? =>
          degree(u.c) >= u.m => [[0, true]]
          [[u.alpha * integrate(u.c) + u.beta, false]]
      [u::SPE]

-- cancellation at infinity is possible, A is assumed nonzero
-- if u.b = 0 then u.a = 1 already, but no degree check is done
-- returns "failed" if a p' + b p = c has no soln of degree at most d,
-- otherwise [B, C, m, \alpha, \beta] such that any soln p of degree at
-- most d of  a p' + b p = c  must be of the form p = \alpha h + \beta,
-- where h' + B h = C and h has degree at most m
    SPDE(aa, bb, cc, d, derivation) ==
      zero? cc => [0, 0, 0, 0, 0]
      d < 0 => "failed"
      (u := cc exquo (g := gcd(aa, bb))) case "failed" => "failed"
      aa := (aa exquo g)::UP
      bb := (bb exquo g)::UP
      cc := u::UP
      (ra := retractIfCan(aa)@Union(F, "failed")) case F =>
        a1 := inv(ra::F)
        [a1 * bb, a1 * cc, d, 1, 0]
      bc := extendedEuclidean(bb, aa, cc)::Record(coef1:UP, coef2:UP)
      qr := divide(bc.coef1, aa)
      r  := qr.remainder         -- z = bc.coef2 + b * qr.quotient
      (v  := SPDE(aa, bb + derivation aa,
                  bc.coef2 + bb * qr.quotient - derivation r,
                   d - degree(aa)::Z, derivation)) case "failed" => "failed"
      [v.b, v.c, v.m, aa * v.alpha, aa * v.beta + r]

-- solves q' + b q = c  with deg(q) <= d
-- case (B <> 0) and (D = d/dt or degree(B) > max(0, degree(Dt) - 1))
-- this implies no cancellation at infinity, BQ term dominates
-- returns [Q, flag] such that Q is a solution if flag is false,
-- a partial solution otherwise.
    SPDEnocancel1(bb, cc, d, derivation) ==
      q:UP := 0
      db := (degree bb)::Z
      lb := leadingCoefficient bb
      while cc ^= 0 repeat
        d < 0 or (n := (degree cc)::Z - db) < 0 or n > d => return [q, true]
        r := monomial((leadingCoefficient cc) / lb, n::N)
        cc := cc - bb * r - derivation r
        d := n - 1
        q := q + r
      [q, false]

-- case (t is a nonlinear monomial) and (B = 0 or degree(B) < degree(Dt) - 1)
-- this implies no cancellation at infinity, DQ term dominates or degree(Q) = 0
-- dtm1 = degree(Dt) - 1
    SPDEnocancel2(bb, cc, d, dtm1, lt, derivation) ==
      q:UP := 0
      while cc ^= 0 repeat
        d < 0 or (n := (degree cc)::Z - dtm1) < 0 or n > d => return [[q,true]]
        if n > 0 then
          r  := monomial((leadingCoefficient cc) / (n * lt), n::N)
          cc := cc - bb * r - derivation r
          d  := n - 1
          q  := q + r
        else        -- n = 0 so solution must have degree 0
          db:N := (zero? bb => 0; degree bb);
          db ^= degree(cc) => return [[q, true]]
          zero? db => return [[bb, cc, 0, 1, q]]
          r  := leadingCoefficient(cc) / leadingCoefficient(bb)
          cc := cc - r * bb - derivation(r::UP)
          d  := - 1
          q := q + r::UP
      [[q, false]]

    monomRDE(f, g, derivation) ==
      gg := gcd(d := normalDenom(f,derivation), e := normalDenom(g,derivation))
      tt := (gcd(e, differentiate e) exquo gcd(gg,differentiate gg))::UP
      (u := ((tt * (aa := d * tt)) exquo e)) case "failed" => "failed"
      [aa, aa * f - (d * derivation tt)::RF, u::UP * e * g, tt]

-- solve y' + f y = g for y in RF
-- assumes that f is weakly normalized (no finite cancellation)
-- base case: F' = 0
    baseRDE(f, g) ==
      (u := monomRDE(f, g, differentiate)) case "failed" => [0, true]
      n := getBound(u.a,bb := retract(u.b)@UP,degree(cc := retract(u.c)@UP)::Z)
      v := polyRDE(u.a, bb, cc, n, differentiate).ans
      [v.ans / u.t, v.nosol]

-- return an a bound on the degree of a solution of A P'+ B P = C,A ^= 0
-- cancellation at infinity is possible
-- base case: F' = 0
    getBound(a, b, dc) ==
      da := (degree a)::Z
      zero? b => max(0, dc - da + 1)
      db := (degree b)::Z
      da > (db + 1) => max(0, dc - da + 1)
      da < (db + 1) => dc - db
      (n := retractIfCan(- leadingCoefficient(b) / leadingCoefficient(a)
                      )@Union(Z, "failed")) case Z => max(n::Z, dc - db)
      dc - db