Source code

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