001/*
002 * $Id: GreatestCommonDivisorSimple.java 5105 2015-02-07 14:35:17Z kredel $
003 */
004
005package edu.jas.fd;
006
007
008import org.apache.log4j.Logger;
009
010import edu.jas.poly.GenPolynomial;
011import edu.jas.poly.GenSolvablePolynomial;
012import edu.jas.poly.PolyUtil;
013import edu.jas.structure.GcdRingElem;
014import edu.jas.structure.RingFactory;
015
016
017/**
018 * (Non-unique) factorization domain greatest common divisor common algorithms
019 * with monic polynomial remainder sequence. If C is a field, then the monic PRS
020 * (on coefficients) is computed otherwise no simplifications in the reduction
021 * are made.
022 * @param <C> coefficient type
023 * @author Heinz Kredel
024 */
025
026public class GreatestCommonDivisorSimple<C extends GcdRingElem<C>> extends GreatestCommonDivisorAbstract<C> {
027
028
029    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorSimple.class);
030
031
032    private final boolean debug = true; //logger.isDebugEnabled();
033
034
035    /**
036     * Constructor.
037     * @param cf coefficient ring.
038     */
039    public GreatestCommonDivisorSimple(RingFactory<C> cf) {
040        super(cf);
041    }
042
043
044    /**
045     * Univariate GenSolvablePolynomial greatest common divisor. Uses
046     * pseudoRemainder for remainder.
047     * @param P univariate GenSolvablePolynomial.
048     * @param S univariate GenSolvablePolynomial.
049     * @return gcd(P,S) with P = P'*gcd(P,S) and S = S'*gcd(P,S).
050     */
051    @Override
052    public GenSolvablePolynomial<C> leftBaseGcd(GenSolvablePolynomial<C> P, GenSolvablePolynomial<C> S) {
053        if (S == null || S.isZERO()) {
054            return P;
055        }
056        if (P == null || P.isZERO()) {
057            return S;
058        }
059        if (P.ring.nvar > 1) {
060            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
061        }
062        boolean field = P.ring.coFac.isField();
063        long e = P.degree(0);
064        long f = S.degree(0);
065        GenSolvablePolynomial<C> q;
066        GenSolvablePolynomial<C> r;
067        if (f > e) {
068            r = P;
069            q = S;
070            long g = f;
071            f = e;
072            e = g;
073        } else {
074            q = P;
075            r = S;
076        }
077        if (debug) {
078            logger.debug("degrees: e = " + e + ", f = " + f);
079        }
080        C c;
081        if (field) {
082            r = r.monic();
083            q = q.monic();
084            c = P.ring.getONECoefficient();
085        } else {
086            r = (GenSolvablePolynomial<C>) r.abs();
087            q = (GenSolvablePolynomial<C>) q.abs();
088            C a = rightBaseContent(r);
089            C b = rightBaseContent(q);
090            r = divide(r, a); // indirection
091            q = divide(q, b); // indirection
092            c = gcd(a, b); // indirection
093        }
094        //System.out.println("baseCont: gcd(cont) = " + b);
095        if (r.isONE()) {
096            return r.multiply(c);
097        }
098        if (q.isONE()) {
099            return q.multiply(c);
100        }
101        GenSolvablePolynomial<C> x;
102        //System.out.println("baseGCD: q = " + q);
103        //System.out.println("baseGCD: r = " + r);
104        while (!r.isZERO()) {
105            x = FDUtil.<C> leftBaseSparsePseudoRemainder(q, r);
106            q = r;
107            if (field) {
108                r = x.monic();
109            } else {
110                r = x;
111            }
112            //System.out.println("baseGCD: q = " + q);
113            //System.out.println("baseGCD: r = " + r);
114        }
115        q = leftBasePrimitivePart(q);
116        return (GenSolvablePolynomial<C>) (q.multiply(c)).abs();
117    }
118
119
120    /**
121     * Univariate GenSolvablePolynomial right greatest common divisor. Uses
122     * pseudoRemainder for remainder.
123     * @param P univariate GenSolvablePolynomial.
124     * @param S univariate GenSolvablePolynomial.
125     * @return gcd(P,S) with P = gcd(P,S)*P' and S = gcd(P,S)*S'.
126     */
127    @Override
128    public GenSolvablePolynomial<C> rightBaseGcd(GenSolvablePolynomial<C> P, GenSolvablePolynomial<C> S) {
129        if (S == null || S.isZERO()) {
130            return P;
131        }
132        if (P == null || P.isZERO()) {
133            return S;
134        }
135        if (P.ring.nvar > 1) {
136            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
137        }
138        boolean field = P.ring.coFac.isField();
139        long e = P.degree(0);
140        long f = S.degree(0);
141        GenSolvablePolynomial<C> q;
142        GenSolvablePolynomial<C> r;
143        if (f > e) {
144            r = P;
145            q = S;
146            long g = f;
147            f = e;
148            e = g;
149        } else {
150            q = P;
151            r = S;
152        }
153        if (debug) {
154            logger.debug("degrees: e = " + e + ", f = " + f);
155        }
156        C c;
157        if (field) {
158            r = r.monic();
159            q = q.monic();
160            c = P.ring.getONECoefficient();
161        } else {
162            r = (GenSolvablePolynomial<C>) r.abs();
163            q = (GenSolvablePolynomial<C>) q.abs();
164            C a = leftBaseContent(r);
165            C b = leftBaseContent(q);
166            r = divide(r, a); // indirection
167            q = divide(q, b); // indirection
168            c = gcd(a, b); // indirection
169        }
170        //System.out.println("baseCont: gcd(cont) = " + b);
171        if (r.isONE()) {
172            return r.multiply(c);
173        }
174        if (q.isONE()) {
175            return q.multiply(c);
176        }
177        GenSolvablePolynomial<C> x;
178        //System.out.println("baseGCD: q = " + q);
179        //System.out.println("baseGCD: r = " + r);
180        while (!r.isZERO()) {
181            x = FDUtil.<C> rightBaseSparsePseudoRemainder(q, r);
182            q = r;
183            if (field) {
184                r = x.monic();
185            } else {
186                r = x;
187            }
188            //System.out.println("baseGCD: q = " + q);
189            //System.out.println("baseGCD: r = " + r);
190        }
191        q = rightBasePrimitivePart(q);
192        return (GenSolvablePolynomial<C>) (q.multiplyLeft(c)).abs();
193    }
194
195
196    /**
197     * Univariate GenSolvablePolynomial left recursive greatest common divisor.
198     * Uses pseudoRemainder for remainder.
199     * @param P univariate recursive GenSolvablePolynomial.
200     * @param S univariate recursive GenSolvablePolynomial.
201     * @return gcd(P,S) with P = P'*gcd(P,S)*p and S = S'*gcd(P,S)*s, where
202     *         deg_main(p) = deg_main(s) == 0.
203     */
204    @Override
205    public GenSolvablePolynomial<GenPolynomial<C>> leftRecursiveUnivariateGcd(
206                    GenSolvablePolynomial<GenPolynomial<C>> P, GenSolvablePolynomial<GenPolynomial<C>> S) {
207        if (S == null || S.isZERO()) {
208            return P;
209        }
210        if (P == null || P.isZERO()) {
211            return S;
212        }
213        if (P.ring.nvar > 1) {
214            throw new IllegalArgumentException("no univariate polynomial");
215        }
216        boolean field = P.leadingBaseCoefficient().ring.coFac.isField();
217        long e = P.degree(0);
218        long f = S.degree(0);
219        GenSolvablePolynomial<GenPolynomial<C>> q, r, x, qs, rs, qp, rp;
220        if (f > e) {
221            r = P;
222            q = S;
223            long g = f;
224            f = e;
225            e = g;
226        } else if (f < e) {
227            q = P;
228            r = S;
229        } else { // f == e
230            if (P.leadingBaseCoefficient().degree() > S.leadingBaseCoefficient().degree()) {
231                q = P;
232                r = S;
233            } else {
234                r = P;
235                q = S;
236            }
237        }
238        if (debug) {
239            logger.debug("degrees: e = " + e + ", f = " + f);
240        }
241        if (field) {
242            r = PolyUtil.<C> monic(r);
243            q = PolyUtil.<C> monic(q);
244        } else {
245            r = (GenSolvablePolynomial<GenPolynomial<C>>) r.abs();
246            q = (GenSolvablePolynomial<GenPolynomial<C>>) q.abs();
247        }
248        GenSolvablePolynomial<C> a = rightRecursiveContent(r);
249        rs = FDUtil.<C> recursiveDivideRightEval(r, a);
250        if (debug) {
251            logger.info("recCont a = " + a + ", r = " + r);
252            logger.info("recCont r/a = " + rs + ", r%a = " + r.subtract(rs.multiply(a)));
253            if (!r.equals(rs.multiply(a))) {
254                if (!rs.multiplyLeft(a).equals(r)) {
255                    System.out.println("recGcd, r         = " + r);
256                    System.out.println("recGcd, cont(r)   = " + a);
257                    System.out.println("recGcd, pp(r)     = " + rs);
258                    System.out.println("recGcd, pp(r)c(r) = " + rs.multiply(a));
259                    System.out.println("recGcd, c(r)pp(r) = " + rs.multiplyLeft(a));
260                    throw new RuntimeException("recGcd, pp: not divisible");
261                }
262            }
263        }
264        r = rs;
265        GenSolvablePolynomial<C> b = rightRecursiveContent(q);
266        qs = FDUtil.<C> recursiveDivideRightEval(q, b);
267        if (debug) {
268            logger.info("recCont b = " + b + ", q = " + q);
269            logger.info("recCont q/b = " + qs + ", q%b = " + q.subtract(qs.multiply(b)));
270            if (!q.equals(qs.multiply(b))) {
271                if (!qs.multiplyLeft(b).equals(q)) {
272                    System.out.println("recGcd, q         = " + q);
273                    System.out.println("recGcd, cont(q)   = " + b);
274                    System.out.println("recGcd, pp(q)     = " + qs);
275                    System.out.println("recGcd, pp(q)c(q) = " + qs.multiply(b));
276                    System.out.println("recGcd, c(q)pp(q) = " + qs.multiplyLeft(b));
277                    throw new RuntimeException("recGcd, pp: not divisible");
278                }
279            }
280        }
281        q = qs;
282        logger.info("Gcd(content).ring = " + a.ring.toScript() + ", a = " + a + ", b = " + b);
283        //no: GenSolvablePolynomial<C> c = leftGcd(a, b); // go to recursion
284        GenSolvablePolynomial<C> c = rightGcd(a, b); // go to recursion
285        logger.info("Gcd(contents) c = " + c + ", a = " + a + ", b = " + b);
286        if (r.isONE()) {
287            return r.multiply(c);
288        }
289        if (q.isONE()) {
290            return q.multiply(c);
291        }
292        rs = r;
293        qs = q;
294        if (debug) {
295            logger.info("r.ring = " + r.ring.toScript());
296            logger.info("gcd-loop, start: q = " + q + ", r = " + r);
297        }
298        while (!r.isZERO()) {
299            x = FDUtil.<C> recursiveSparsePseudoRemainder(q, r);
300            q = r;
301            if (field) {
302                r = PolyUtil.<C> monic(x);
303            } else {
304                r = x;
305            }
306            if (debug) {
307                logger.info("gcd-loop, rem: q = " + q + ", r = " + r);
308            }
309        }
310        if (debug) {
311            logger.info("gcd(div) = " + q + ", rs = " + rs + ", qs = " + qs);
312            rp = FDUtil.<C> recursiveSparsePseudoRemainder(rs, q);
313            qp = FDUtil.<C> recursiveSparsePseudoRemainder(qs, q);
314            if (!qp.isZERO() || !rp.isZERO()) {
315                logger.info("gcd(div): rem(r,g) = " + rp + ", rem(q,g) = " + qp);
316                rp = FDUtil.<C> recursivePseudoQuotient(rs, q);
317                qp = FDUtil.<C> recursivePseudoQuotient(qs, q);
318                logger.info("gcd(div): r/g = " + rp + ", q/g = " + qp);
319                throw new RuntimeException("recGcd, div: not divisible");
320            }
321        }
322        qp = q;
323        if (false) { // not checkable
324            System.out.println("recGcd, qs = " + qs + ", P = " + P);
325            System.out.println("recGcd, rs = " + rs + ", S = " + S);
326            qs = FDUtil.<C> recursiveSparsePseudoRemainder(P, qp); // P
327            rs = FDUtil.<C> recursiveSparsePseudoRemainder(S, qp); // S
328            if (!qs.isZERO() || !rs.isZERO()) {
329                System.out.println("recGcd, P  = " + P);
330                System.out.println("recGcd, S  = " + S);
331                System.out.println("recGcd, qp = " + qp);
332                System.out.println("recGcd, qs = " + qs);
333                System.out.println("recGcd, rs = " + rs);
334                rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, qp);
335                System.out.println("recGcd, rs = " + rs);
336                if (!rs.isZERO()) {
337                    throw new RuntimeException("recGcd, pre pp: not divisible");
338                }
339            }
340            logger.info("recGcd(P,S) okay: qp = " + qp);
341        }
342        q = leftRecursivePrimitivePart(q);
343        if (!qp.equals(q)) {
344            logger.info("gcd(pp) = " + q + ", qp = " + qp);
345        }
346        q = (GenSolvablePolynomial<GenPolynomial<C>>) q.multiply(c).abs();
347        if (false) { // not checkable
348            logger.info("gcd(pp) = " + q + ", c = " + c);
349            qs = FDUtil.<C> recursiveSparsePseudoRemainder(P, q);
350            rs = FDUtil.<C> recursiveSparsePseudoRemainder(S, q);
351            if (!qs.isZERO() || !rs.isZERO()) {
352                System.out.println("recGcd-pp, P  = " + P);
353                System.out.println("recGcd-pp, S  = " + S);
354                System.out.println("recGcd-pp, q  = " + q);
355                System.out.println("recGcd-pp, qs = " + qs);
356                System.out.println("recGcd-pp, rs = " + rs);
357                rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, q);
358                System.out.println("recGcd-pp, rs = " + rs);
359                if (!rs.isZERO()) {
360                    throw new RuntimeException("recGcd: not divisible");
361                }
362            }
363            logger.info("recGcd(P,S) okay: q = " + q);
364        }
365        return q;
366    }
367
368
369    /**
370     * Univariate GenSolvablePolynomial right recursive greatest common divisor.
371     * Uses pseudoRemainder for remainder.
372     * @param P univariate recursive GenSolvablePolynomial.
373     * @param S univariate recursive GenSolvablePolynomial.
374     * @return gcd(P,S) with P = p*gcd(P,S)*P' and S = s*gcd(P,S)*S', where
375     *         deg_main(p) = deg_main(s) == 0.
376     */
377    @Override
378    public GenSolvablePolynomial<GenPolynomial<C>> rightRecursiveUnivariateGcd(
379                    GenSolvablePolynomial<GenPolynomial<C>> P, GenSolvablePolynomial<GenPolynomial<C>> S) {
380        if (S == null || S.isZERO()) {
381            return P;
382        }
383        if (P == null || P.isZERO()) {
384            return S;
385        }
386        if (P.ring.nvar > 1) {
387            throw new IllegalArgumentException("no univariate polynomial");
388        }
389        boolean field = P.leadingBaseCoefficient().ring.coFac.isField();
390        long e = P.degree(0);
391        long f = S.degree(0);
392        GenSolvablePolynomial<GenPolynomial<C>> q, r, x, qs, rs, qp, rp;
393        if (f > e) {
394            r = P;
395            q = S;
396            long g = f;
397            f = e;
398            e = g;
399        } else if (f < e) {
400            q = P;
401            r = S;
402        } else { // f == e
403            if (P.leadingBaseCoefficient().degree() > S.leadingBaseCoefficient().degree()) {
404                q = P;
405                r = S;
406            } else {
407                r = P;
408                q = S;
409            }
410        }
411        if (debug) {
412            logger.debug("RI-degrees: e = " + e + ", f = " + f);
413        }
414        if (field) {
415            r = PolyUtil.<C> monic(r);
416            q = PolyUtil.<C> monic(q);
417        } else {
418            r = (GenSolvablePolynomial<GenPolynomial<C>>) r.abs();
419            q = (GenSolvablePolynomial<GenPolynomial<C>>) q.abs();
420        }
421        GenSolvablePolynomial<C> a = leftRecursiveContent(r);
422        rs = FDUtil.<C> recursiveDivide(r, a);
423        if (debug) {
424            logger.info("RI-recCont a = " + a + ", r = " + r);
425            logger.info("RI-recCont r/a = " + r + ", r%a = " + r.subtract(rs.multiplyLeft(a)));
426            if (!r.equals(rs.multiplyLeft(a))) {
427                System.out.println("RI-recGcd, r         = " + r);
428                System.out.println("RI-recGcd, cont(r)   = " + a);
429                System.out.println("RI-recGcd, pp(r)     = " + rs);
430                System.out.println("RI-recGcd, pp(r)c(r) = " + rs.multiplyLeft(a));
431                throw new RuntimeException("RI-recGcd, pp: not divisible");
432            }
433        }
434        r = rs;
435        GenSolvablePolynomial<C> b = leftRecursiveContent(q);
436        qs = FDUtil.<C> recursiveDivide(q, b);
437        if (debug) {
438            logger.info("RI-recCont b = " + b + ", q = " + q);
439            logger.info("RI-recCont q/b = " + qs + ", q%b = " + q.subtract(qs.multiplyLeft(b)));
440            if (!q.equals(qs.multiplyLeft(b))) {
441                System.out.println("RI-recGcd, q         = " + q);
442                System.out.println("RI-recGcd, cont(q)   = " + b);
443                System.out.println("RI-recGcd, pp(q)     = " + qs);
444                System.out.println("RI-recGcd, pp(q)c(q) = " + qs.multiplyLeft(b));
445                throw new RuntimeException("RI-recGcd, pp: not divisible");
446            }
447        }
448        q = qs;
449        //no: GenSolvablePolynomial<C> c = rightGcd(a, b); // go to recursion
450        GenSolvablePolynomial<C> c = leftGcd(a, b); // go to recursion
451        logger.info("RI-Gcd(contents) c = " + c + ", a = " + a + ", b = " + b);
452        if (r.isONE()) {
453            return r.multiplyLeft(c);
454        }
455        if (q.isONE()) {
456            return q.multiplyLeft(c);
457        }
458        rs = r;
459        qs = q;
460        if (debug) {
461            logger.info("RI-r.ring = " + r.ring.toScript());
462            logger.info("RI-gcd-loop, start: q = " + q + ", r = " + r);
463        }
464        while (!r.isZERO()) {
465            x = FDUtil.<C> recursiveRightSparsePseudoRemainder(q, r);
466            q = r;
467            if (field) {
468                r = PolyUtil.<C> monic(x);
469            } else {
470                r = x;
471            }
472            if (debug) {
473                logger.info("RI-gcd-loop, rem: q = " + q + ", r = " + r);
474            }
475        }
476        if (debug) {
477            logger.info("RI-gcd(div) = " + q + ", rs = " + rs + ", qs = " + qs);
478            rp = FDUtil.<C> recursiveRightSparsePseudoRemainder(rs, q);
479            qp = FDUtil.<C> recursiveRightSparsePseudoRemainder(qs, q);
480            if (!qp.isZERO() || !rp.isZERO()) {
481                logger.info("RI-gcd(div): rem(r,g) = " + rp + ", rem(q,g) = " + qp);
482                rp = FDUtil.<C> recursivePseudoQuotient(rs, q);
483                qp = FDUtil.<C> recursivePseudoQuotient(qs, q);
484                logger.info("RI-gcd(div): r/g = " + rp + ", q/g = " + qp);
485                //logger.info("gcd(div): rp*g = " + rp.multiply(q) + ", qp*g = " + qp.multiply(q));
486                throw new RuntimeException("recGcd, div: not divisible");
487            }
488        }
489        qp = q;
490        if (false) { // not checkable
491            System.out.println("RI-recGcd, qs = " + qs + ", P = " + P);
492            System.out.println("RI-recGcd, rs = " + rs + ", S = " + S);
493            qs = FDUtil.<C> recursiveRightSparsePseudoRemainder(P, qp); // P
494            rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, qp); // S
495            if (!qs.isZERO() || !rs.isZERO()) {
496                System.out.println("RI-recGcd, P  = " + P);
497                System.out.println("RI-recGcd, S  = " + S);
498                System.out.println("RI-recGcd, qp = " + qp);
499                System.out.println("RI-recGcd, qs = " + qs);
500                System.out.println("RI-recGcd, rs = " + rs);
501                throw new RuntimeException("RI-recGcd, pre pp: not divisible");
502            }
503            logger.info("RI-recGcd(P,S) pre pp okay: qp = " + qp);
504        }
505        q = rightRecursivePrimitivePart(q);
506        if (!qp.equals(q)) {
507            logger.info("RI-gcd(pp) = " + q + ", qp = " + qp); // + ", ring = " + P.ring.toScript());
508        }
509        q = (GenSolvablePolynomial<GenPolynomial<C>>) q.multiplyLeft(c).abs();
510        if (false) { // not checkable
511            qs = FDUtil.<C> recursiveRightSparsePseudoRemainder(P, q);
512            rs = FDUtil.<C> recursiveRightSparsePseudoRemainder(S, q);
513            if (!qs.isZERO() || !rs.isZERO()) {
514                System.out.println("RI-recGcd-pp, P  = " + P);
515                System.out.println("RI-recGcd-pp, S  = " + S);
516                System.out.println("RI-recGcd-pp, q  = " + q);
517                System.out.println("RI-recGcd-pp, qs = " + qs);
518                System.out.println("RI-recGcd-pp, rs = " + rs);
519                throw new RuntimeException("RI-recGcd: not divisible");
520            }
521            logger.info("RI-recGcd(P,S) post pp okay: q = " + q);
522        }
523        return q;
524    }
525
526}