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