001 /*
002 * $Id: ElementaryIntegration.java 3356 2010-10-23 16:41:01Z kredel $
003 */
004
005 package edu.jas.integrate;
006
007
008 import java.util.ArrayList;
009 import java.util.List;
010 import java.util.SortedMap;
011
012 import org.apache.log4j.Logger;
013
014 import edu.jas.poly.AlgebraicNumber;
015 import edu.jas.poly.AlgebraicNumberRing;
016 import edu.jas.poly.GenPolynomial;
017 import edu.jas.poly.GenPolynomialRing;
018 import edu.jas.poly.PolyUtil;
019 import edu.jas.structure.GcdRingElem;
020 import edu.jas.structure.Power;
021 import edu.jas.structure.RingFactory;
022 import edu.jas.ufd.FactorAbstract;
023 import edu.jas.ufd.FactorFactory;
024 import edu.jas.ufd.GCDFactory;
025 import edu.jas.ufd.GreatestCommonDivisorAbstract;
026 import edu.jas.ufd.GreatestCommonDivisorSubres;
027 import edu.jas.ufd.PolyUfdUtil;
028 import edu.jas.ufd.Quotient;
029 import edu.jas.ufd.QuotientRing;
030 import edu.jas.ufd.SquarefreeAbstract;
031 import edu.jas.ufd.SquarefreeFactory;
032
033
034 /**
035 * Methods related to elementary integration. In particular there are methods
036 * for Hermite reduction and Rothstein-Trager integration of the logarithmic
037 * part.
038 *
039 * @author Axel Kramer
040 * @author Heinz Kredel
041 * @param <C> coefficient type
042 */
043
044 public class ElementaryIntegration<C extends GcdRingElem<C>> {
045
046
047 private static final Logger logger = Logger.getLogger(ElementaryIntegration.class);
048
049
050 private final boolean debug = logger.isDebugEnabled();
051
052
053 /**
054 * Engine for factorization.
055 */
056 public final FactorAbstract<C> irr;
057
058
059 /**
060 * Engine for squarefree decomposition.
061 */
062 public final SquarefreeAbstract<C> sqf;
063
064
065 /**
066 * Engine for greatest common divisors.
067 */
068 public final GreatestCommonDivisorAbstract<C> ufd;
069
070
071 /**
072 * Constructor.
073 */
074 public ElementaryIntegration(RingFactory<C> br) {
075 ufd = GCDFactory.<C> getProxy(br);
076 sqf = SquarefreeFactory.<C> getImplementation(br);
077 irr = /*(FactorAbsolute<C>)*/FactorFactory.<C> getImplementation(br);
078 }
079
080
081 /**
082 * Integration of a rational function.
083 *
084 * @param r rational function
085 * @return Integral container, such that integrate(r) = sum_i(g_i) + sum_j(
086 * an_j log(hd_j) )
087 */
088 public QuotIntegral<C> integrate(Quotient<C> r) {
089 Integral<C> integral = integrate(r.num, r.den);
090 return new QuotIntegral<C>(r.ring, integral);
091 }
092
093
094 /**
095 * Integration of a rational function.
096 *
097 * @param a numerator
098 * @param d denominator
099 * @return Integral container, such that integrate(a/d) = sum_i(gn_i/gd_i) +
100 * integrate(h0) + sum_j( an_j log(hd_j) )
101 */
102 public Integral<C> integrate(GenPolynomial<C> a, GenPolynomial<C> d) {
103 if (d == null || a == null || d.isZERO()) {
104 throw new IllegalArgumentException("zero or null not allowed");
105 }
106 if (a.isZERO()) {
107 return new Integral<C>(a, d, a);
108 }
109 if (d.isONE()) {
110 GenPolynomial<C> pi = PolyUtil.<C> baseIntegral(a);
111 return new Integral<C>(a, d, pi);
112 }
113 GenPolynomialRing<C> pfac = d.ring;
114 if (pfac.nvar > 1) {
115 throw new IllegalArgumentException("only for univariate polynomials " + pfac);
116 }
117 if (!pfac.coFac.isField()) {
118 throw new IllegalArgumentException("only for field coefficients " + pfac);
119 }
120
121 GenPolynomial<C>[] qr = PolyUtil.<C> basePseudoQuotientRemainder(a, d);
122 GenPolynomial<C> p = qr[0];
123 GenPolynomial<C> r = qr[1];
124
125 GenPolynomial<C> c = ufd.gcd(r, d);
126 if (!c.isONE()) {
127 r = PolyUtil.<C> basePseudoQuotientRemainder(r, c)[0];
128 d = PolyUtil.<C> basePseudoQuotientRemainder(d, c)[0];
129 }
130 List<GenPolynomial<C>>[] ih = integrateHermite(r, d);
131 List<GenPolynomial<C>> rat = ih[0];
132 List<GenPolynomial<C>> log = ih[1];
133
134 GenPolynomial<C> pp = log.remove(0);
135 p = p.sum(pp);
136 GenPolynomial<C> pi = PolyUtil.<C> baseIntegral(p);
137
138 if (debug) {
139 logger.debug("pi = " + pi);
140 logger.debug("rat = " + rat);
141 logger.debug("log = " + log);
142 }
143 if (log.size() == 0) {
144 return new Integral<C>(a, d, pi, rat);
145 }
146
147 List<LogIntegral<C>> logi = new ArrayList<LogIntegral<C>>(log.size() / 2);
148 for (int i = 0; i < log.size(); i++) {
149 GenPolynomial<C> ln = log.get(i++);
150 GenPolynomial<C> ld = log.get(i);
151 LogIntegral<C> pf = integrateLogPart(ln, ld);
152 logi.add(pf);
153 }
154 if (debug) {
155 logger.debug("logi = " + logi);
156 }
157 return new Integral<C>(a, d, pi, rat, logi);
158 }
159
160
161 /**
162 * Integration of the rational part, Hermite reduction step.
163 *
164 * @param a numerator
165 * @param d denominator, gcd(a,d) == 1
166 * @return [ [ gn_i, gd_i ], [ h0, hn_j, hd_j ] ] such that integrate(a/d) =
167 * sum_i(gn_i/gd_i) + integrate(h0) + sum_j( integrate(hn_j/hd_j) )
168 */
169 public List<GenPolynomial<C>>[] integrateHermite(GenPolynomial<C> a, GenPolynomial<C> d) {
170 if (d == null || d.isZERO()) {
171 throw new IllegalArgumentException("d == null or d == 0");
172 }
173 if (a == null || a.isZERO()) {
174 throw new IllegalArgumentException("a == null or a == 0");
175 }
176
177 // get squarefree decomposition
178 SortedMap<GenPolynomial<C>, Long> sfactors = sqf.squarefreeFactors(d);
179
180 List<GenPolynomial<C>> D = new ArrayList<GenPolynomial<C>>(sfactors.keySet());
181 List<GenPolynomial<C>> DP = new ArrayList<GenPolynomial<C>>();
182 for (GenPolynomial<C> f : D) {
183 long e = sfactors.get(f);
184 GenPolynomial<C> dp = Power.<GenPolynomial<C>> positivePower(f, e);
185 DP.add(dp);
186 }
187 //System.out.println("D: " + D);
188 //System.out.println("DP: " + DP);
189
190 // get partial fraction decompostion
191 List<GenPolynomial<C>> Ai = ufd.basePartialFraction(a, DP);
192 //System.out.println("Ai: " + Ai);
193
194 List<GenPolynomial<C>> G = new ArrayList<GenPolynomial<C>>();
195 List<GenPolynomial<C>> H = new ArrayList<GenPolynomial<C>>();
196 H.add(Ai.remove(0)); // P
197
198 GenPolynomialRing<C> fac = d.ring;
199 int i = 0;
200 for (GenPolynomial<C> v : D) {
201 //System.out.println("V:" + v.toString());
202 GenPolynomial<C> Ak = Ai.get(i++);
203 //System.out.println("Ak: " + Ak.toString());
204 int k = sfactors.get(v).intValue(); // assert low power
205 for (int j = k - 1; j >= 1; j--) {
206 //System.out.println("Step(" + k + "," + j + ")");
207 GenPolynomial<C> DV_dx = PolyUtil.<C> baseDeriviative(v);
208 GenPolynomial<C> Aik = Ak.divide(fac.fromInteger(-j));
209 GenPolynomial<C>[] BC = ufd.baseGcdDiophant(DV_dx, v, Aik);
210 GenPolynomial<C> b = BC[0];
211 GenPolynomial<C> c = BC[1];
212 GenPolynomial<C> vj = Power.<GenPolynomial<C>> positivePower(v, j);
213 G.add(b); // B
214 G.add(vj); // v^j
215 Ak = fac.fromInteger(-j).multiply(c).subtract(PolyUtil.<C> baseDeriviative(b));
216 //System.out.println("B: " + b.toString());
217 //System.out.println("C: " + c.toString());
218 }
219 //System.out.println("V:" + v.toString());
220 //System.out.println("Ak: " + Ak.toString());
221 if (!Ak.isZERO()) {
222 H.add(Ak); // A_k
223 H.add(v); // v
224 }
225 }
226 List<GenPolynomial<C>>[] ret = (List<GenPolynomial<C>>[]) new List[2];
227 ret[0] = G;
228 ret[1] = H;
229 return ret;
230 }
231
232
233 /**
234 * Univariate GenPolynomial integration of the logaritmic part,
235 * Rothstein-Trager algorithm.
236 * @param A univariate GenPolynomial, deg(A) < deg(P).
237 * @param P univariate squarefree GenPolynomial, gcd(A,P) == 1.
238 * @return logarithmic part container.
239 */
240 public LogIntegral<C> integrateLogPart(GenPolynomial<C> A, GenPolynomial<C> P) {
241 if (P == null || P.isZERO()) {
242 throw new IllegalArgumentException(" P == null or P == 0");
243 }
244 if (A == null || A.isZERO()) {
245 throw new IllegalArgumentException(" A == null or A == 0");
246 }
247 //System.out.println("\nP_base_algeb_part = " + P);
248 GenPolynomialRing<C> pfac = P.ring; // K[x]
249 if (pfac.nvar > 1) {
250 throw new IllegalArgumentException("only for univariate polynomials " + pfac);
251 }
252 if (!pfac.coFac.isField()) {
253 throw new IllegalArgumentException("only for field coefficients " + pfac);
254 }
255 List<C> cfactors = new ArrayList<C>();
256 List<GenPolynomial<C>> cdenom = new ArrayList<GenPolynomial<C>>();
257 List<AlgebraicNumber<C>> afactors = new ArrayList<AlgebraicNumber<C>>();
258 List<GenPolynomial<AlgebraicNumber<C>>> adenom = new ArrayList<GenPolynomial<AlgebraicNumber<C>>>();
259
260 // P linear
261 if (P.degree(0) <= 1) {
262 cfactors.add(A.leadingBaseCoefficient());
263 cdenom.add(P);
264 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom);
265 }
266 List<GenPolynomial<C>> Pfac = irr.baseFactorsSquarefree(P);
267 //System.out.println("\nPfac = " + Pfac);
268
269 List<GenPolynomial<C>> Afac = ufd.basePartialFraction(A, Pfac);
270
271 GenPolynomial<C> A0 = Afac.remove(0);
272 if (!A0.isZERO()) {
273 throw new RuntimeException(" A0 != 0: deg(A)>= deg(P)");
274 }
275
276 // algebraic and linear factors
277 int i = 0;
278 for (GenPolynomial<C> pi : Pfac) {
279 GenPolynomial<C> ai = Afac.get(i++);
280 if (pi.degree(0) <= 1) {
281 cfactors.add(ai.leadingBaseCoefficient());
282 cdenom.add(pi);
283 continue;
284 }
285 LogIntegral<C> pf = integrateLogPartIrreducible(ai, pi);
286 cfactors.addAll(pf.cfactors);
287 cdenom.addAll(pf.cdenom);
288 afactors.addAll(pf.afactors);
289 adenom.addAll(pf.adenom);
290 }
291 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom);
292 }
293
294
295 /**
296 * Univariate GenPolynomial integration of the logaritmic part,
297 * Rothstein-Trager algorithm.
298 * @param A univariate GenPolynomial, deg(A) < deg(P).
299 * @param P univariate irreducible GenPolynomial. // gcd(A,P) == 1 automatic
300 * @return logarithmic part container.
301 */
302 public LogIntegral<C> integrateLogPartIrreducible(GenPolynomial<C> A, GenPolynomial<C> P) {
303 if (P == null || P.isZERO()) {
304 throw new IllegalArgumentException("P == null or P == 0");
305 }
306 //System.out.println("\nP_base_algeb_part = " + P);
307 GenPolynomialRing<C> pfac = P.ring; // K[x]
308 if (pfac.nvar > 1) {
309 throw new IllegalArgumentException("only for univariate polynomials " + pfac);
310 }
311 if (!pfac.coFac.isField()) {
312 throw new IllegalArgumentException("only for field coefficients " + pfac);
313 }
314 List<C> cfactors = new ArrayList<C>();
315 List<GenPolynomial<C>> cdenom = new ArrayList<GenPolynomial<C>>();
316 List<AlgebraicNumber<C>> afactors = new ArrayList<AlgebraicNumber<C>>();
317 List<GenPolynomial<AlgebraicNumber<C>>> adenom = new ArrayList<GenPolynomial<AlgebraicNumber<C>>>();
318
319 // P linear
320 if (P.degree(0) <= 1) {
321 cfactors.add(A.leadingBaseCoefficient());
322 cdenom.add(P);
323 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom);
324 }
325
326 // deriviative
327 GenPolynomial<C> Pp = PolyUtil.<C> baseDeriviative(P);
328 //no: Pp = Pp.monic();
329 //System.out.println("\nP = " + P);
330 //System.out.println("Pp = " + Pp);
331
332 // Q[t]
333 String[] vars = new String[] { "t" };
334 GenPolynomialRing<C> cfac = new GenPolynomialRing<C>(pfac.coFac, 1, pfac.tord, vars);
335 GenPolynomial<C> t = cfac.univariate(0);
336 //System.out.println("t = " + t);
337
338 // Q[x][t]
339 GenPolynomialRing<GenPolynomial<C>> rfac = new GenPolynomialRing<GenPolynomial<C>>(pfac, cfac); // sic
340 //System.out.println("rfac = " + rfac.toScript());
341
342 // transform polynomials to bi-variate polynomial
343 GenPolynomial<GenPolynomial<C>> Ac = PolyUfdUtil.<C> introduceLowerVariable(rfac, A);
344 //System.out.println("Ac = " + Ac);
345 GenPolynomial<GenPolynomial<C>> Pc = PolyUfdUtil.<C> introduceLowerVariable(rfac, P);
346 //System.out.println("Pc = " + Pc);
347 GenPolynomial<GenPolynomial<C>> Pcp = PolyUfdUtil.<C> introduceLowerVariable(rfac, Pp);
348 //System.out.println("Pcp = " + Pcp);
349
350 // Q[t][x]
351 GenPolynomialRing<GenPolynomial<C>> rfac1 = Pc.ring;
352 //System.out.println("rfac1 = " + rfac1.toScript());
353
354 // A - t P'
355 GenPolynomial<GenPolynomial<C>> tc = rfac1.getONE().multiply(t);
356 //System.out.println("tc = " + tc);
357 GenPolynomial<GenPolynomial<C>> At = Ac.subtract(tc.multiply(Pcp));
358 //System.out.println("At = " + At);
359
360 GreatestCommonDivisorSubres<C> engine = new GreatestCommonDivisorSubres<C>();
361 // = GCDFactory.<C>getImplementation( cfac.coFac );
362 GreatestCommonDivisorAbstract<AlgebraicNumber<C>> aengine = null;
363
364 GenPolynomial<GenPolynomial<C>> Rc = engine.recursiveResultant(Pc, At);
365 //System.out.println("Rc = " + Rc);
366 GenPolynomial<C> res = Rc.leadingBaseCoefficient();
367 //no: res = res.monic();
368 //System.out.println("\nres = " + res);
369
370 SortedMap<GenPolynomial<C>, Long> resfac = irr.baseFactors(res);
371 //System.out.println("resfac = " + resfac + "\n");
372
373 for (GenPolynomial<C> r : resfac.keySet()) {
374 //System.out.println("\nr(t) = " + r);
375 if (r.isConstant()) {
376 continue;
377 }
378 //vars = new String[] { "z_" + Math.abs(r.hashCode() % 1000) };
379 vars = pfac.newVars("z_");
380 pfac = pfac.clone();
381 vars = pfac.setVars(vars);
382 r = pfac.copy(r); // hack to exchange the variables
383 //System.out.println("r(z_) = " + r);
384 AlgebraicNumberRing<C> afac = new AlgebraicNumberRing<C>(r, true); // since irreducible
385 logger.debug("afac = " + afac.toScript());
386 AlgebraicNumber<C> a = afac.getGenerator();
387 //no: a = a.negate();
388 //System.out.println("a = " + a);
389
390 // K(alpha)[x]
391 GenPolynomialRing<AlgebraicNumber<C>> pafac = new GenPolynomialRing<AlgebraicNumber<C>>(afac,
392 Pc.ring);
393 //System.out.println("pafac = " + pafac.toScript());
394
395 // convert to K(alpha)[x]
396 GenPolynomial<AlgebraicNumber<C>> Pa = PolyUtil.<C> convertToAlgebraicCoefficients(pafac, P);
397 //System.out.println("Pa = " + Pa);
398 GenPolynomial<AlgebraicNumber<C>> Pap = PolyUtil.<C> convertToAlgebraicCoefficients(pafac, Pp);
399 //System.out.println("Pap = " + Pap);
400 GenPolynomial<AlgebraicNumber<C>> Aa = PolyUtil.<C> convertToAlgebraicCoefficients(pafac, A);
401 //System.out.println("Aa = " + Aa);
402
403 // A - a P'
404 GenPolynomial<AlgebraicNumber<C>> Ap = Aa.subtract(Pap.multiply(a));
405 //System.out.println("Ap = " + Ap);
406
407 if (aengine == null) {
408 aengine = GCDFactory.<AlgebraicNumber<C>> getImplementation(afac);
409 }
410 GenPolynomial<AlgebraicNumber<C>> Ga = aengine.baseGcd(Pa, Ap);
411 //System.out.println("Ga = " + Ga);
412 if (Ga.isConstant()) {
413 //System.out.println("warning constant gcd ignored");
414 continue;
415 }
416 afactors.add(a);
417 adenom.add(Ga);
418 // special quadratic case
419 if (P.degree(0) == 2 && Ga.degree(0) == 1) {
420 GenPolynomial<AlgebraicNumber<C>>[] qra = PolyUtil
421 .<AlgebraicNumber<C>> basePseudoQuotientRemainder(Pa, Ga);
422 GenPolynomial<AlgebraicNumber<C>> Qa = qra[0];
423 if (!qra[1].isZERO()) {
424 throw new ArithmeticException("remainder not zero");
425 }
426 //System.out.println("Qa = " + Qa);
427 afactors.add(a.negate());
428 adenom.add(Qa);
429 }
430 // todo: eventually implement special cases deg = 3, 4
431 }
432 return new LogIntegral<C>(A, P, cfactors, cdenom, afactors, adenom);
433 }
434
435
436 /**
437 * Derivation of a univariate rational function.
438 *
439 * @param r rational function
440 * @return dr/dx
441 */
442 public Quotient<C> deriviative(Quotient<C> r) {
443 GenPolynomial<C> num = r.num;
444 GenPolynomial<C> den = r.den;
445 GenPolynomial<C> nump = PolyUtil.<C> baseDeriviative(num);
446 if (den.isONE()) {
447 return new Quotient<C>(r.ring, nump, den);
448 }
449 GenPolynomial<C> denp = PolyUtil.<C> baseDeriviative(den);
450
451 GenPolynomial<C> n = den.multiply(nump).subtract(num.multiply(denp));
452 GenPolynomial<C> d = den.multiply(den);
453
454 Quotient<C> der = new Quotient<C>(r.ring, n, d);
455 return der;
456 }
457
458
459 /**
460 * Test of integration of a rational function.
461 *
462 * @param ri integral
463 * @return true, if ri is an integral, else false.
464 */
465 public boolean isIntegral(QuotIntegral<C> ri) {
466 Quotient<C> r = ri.quot;
467 QuotientRing<C> qr = r.ring;
468 Quotient<C> i = r.ring.getZERO();
469 for (Quotient<C> q : ri.rational) {
470 Quotient<C> qd = deriviative(q);
471 i = i.sum(qd);
472 }
473 if (ri.logarithm.size() == 0) {
474 return r.equals(i);
475 }
476 for (LogIntegral<C> li : ri.logarithm) {
477 Quotient<C> q = new Quotient<C>(qr, li.num, li.den);
478 i = i.sum(q);
479 }
480 boolean t = r.equals(i);
481 if (!t) {
482 return false;
483 }
484 for (LogIntegral<C> li : ri.logarithm) {
485 t = isIntegral(li);
486 if (!t) {
487 return false;
488 }
489 }
490 return true;
491 }
492
493
494 /**
495 * Test of integration of the logarithmic part of a rational function.
496 *
497 * @param rl logarithmic part of an integral
498 * @return true, if rl is an integral, else false.
499 */
500 public boolean isIntegral(LogIntegral<C> rl) {
501 QuotientRing<C> qr = new QuotientRing<C>(rl.den.ring);
502 Quotient<C> r = new Quotient<C>(qr, rl.num, rl.den);
503
504 Quotient<C> i = qr.getZERO();
505 int j = 0;
506 for (GenPolynomial<C> d : rl.cdenom) {
507 GenPolynomial<C> dp = PolyUtil.<C> baseDeriviative(d);
508 dp = dp.multiply(rl.cfactors.get(j++));
509 Quotient<C> f = new Quotient<C>(qr, dp, d);
510 i = i.sum(f);
511 }
512 if (rl.afactors.size() == 0) {
513 return r.equals(i);
514 }
515 r = r.subtract(i);
516 QuotientRing<AlgebraicNumber<C>> aqr = new QuotientRing<AlgebraicNumber<C>>(rl.adenom.get(0).ring);
517 Quotient<AlgebraicNumber<C>> ai = aqr.getZERO();
518
519 GenPolynomial<AlgebraicNumber<C>> aqn = PolyUtil.<C> convertToAlgebraicCoefficients(aqr.ring, r.num);
520 GenPolynomial<AlgebraicNumber<C>> aqd = PolyUtil.<C> convertToAlgebraicCoefficients(aqr.ring, r.den);
521 Quotient<AlgebraicNumber<C>> ar = new Quotient<AlgebraicNumber<C>>(aqr, aqn, aqd);
522
523 j = 0;
524 for (GenPolynomial<AlgebraicNumber<C>> d : rl.adenom) {
525 GenPolynomial<AlgebraicNumber<C>> dp = PolyUtil.<AlgebraicNumber<C>> baseDeriviative(d);
526 dp = dp.multiply(rl.afactors.get(j++));
527 Quotient<AlgebraicNumber<C>> f = new Quotient<AlgebraicNumber<C>>(aqr, dp, d);
528 ai = ai.sum(f);
529 }
530 boolean t = ar.equals(ai);
531 if (t) {
532 return true;
533 }
534 logger.warn("log integral not verified");
535 //System.out.println("r = " + r);
536 //System.out.println("afactors = " + rl.afactors);
537 //System.out.println("adenom = " + rl.adenom);
538 //System.out.println("ar = " + ar);
539 //System.out.println("ai = " + ai);
540 return true;
541 }
542
543 }