Source code

001/*
002 * $Id: GreatestCommonDivisorSubres.java 5918 2018-08-29 20:32:55Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.Logger;
012import org.apache.logging.log4j.LogManager; 
013
014import edu.jas.poly.GenPolynomial;
015import edu.jas.poly.GenPolynomialRing;
016import edu.jas.poly.PolyUtil;
017import edu.jas.structure.GcdRingElem;
018import edu.jas.structure.RingFactory;
019
020
021/**
022 * Greatest common divisor algorithms with subresultant polynomial remainder
023 * sequence.
024 * @author Heinz Kredel
025 * @author Youssef Elbarbary
026 */
027
028public class GreatestCommonDivisorSubres<C extends GcdRingElem<C>> extends GreatestCommonDivisorAbstract<C> {
029
030
031    private static final Logger logger = LogManager.getLogger(GreatestCommonDivisorSubres.class);
032
033
034    private static final boolean debug = logger.isDebugEnabled();
035
036
037    /**
038     * GenPolynomial pseudo remainder. For univariate polynomials.
039     * @param P GenPolynomial.
040     * @param S nonzero GenPolynomial.
041     * @return remainder with ldcf(S)<sup>m</sup> P = quotient * S + remainder.
042     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
043     * @deprecated Use
044     *             {@link edu.jas.poly.PolyUtil#baseDensePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)}
045     *             instead
046     */
047    @Deprecated
048    public GenPolynomial<C> basePseudoRemainder(GenPolynomial<C> P, GenPolynomial<C> S) {
049        return PolyUtil.<C> baseDensePseudoRemainder(P, S);
050    }
051
052
053    /**
054     * GenPolynomial pseudo remainder. For recursive polynomials.
055     * @param P recursive GenPolynomial.
056     * @param S nonzero recursive GenPolynomial.
057     * @return remainder with ldcf(S)<sup>m</sup> P = quotient * S + remainder.
058     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
059     * @deprecated Use
060     *             {@link edu.jas.poly.PolyUtil#recursiveDensePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)}
061     *             instead
062     */
063    @Deprecated
064    public GenPolynomial<GenPolynomial<C>> recursivePseudoRemainder(GenPolynomial<GenPolynomial<C>> P,
065                    GenPolynomial<GenPolynomial<C>> S) {
066        return PolyUtil.<C> recursiveDensePseudoRemainder(P, S);
067    }
068
069
070    /**
071     * Univariate GenPolynomial greatest comon divisor. Uses pseudoRemainder for
072     * remainder.
073     * @param P univariate GenPolynomial.
074     * @param S univariate GenPolynomial.
075     * @return gcd(P,S).
076     */
077    @Override
078    public GenPolynomial<C> baseGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
079        if (S == null || S.isZERO()) {
080            return P;
081        }
082        if (P == null || P.isZERO()) {
083            return S;
084        }
085        if (P.ring.nvar > 1) {
086            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
087        }
088        long e = P.degree(0);
089        long f = S.degree(0);
090        GenPolynomial<C> q;
091        GenPolynomial<C> r;
092        if (f > e) {
093            r = P;
094            q = S;
095            long g = f;
096            f = e;
097            e = g;
098        } else {
099            q = P;
100            r = S;
101        }
102        if (debug) {
103            logger.debug("degrees: e = " + e + ", f = " + f);
104        }
105        r = r.abs();
106        q = q.abs();
107        C a = baseContent(r);
108        C b = baseContent(q);
109        C c = gcd(a, b); // indirection
110        r = divide(r, a); // indirection
111        q = divide(q, b); // indirection
112        if (r.isONE()) {
113            return r.multiply(c);
114        }
115        if (q.isONE()) {
116            return q.multiply(c);
117        }
118        C g = r.ring.getONECoefficient();
119        C h = r.ring.getONECoefficient();
120        GenPolynomial<C> x;
121        C z;
122        while (!r.isZERO()) {
123            long delta = q.degree(0) - r.degree(0);
124            //System.out.println("delta    = " + delta);
125            x = PolyUtil.<C> baseDensePseudoRemainder(q, r);
126            q = r;
127            if (!x.isZERO()) {
128                z = g.multiply(h.power(delta)); //power(P.ring.coFac, h, delta));
129                //System.out.println("z  = " + z);
130                r = x.divide(z);
131                g = q.leadingBaseCoefficient();
132                z = g.power(delta); //power(P.ring.coFac, g, delta);
133                h = z.divide(h.power(delta - 1)); //power(P.ring.coFac, h, delta - 1));
134                //System.out.println("g  = " + g);
135                //System.out.println("h  = " + h);
136            } else {
137                r = x;
138            }
139        }
140        q = basePrimitivePart(q);
141        return (q.multiply(c)).abs();
142    }
143
144
145    /**
146     * Univariate GenPolynomial recursive greatest comon divisor. Uses
147     * pseudoRemainder for remainder.
148     * @param P univariate recursive GenPolynomial.
149     * @param S univariate recursive GenPolynomial.
150     * @return gcd(P,S).
151     */
152    @Override
153    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateGcd(GenPolynomial<GenPolynomial<C>> P,
154                    GenPolynomial<GenPolynomial<C>> S) {
155        if (S == null || S.isZERO()) {
156            return P;
157        }
158        if (P == null || P.isZERO()) {
159            return S;
160        }
161        if (P.ring.nvar > 1) {
162            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
163        }
164        long e = P.degree(0);
165        long f = S.degree(0);
166        GenPolynomial<GenPolynomial<C>> q;
167        GenPolynomial<GenPolynomial<C>> r;
168        if (f > e) {
169            r = P;
170            q = S;
171            long g = f;
172            f = e;
173            e = g;
174        } else {
175            q = P;
176            r = S;
177        }
178        if (debug) {
179            logger.debug("degrees: e = " + e + ", f = " + f);
180        }
181        r = r.abs();
182        q = q.abs();
183        GenPolynomial<C> a = recursiveContent(r);
184        GenPolynomial<C> b = recursiveContent(q);
185
186        GenPolynomial<C> c = gcd(a, b); // go to recursion
187        //System.out.println("rgcd c = " + c);
188        r = PolyUtil.<C> recursiveDivide(r, a);
189        q = PolyUtil.<C> recursiveDivide(q, b);
190        if (r.isONE()) {
191            return r.multiply(c);
192        }
193        if (q.isONE()) {
194            return q.multiply(c);
195        }
196        GenPolynomial<C> g = r.ring.getONECoefficient();
197        GenPolynomial<C> h = r.ring.getONECoefficient();
198        GenPolynomial<GenPolynomial<C>> x;
199        GenPolynomial<C> z = null;
200        while (!r.isZERO()) {
201            long delta = q.degree(0) - r.degree(0);
202            //System.out.println("rgcd delta = " + delta);
203            x = PolyUtil.<C> recursiveDensePseudoRemainder(q, r);
204            if (logger.isDebugEnabled()) {
205                logger.info("recursiveDensePseudoRemainder.bits = " + x.bitLength());
206            }
207            q = r;
208            if (!x.isZERO()) {
209                z = g.multiply(h.power(delta)); //power(P.ring.coFac, h, delta));
210                r = PolyUtil.<C> recursiveDivide(x, z);
211                g = q.leadingBaseCoefficient();
212                z = g.power(delta); //power(P.ring.coFac, g, delta);
213                h = PolyUtil.<C> basePseudoDivide(z, h.power(delta - 1)); // power(P.ring.coFac, h, delta - 1)
214            } else {
215                r = x;
216            }
217        }
218        q = recursivePrimitivePart(q);
219        return q.abs().multiply(c); //.abs();
220    }
221
222
223    /**
224     * Univariate GenPolynomial resultant. Uses pseudoRemainder for remainder.
225     * @param P univariate GenPolynomial.
226     * @param S univariate GenPolynomial.
227     * @return res(P,S).
228     */
229    @Override
230    public GenPolynomial<C> baseResultant(GenPolynomial<C> P, GenPolynomial<C> S) {
231        if (S == null || S.isZERO()) {
232            return S;
233        }
234        if (P == null || P.isZERO()) {
235            return P;
236        }
237        if (P.ring.nvar > 1) {
238            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
239        }
240        long e = P.degree(0);
241        long f = S.degree(0);
242        GenPolynomial<C> q;
243        GenPolynomial<C> r;
244        if (f > e) {
245            r = P;
246            q = S;
247            long g = f;
248            f = e;
249            e = g;
250        } else {
251            q = P;
252            r = S;
253        }
254        //r = r.abs();
255        //q = q.abs();
256        C a = baseContent(r);
257        C b = baseContent(q);
258        r = divide(r, a); // indirection
259        q = divide(q, b); // indirection
260        RingFactory<C> cofac = P.ring.coFac;
261        C g = cofac.getONE();
262        C h = cofac.getONE();
263        C t = a.power(e); //power(cofac, a, e);
264        t = t.multiply(b.power(f)); //power(cofac, b, f));
265        long s = 1;
266        GenPolynomial<C> x;
267        C z;
268        while (r.degree(0) > 0) {
269            long delta = q.degree(0) - r.degree(0);
270            //System.out.println("delta    = " + delta);
271            if ((q.degree(0) % 2 != 0) && (r.degree(0) % 2 != 0)) {
272                s = -s;
273            }
274            x = PolyUtil.<C> baseDensePseudoRemainder(q, r);
275            //System.out.println("x  = " + x);
276            q = r;
277            if (x.degree(0) > 0) {
278                z = g.multiply(h.power(delta)); //power(cofac, h, delta));
279                //System.out.println("z  = " + z);
280                r = x.divide(z);
281                g = q.leadingBaseCoefficient();
282                z = g.power(delta); //power(cofac, g, delta);
283                h = z.divide(h.power(delta - 1));
284            } else {
285                r = x;
286            }
287        }
288        z = r.leadingBaseCoefficient().power(q.degree(0)); //power(cofac, r.leadingBaseCoefficient(), q.degree(0));
289        h = z.divide(h.power(q.degree() - 1)); //power(cofac, h, q.degree(0) - 1));
290        z = h.multiply(t);
291        if (s < 0) {
292            z = z.negate();
293        }
294        x = P.ring.getONE().multiply(z);
295        return x;
296    }
297
298
299    /**
300     * Univariate GenPolynomial recursive resultant. Uses pseudoRemainder for
301     * remainder.
302     * @param P univariate recursive GenPolynomial.
303     * @param S univariate recursive GenPolynomial.
304     * @return res(P,S).
305     */
306    @Override
307    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<C>> P,
308                    GenPolynomial<GenPolynomial<C>> S) {
309        if (S == null || S.isZERO()) {
310            return S;
311        }
312        if (P == null || P.isZERO()) {
313            return P;
314        }
315        if (P.ring.nvar > 1) {
316            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
317        }
318        long e = P.degree(0);
319        long f = S.degree(0);
320        GenPolynomial<GenPolynomial<C>> q;
321        GenPolynomial<GenPolynomial<C>> r;
322        if (f > e) {
323            r = P;
324            q = S;
325            long g = f;
326            f = e;
327            e = g;
328        } else {
329            q = P;
330            r = S;
331        }
332        r = r.abs();
333        q = q.abs();
334        GenPolynomial<C> a = recursiveContent(r);
335        GenPolynomial<C> b = recursiveContent(q);
336        r = PolyUtil.<C> recursiveDivide(r, a);
337        q = PolyUtil.<C> recursiveDivide(q, b);
338        RingFactory<GenPolynomial<C>> cofac = P.ring.coFac;
339        GenPolynomial<C> g = cofac.getONE();
340        GenPolynomial<C> h = cofac.getONE();
341        GenPolynomial<GenPolynomial<C>> x;
342        GenPolynomial<C> t;
343        if (f == 0 && e == 0 && g.ring.nvar > 0) {
344            // if coeffs are multivariate (and non constant)
345            // otherwise it would be 1
346            t = resultant(a, b);
347            x = P.ring.getONE().multiply(t);
348            return x;
349        }
350        t = a.power(e); //power(cofac, a, e);
351        t = t.multiply(b.power(f)); //power(cofac, b, f));
352        long s = 1;
353        GenPolynomial<C> z;
354        while (r.degree(0) > 0) {
355            long delta = q.degree(0) - r.degree(0);
356            //System.out.println("delta    = " + delta);
357            if ((q.degree(0) % 2 != 0) && (r.degree(0) % 2 != 0)) {
358                s = -s;
359            }
360            x = PolyUtil.<C> recursiveDensePseudoRemainder(q, r);
361            //System.out.println("x  = " + x);
362            q = r;
363            if (x.degree(0) > 0) {
364                z = g.multiply(h.power(delta)); //power(P.ring.coFac, h, delta));
365                r = PolyUtil.<C> recursiveDivide(x, z);
366                g = q.leadingBaseCoefficient();
367                z = g.power(delta); //power(cofac, g, delta);
368                h = PolyUtil.<C> basePseudoDivide(z, h.power(delta - 1)); //power(cofac, h, delta - 1));
369            } else {
370                r = x;
371            }
372        }
373        z = r.leadingBaseCoefficient().power(q.degree(0)); //power(cofac, r.leadingBaseCoefficient(), q.degree(0));
374        h = PolyUtil.<C> basePseudoDivide(z, h.power(q.degree() - 1)); //power(cofac, h, q.degree(0) - 1));
375        z = h.multiply(t);
376        if (s < 0) {
377            z = z.negate();
378        }
379        x = P.ring.getONE().multiply(z);
380        return x;
381    }
382
383
384    /**
385     * Univariate GenPolynomial recursive Subresultant list. Uses
386     * pseudoRemainder for remainder.
387     * <b>Author:</b> Youssef Elbarbary
388     * @param P univariate recursive GenPolynomial.
389     * @param S univariate recursive GenPolynomial.
390     * @return subResList(P,S).
391     */
392    public List<GenPolynomial<GenPolynomial<C>>> recursiveUnivariateSubResultantList(
393                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S) {
394        List<GenPolynomial<GenPolynomial<C>>> myList = new ArrayList<GenPolynomial<GenPolynomial<C>>>();
395        if (S == null || S.isZERO()) {
396            myList.add(S);
397            return myList;
398        }
399        if (P == null || P.isZERO()) {
400            myList.add(P);
401            return myList;
402        }
403        if (P.ring.nvar > 1) {
404            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
405        }
406        long e = P.degree(0);
407        long f = S.degree(0);
408        GenPolynomial<GenPolynomial<C>> q;
409        GenPolynomial<GenPolynomial<C>> r;
410        if (f > e) {
411            r = P;
412            q = S;
413            long g = f;
414            f = e;
415            e = g;
416        } else {
417            q = P;
418            r = S;
419        }
420        r = r.abs();
421        q = q.abs();
422        GenPolynomial<C> a = recursiveContent(r);
423        GenPolynomial<C> b = recursiveContent(q);
424        r = PolyUtil.<C> recursiveDivide(r, a);
425        q = PolyUtil.<C> recursiveDivide(q, b);
426        RingFactory<GenPolynomial<C>> cofac = P.ring.coFac;
427        GenPolynomial<C> g = cofac.getONE();
428        GenPolynomial<C> h = cofac.getONE();
429        GenPolynomial<GenPolynomial<C>> x;
430        GenPolynomial<C> t;
431        if (f == 0 && e == 0 && g.ring.nvar > 0) {
432            // if coeffs are multivariate (and non constant)
433            // otherwise it would be 1
434            t = resultant(a, b);
435            x = P.ring.getONE().multiply(t);
436            myList.add(x);
437            return myList;
438        }
439        t = a.power(e); // power(cofac, a, e);
440        t = t.multiply(b.power(f)); // power(cofac, b, f));
441        long s = 1;
442        GenPolynomial<C> z;
443        myList.add(P); // adding R0
444        myList.add(S); // adding R1
445        while (r.degree(0) > 0) {
446            long delta = q.degree(0) - r.degree(0);
447            if ((q.degree(0) % 2 != 0) && (r.degree(0) % 2 != 0)) {
448                s = -s;
449            }
450            x = PolyUtil.<C> recursiveDensePseudoRemainder(q, r);
451            q = r;
452            if (x.degree(0) >= 0) { // fixed: this was changed from > to >=
453                z = g.multiply(h.power(delta)); // power(P.ring.coFac, h, delta));
454                r = PolyUtil.<C> recursiveDivide(x, z);
455                myList.add(r);
456                g = q.leadingBaseCoefficient();
457                z = g.power(delta); // power(cofac, g, delta);
458                h = PolyUtil.<C> basePseudoDivide(z, h.power(delta - 1)); // power(cofac, h, delta - 1));
459            } else {
460                r = x;
461                myList.add(r);
462            }
463        }
464        z = r.leadingBaseCoefficient().power(q.degree(0)); // power(cofac, r.leadingBaseCoefficient(), q.degree(0));
465        h = PolyUtil.<C> basePseudoDivide(z, h.power(q.degree() - 1)); // power(cofac, h, q.degree(0) - 1));
466        z = h.multiply(t);
467        if (s < 0) {
468            z = z.negate();
469        }
470        x = P.ring.getONE().multiply(z);
471        myList.add(x);
472        // Printing the Subresultant List
473        //System.out.println("Liste von den SubResultanten(A - tD'):");
474        //for (int i = 0; i < myList.size(); i++) { // just for printing the list
475        //    System.out.println(myList.get(i));
476        //}
477        if (logger.isInfoEnabled()) {
478            System.out.println("subResCoeffs: " + myList);
479            //logger.info("subResCoeffs: " + myList);
480        }
481        return myList;
482    }
483
484
485    /**
486     * GenPolynomial base coefficient discriminant.
487     * @param P GenPolynomial.
488     * @return discriminant(P).
489     */
490    public GenPolynomial<C> baseDiscriminant(GenPolynomial<C> P) {
491        if (P == null) {
492            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
493        }
494        if (P.isZERO()) {
495            return P;
496        }
497        GenPolynomialRing<C> pfac = P.ring;
498        if (pfac.nvar > 1) {
499            throw new IllegalArgumentException(this.getClass().getName() + " P not univariate");
500        }
501        C a = P.leadingBaseCoefficient();
502        a = a.inverse();
503        GenPolynomial<C> Pp = PolyUtil.<C> baseDeriviative(P);
504        GenPolynomial<C> res = baseResultant(P, Pp);
505        GenPolynomial<C> disc = res.multiply(a);
506        long n = P.degree(0);
507        n = n * (n - 1);
508        n = n / 2;
509        if (n % 2L != 0L) {
510            disc = disc.negate();
511        }
512        return disc;
513    }
514
515}