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