share/integration/hermite_reduce.mac

   1 /* nicked from mailing list; here is the email message + headers:
   2  *
   3 From: toy.raymond at gmail.com (Raymond Toy)
   4 Date: Mon, 03 May 2010 09:41:08 -0400
   5 Subject: [Maxima] Rothstein-Trager algorithm.
   6
   7 Here is a very rudimentary implementation of the Rothstein-Trager
   8 algorithm.  Seems to work for the few tests that I've run.
   9
  10 f:(x^4-3*x^2+6)/(x^6-5*x^4+5*x^2+4)$
  11 int_rational_function(f,x);
  12 -> lsum(t*log(x^3+t*x^2-3*x-2*t)/2,t,rootsof(t^2+1,t))
  13
  14 One refinement that needs to be done is to try to solve for the roots
  15 and substitute them in if possible, instead of leaving the result as a
  16 sum over the roots.
  17
  18 We could also implement the Lazard-Rioboo-Trager algorithm, but maxima
  19 doesn't seem to provide the subresultant prs sequence.  That would have
  20 to be implemented.
  21
  22 Ray
  23  *
  24  */
  25
  26 /*
  27  * Simple implementation of Rothstein-Trager algorithm for integrating
  28  * rational functions.
  29  *
  30  * $Id: hermite_reduce.mac,v 1.4 2010/05/03 13:28:46 toy Exp $
  31  */
  32
  33 /* Extended Euclidean algorithm, diophantine version
  34  *
  35  * Find s and t such that s*a + t*b = c and either s = 0 or deg(s) < deg(b)
  36  *
  37  * var is the variable of the polynomials.
  38  */
  39 extended_euclidean(a, b, c, var) :=
  40   block([s, t, g, q, r],
  41     [s, t, g] : gcdex(a, b, var),
  42     [q, r] : divide(c, g),
  43     s : q*s,
  44     t : q*t,
  45     if is(s # 0) and is(hipow(s, var) > hipow(b, var)) then
  46       ([q, r] : divide(s, b), s : r, t : t + q*a),
  47     [s, t]);
  48
  49 /* Hermite Reduction, Mack's linear version
  50  *
  51  * Find g and h such that A/D = diff(g,x) + h and return g and h.
  52  * h has a square-free denominator.
  53  */
  54 hermite_reduce(a, d, var) :=
  55   block([g : 0, d_minus: gcd(d, diff(d, var)), d_star, deg],
  56     d_star : quotient(d, d_minus),
  57     /*
  58     print("d- = ", d_minus),
  59     print("d* = ", d_star),
  60     */
  61     for i : 1 while hipow(d_minus, var) > 0 do
  62       block([d_minus2 : gcd(d_minus, diff(d_minus, var)), d_minus_s, b, c],
  63         d_minus_s : quotient(d_minus, d_minus2),
  64     /*
  65         print("d-2 = ", d_minus2),
  66         print("d-* = ", d_minus_s),
  67     */
  68         [b,c] : extended_euclidean(quotient(-d_star * diff(d_minus,var),
  69 d_minus), d_minus_s, a, var),
  70         /*print("[b, c] = ", [b, c]),*/
  71         a : c - quotient(diff(b, var) * d_star, d_minus_s),
  72         /*print("a = ", a),*/
  73         g : g + b/d_minus,
  74         /*print("g = ", g),*/
  75         d_minus : d_minus2),
  76     [g, a/d_star]);
  77
  78 /* Make p a monic polynomial.
  79  *
  80  * Let p(x) = a[0]*x^n + a[1]*x^(n-1) ...
  81  *
  82  * Use the substitution x = y/a[0]^(1/n) to get
  83  *
  84  * p(y) = y^n + a[1]/a[0]^((n-1)/n)*y^(n-1) ...
  85  *
  86  * Return the new polynomial and a[0]^(1/n)
  87  */
  88 monic(p, var) :=
  89   block([e : hipow(p, var), c, temp],
  90     c : coeff(p, var, e),
  91     [subst([var = var/(c^(1/e))], p), c^(1/e)]);
  92
  93 /* Algebraic GCD */
  94 algebraic_gcd(a,b) :=
  95   block([algebraic:true, gcd : subres],
  96     gcd(a,b));
  97
  98 /* Rothstein-Trager algorithm
  99  *
 100  * The numerator is a, the denominator is d, and the variable of
 101  * integration is var.
 102  */
 103 int_rational_log_part(a, d, var) :=
 104   block([resultant : subres, c, t, r, rr, s : 0],
 105     r : resultant(d, a - t*diff(d,var), var),
 106     r : sqfr(r),
 107     /* Can this be done better?  Say we have -c*R1(x)*R2(x)....  We
 108      * really just want the term product terms, without the leading
 109      * coefficient.  Just look for a unary minus and negate if it is.
 110      */
 111     if op(r) = "-" and length(args(r)) = 1 then
 112       r : -r,
 113     for i : 1 thru length(r) do
 114       block([f : part(r, i)],
 115         /* Skip any expressions that don't involve t */
 116         if not freeof(t, f) then
 117       block([g, p, e],
 118         /* Look for a expression like (t^2-16)^3.  Handle (t^2-16)
 119 carefully. */
 120         if op(f) = "^" then
 121           [p, e] : args(f)
 122             else
 123               [p, e] : [f, 1],
 124         [rr, c] : monic(p, t),
 125         tellrat(rr),
 126         g : algebraic_gcd(d, a - t*expand(diff(d, var)/c)),
 127         untellrat(t),
 128         g : log(g),
 129             s : s + apply('lsum, [t*g/c, t, rootsof(rr, t)]))),
 130       s);
 131
 132
 133 /* Integrate a rational function f */
 134 int_rational_function(f, var) :=
 135   block([g,h,q,r],
 136     [g, h] : hermite_reduce(num(f), denom(f), var),
 137     [q, r] : divide(num(h), denom(h)),
 138     if r = 0 then
 139       g + integrate(q, var)
 140     else
 141       g + integrate(q, var) + int_rational_log_part(r, denom(h), var));