001/*
002 * $Id$
003 */
004
005package edu.jas.ufd;
006
007
008import java.util.ArrayList;
009import java.util.List;
010
011import org.apache.logging.log4j.LogManager;
012import org.apache.logging.log4j.Logger;
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(forRemoval=true) 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(forRemoval=true) 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 common 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 = {}, f = {}", e, 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 common 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 = {}, f = {}", e, 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. <b>Author:</b> Youssef Elbarbary
387     * @param P univariate recursive GenPolynomial.
388     * @param S univariate recursive GenPolynomial.
389     * @return subResList(P,S).
390     */
391    public List<GenPolynomial<GenPolynomial<C>>> recursiveUnivariateSubResultantList(
392                    GenPolynomial<GenPolynomial<C>> P, GenPolynomial<GenPolynomial<C>> S) {
393        List<GenPolynomial<GenPolynomial<C>>> myList = new ArrayList<GenPolynomial<GenPolynomial<C>>>();
394        if (S == null || S.isZERO()) {
395            myList.add(S);
396            return myList;
397        }
398        if (P == null || P.isZERO()) {
399            myList.add(P);
400            return myList;
401        }
402        if (P.ring.nvar > 1) {
403            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
404        }
405        long e = P.degree(0);
406        long f = S.degree(0);
407        GenPolynomial<GenPolynomial<C>> q;
408        GenPolynomial<GenPolynomial<C>> r;
409        if (f > e) {
410            r = P;
411            q = S;
412            long g = f;
413            f = e;
414            e = g;
415        } else {
416            q = P;
417            r = S;
418        }
419        r = r.abs();
420        q = q.abs();
421        GenPolynomial<C> a = recursiveContent(r);
422        GenPolynomial<C> b = recursiveContent(q);
423        r = PolyUtil.<C> recursiveDivide(r, a);
424        q = PolyUtil.<C> recursiveDivide(q, b);
425        RingFactory<GenPolynomial<C>> cofac = P.ring.coFac;
426        GenPolynomial<C> g = cofac.getONE();
427        GenPolynomial<C> h = cofac.getONE();
428        GenPolynomial<GenPolynomial<C>> x;
429        GenPolynomial<C> t;
430        if (f == 0 && e == 0 && g.ring.nvar > 0) {
431            // if coeffs are multivariate (and non constant)
432            // otherwise it would be 1
433            t = resultant(a, b);
434            x = P.ring.getONE().multiply(t);
435            myList.add(x);
436            return myList;
437        }
438        t = a.power(e); // power(cofac, a, e);
439        t = t.multiply(b.power(f)); // power(cofac, b, f));
440        long s = 1;
441        GenPolynomial<C> z;
442        myList.add(P); // adding R0
443        myList.add(S); // adding R1
444        while (r.degree(0) > 0) {
445            long delta = q.degree(0) - r.degree(0);
446            if ((q.degree(0) % 2 != 0) && (r.degree(0) % 2 != 0)) {
447                s = -s;
448            }
449            x = PolyUtil.<C> recursiveDensePseudoRemainder(q, r);
450            q = r;
451            if (x.degree(0) >= 0) { // fixed: this was changed from > to >=
452                z = g.multiply(h.power(delta)); // power(P.ring.coFac, h, delta));
453                r = PolyUtil.<C> recursiveDivide(x, z);
454                myList.add(r);
455                g = q.leadingBaseCoefficient();
456                z = g.power(delta); // power(cofac, g, delta);
457                h = PolyUtil.<C> basePseudoDivide(z, h.power(delta - 1)); // power(cofac, h, delta - 1));
458            } else {
459                r = x;
460                myList.add(r);
461            }
462        }
463        z = r.leadingBaseCoefficient().power(q.degree(0)); // power(cofac, r.leadingBaseCoefficient(), q.degree(0));
464        h = PolyUtil.<C> basePseudoDivide(z, h.power(q.degree() - 1)); // power(cofac, h, q.degree(0) - 1));
465        z = h.multiply(t);
466        if (s < 0) {
467            z = z.negate();
468        }
469        x = P.ring.getONE().multiply(z);
470        myList.add(x);
471        // Printing the Subresultant List
472        //System.out.println("Liste von den SubResultanten(A - tD'):");
473        //for (int i = 0; i < myList.size(); i++) { // just for printing the list
474        //    System.out.println(myList.get(i));
475        //}
476        if (logger.isInfoEnabled()) {
477            System.out.println("subResCoeffs: " + myList);
478            //logger.info("subResCoeffs: {}", myList);
479        }
480        return myList;
481    }
482
483
484    /**
485     * GenPolynomial base coefficient discriminant.
486     * @param P GenPolynomial.
487     * @return discriminant(P).
488     */
489    public GenPolynomial<C> baseDiscriminant(GenPolynomial<C> P) {
490        if (P == null) {
491            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
492        }
493        if (P.isZERO()) {
494            return P;
495        }
496        GenPolynomialRing<C> pfac = P.ring;
497        if (pfac.nvar > 1) {
498            throw new IllegalArgumentException(this.getClass().getName() + " P not univariate");
499        }
500        C a = P.leadingBaseCoefficient();
501        a = a.inverse();
502        GenPolynomial<C> Pp = PolyUtil.<C> baseDerivative(P);
503        GenPolynomial<C> res = baseResultant(P, Pp);
504        GenPolynomial<C> disc = res.multiply(a);
505        long n = P.degree(0);
506        n = n * (n - 1);
507        n = n / 2;
508        if (n % 2L != 0L) {
509            disc = disc.negate();
510        }
511        return disc;
512    }
513
514}