Source code

001/*
002 * $Id: ElementaryIntegration.java 4125 2012-08-19 19:05:22Z kredel $
003 */
004
005package edu.jas.integrate;
006
007
008import java.util.ArrayList;
009import java.util.List;
010import java.util.SortedMap;
011
012import org.apache.log4j.Logger;
013
014import edu.jas.poly.AlgebraicNumber;
015import edu.jas.poly.AlgebraicNumberRing;
016import edu.jas.poly.GenPolynomial;
017import edu.jas.poly.GenPolynomialRing;
018import edu.jas.poly.PolyUtil;
019import edu.jas.structure.GcdRingElem;
020import edu.jas.structure.Power;
021import edu.jas.structure.RingFactory;
022import edu.jas.ufd.FactorAbstract;
023import edu.jas.ufd.FactorFactory;
024import edu.jas.ufd.GCDFactory;
025import edu.jas.ufd.GreatestCommonDivisorAbstract;
026import edu.jas.ufd.GreatestCommonDivisorSubres;
027import edu.jas.ufd.PolyUfdUtil;
028import edu.jas.ufd.Quotient;
029import edu.jas.ufd.QuotientRing;
030import edu.jas.ufd.SquarefreeAbstract;
031import 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
044public 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.recursiveUnivariateResultant(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.copy();
381            String[] unused = 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}