1/*
2  Copyright 2009 by Barton Willis
3  Maxima code for integration of some algebraic functions.
4
5  This is free software; you can redistribute it and/or
6  modify it under the terms of the GNU General Public License,
7  http://www.gnu.org/copyleft/gpl.html.
8
9 This software has NO WARRANTY, not even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
11
12*/
13
14
15load("basic");
16if get('abs_integrate,'version) = false then load("abs_integrate.mac")$
17
18/* The function radcan isn't idempotent for some expressions that you might think
19   that it should be; for example
20
21     ((-x-1)^(-bb-aa-2)*(x+2)^aa)/(((x+2)/(x+1))^aa*(1/(x+1)-(x+2)/(x +1)^2)*(1-(x+2)/(x+1))^bb)
22
23   The function infapply works around this by doing fixed point iteration on a function f. */
24
25infapply(e,f) := block([ee],
26  while e # (ee : apply(f,[e])) do e : ee, e);
27
28
29/* Return a list of the coefficients of a polynomial p[x]. This function doesn't check
30   that p is a polynomial. */
31
32all_poly_coefs(p,x) := (p : ratexpand(p), makelist(ratcoef(p,x,k),k,0, hipow(p,x)));
33
34collect_poly_factors(e,x) := block([l],
35  e : ratsimp(diff(e,x) / e),
36  e : factor(ratdenom(e)),
37  l : if safe_op(e) = "*" then args(e) else [e],
38  l : sublist(l, lambda([s], not freeof(x,s))),
39  if every(lambda([s], polynomialp(s,[x], lambda([u], freeof(x,u)), 'integerp)), l) then
40  setify(map(lambda([s], s / ratcoef(s,x,hipow(s,x))), l)) else false);
41
42generate_hyper_mu(p, r, a, b) := block([l : [], f, n, i, pt, rt,k],
43  n : 3^length(p) - 1,
44  while n > 0 do (
45    f : 0,
46    pt : p, /* pt = tail of p */
47    rt : r,
48    k : n,
49    i : 0, /* i counts the number of nonzero terms in f */
50    while k > 0 do (
51      k : divide(k,3),
52      if second(k) = 1 then (
53        i : i + 1,
54        f : f + first(pt) * (first(rt) + 1)/(a + 1))
55      else if second(k) = 2 then (
56        i : i + 1,
57        f : f + first(pt) * (first(rt) + 1)/(a + b + 1)),
58      k : first(k),
59      rt : rest(rt),
60      pt : rest(pt)),
61    n : n - 1,
62    l : push([i,f],l)),
63  l : sort(l, lambda([a,b], first(a) < first(b))),
64  l : map('second,l));
65
66/* Express sigma as a linear combination of the members of kern. */
67
68xresidue(sigma, kern, x) := block([zip, %g, i, n, k],
69  %g : new_variable('general),
70  n : length(kern),
71  zip : partfrac(sigma - kern . makelist(concat(%g,i),i,1,n),x),
72  zip : if safe_op(zip) = "+" then args(zip) else [zip],
73  map('rhs, linsolve(map('ratnumer,zip), makelist(concat(%g,k),k,1,n))));
74
75/* Return true if the function x |--> e is piecewise constant. We'll assume that
76if g is piecewise constant, so is the composition f(g). This isn't true for all
77functions f, (say f is 1 on the rational numbers and 0 otherwise). */
78
79piecewise_constant_p(e,x,[z]) := (
80  if z = [] then (
81    e : infapply(e,lambda([s], convert_to_signum(radcan(s))))),
82  if mapatom(e) then freeof(x,e)
83  else (
84    numberp(e) or
85    member(safe_op(e), ['floor, 'signum]) or
86    every(lambda([s], piecewise_constant_p(s,x, false)), args(e))));
87
88hyper_int(e,x) := block([kern, sigma, a, b, r, dkern, mu, eq, aa, l,
89                         sol, f, cnst, ie : false, %fo, %so, algexact : true],
90  kern : collect_poly_factors(e,x),
91  if kern # false then (
92    kern : listify(kern),
93    sigma : partfrac(ratsimp(diff(e,x) / e),x),
94    a : new_variable('general),
95    b : new_variable('general),
96
97    dkern : map(lambda([s], diff(s,x) / s), kern),
98    r : xresidue(sigma, dkern, x),
99    l : if r = [ ] then [ ] else generate_hyper_mu(dkern, r, a, b),
100    for mu in l while ie = false do (
101
102      eq : mu^2*(b*diff(sigma,x,1)+sigma^2+b*diff(mu,x,1)+2*a*diff(mu,x,1)+2
103        *diff(mu,x,1))-mu*(b*diff(mu,x,1)*sigma+2*diff(mu,x,1)*sigma+b
104        *diff(mu,x,2))-(b+2*a+2)*mu^3*sigma+(2*b+1)*(diff(mu,x,1))^2+(a+1) *(b+a+1)*mu^4,
105
106        eq : partfrac(eq,x),
107        eq : if safe_op(eq) = "+" then args(eq) else [eq],
108        eq : xreduce('append, map(lambda([s], all_poly_coefs(ratnumer(s),x)),eq)),
109        eq : listify(setify(eq)),
110        sol : algsys(eq, [a,b]),
111        for sk in sol while ie = false do (
112          f : radcan(exp(logcontract(integrate(mu,x)))),
113          f : block([?errorsw : true], errcatch(subst(sk, f))),
114          if f #  [] then (
115            f : first(f),
116            %fo : new_variable('general),
117            f : %fo * f,
118            cnst : subst(sk, diff(f,x) * f^a * (1 - f)^b / e),
119            cnst : infapply(cnst,lambda([s], radcan(rootscontract(s)))),
120            if cnst # 0 then (
121              cnst : ratnumer(ratsimp(diff(cnst,x) / cnst)),
122              cnst : solve(cnst, %fo),
123              aa : subst(sk,a))
124            else cnst : [ ],
125            if cnst # [] and not(integerp(aa) and aa < 0) then (
126              f : subst(cnst,f),
127              %so : ratsimp(subst(sk, diff(f,x) * f^a *  (1 - f)^b / e)),
128              if %so # 0 and piecewise_constant_p(%so,x) = true then (
129                ie : subst(sk, (1/%so) * f^(a+1) * hypergeometric([a+1,-b],[2+a], f)/(a+1)))))))),
130    ie);
131
132generate_elliptic_candidates(p, r, %k) := block([l : [], f, n, i, pt, rt, k],
133  n : 5^length(p) - 1,
134  while n > 0 do (
135    f : 0,
136    pt : p, /* pt = tail of p */
137    rt : r,
138    k : n,
139    i : 0, /* i counts the number of terms in f that involve %k */
140    while k > 0 do (
141      k : divide(k,5),
142      if second(k) = 1 then (
143        f : f + first(pt) * (first(rt) + 1))
144      else if second(k) = 2 then (
145        f : f - first(pt) * (first(rt) + 1))
146      else if second(k) = 3 then (
147        i : i + 1,
148        f : f + first(pt) * %i * (%k^2 - 1)*(first(rt) + 1) /(2 * %k))
149      else if second(k) = 4 then (
150        i : i + 1,
151        f : f - first(pt) * %i * (%k^2 - 1)*(first(rt) + 1) /(2 * %k)),
152      k : first(k),
153      rt : rest(rt),
154      pt : rest(pt)),
155    n : n - 1,
156    l : push([i,f],l)),
157  l : sort(l, lambda([a,b], first(a) < first(b))),
158  l : map('second,l));
159
160inverse_jacobi_int(e,x) := block([kern, sigma, %k, r, dkern, mu, eq, sol, f, cnst, ie : false,l,
161                              %fo, %so, algexact : true],
162  kern : collect_poly_factors(e,x),
163  if kern # false then (
164    kern : listify(kern),
165    sigma : partfrac(ratsimp(diff(e,x) / e),x),
166    %k : new_variable('general),
167    dkern : map(lambda([s], diff(s,x) / s), kern),
168    r : xresidue(sigma, dkern, x),
169    l : if r = [ ] then []  else generate_elliptic_candidates(dkern, r, %k),
170    for mu in l while ie = false do (
171
172      eq : -2*%k^2*(mu^2*(diff(sigma,x,1))^2-4*mu^2*sigma^2*(diff(sigma,x,1))+
173        6*mu*(diff(mu,x,1))*sigma*(diff(sigma,x,1))-2*mu*(diff(mu,x,2))*(diff(sigma,x,1))+
174        8*mu^4*(diff(sigma,x,1))+4*mu^2*sigma^4-12*mu*(diff(mu,x,1))*sigma^3+
175        4*mu*(diff(mu,x,2))*sigma^2+9* (diff(mu,x,1))^2*sigma^2-12*mu^4*sigma^2-
176        6*(diff(mu,x,1))*(diff(mu,x,2))*sigma+16*mu^3*(diff(mu,x,1))*sigma+
177        (diff(mu,x,2))^2-8*mu^3*(diff(mu,x,2))+4*mu^2*(diff(mu,x,1))^2+8*mu^6)+
178      %k^4* (mu*(diff(sigma,x,1))-2*mu*sigma^2+3*(diff(mu,x,1))*sigma-2*mu^2*sigma-
179        diff(mu,x,2)+2*mu*(diff(mu,x,1)))*(mu*(diff(sigma,x,1))-2*mu*sigma^2+
180        3*(diff(mu,x,1))*sigma+2*mu^2*sigma-diff(mu,x,2)-2*mu*(diff(mu,x,1)))+
181      (mu*(diff(sigma,x,1))-2*mu*sigma^2+3*(diff(mu,x,1))*sigma-2*mu^2*sigma-
182        diff(mu,x,2)+2*mu*(diff(mu,x,1)))*(mu*(diff(sigma,x,1))-2*mu*sigma^2+
183        3*(diff(mu,x,1))*sigma+2*mu^2*sigma-diff(mu,x,2)-2*mu*(diff(mu,x,1))),
184
185      eq : partfrac(eq,x),
186      eq : if safe_op(eq) = "+" then args(eq) else [eq],
187      eq : xreduce('append, map(lambda([s], all_poly_coefs(ratnumer(s),x)),eq)),
188      eq : listify(setify(eq)),
189      sol : algsys(eq, [%k]),
190      for sk in sol while ie = false do (
191        f : radcan(exp(logcontract(integrate(mu,x)))),
192        f : block([?errorsw : true], errcatch(subst(sk, f))),
193        if f #  [] then (
194          f : first(f),
195          %fo : new_variable('general),
196          f : %fo * f,
197          cnst : subst(sk, (diff(f,x) / (sqrt(1-f^2) * (sqrt(1-%k^2 * f^2)))) / e),
198          cnst : infapply(cnst,lambda([s], radcan(rootscontract(s)))),
199          if cnst # 0 then (
200            cnst : ratnumer(ratsimp(diff(cnst,x) / cnst)),
201            cnst : algsys([cnst], [%fo]),
202            cnst : first(sublist(cnst, lambda([s], freeof(x, rhs(first(s)))))))
203          else cnst : [ ],
204
205          if cnst # [] then (
206            f : subst(cnst,f),
207            %so : ratsimp(subst(sk, diff(f,x) / (sqrt(1-f^2) * (sqrt(1-%k^2 * f^2))) / e)),
208            if %so # 0 and piecewise_constant_p(%so,x) = true then (
209              ie : subst(sk, (1/%so) * inverse_jacobi_sn(f, %k^2)))))))),
210    ie);
211