Source code

001/*
002 * $Id: GreatestCommonDivisorSimple.java 5871 2018-07-20 15:58:45Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import org.apache.logging.log4j.Logger;
009import org.apache.logging.log4j.LogManager; 
010
011import edu.jas.poly.GenPolynomial;
012import edu.jas.poly.PolyUtil;
013import edu.jas.structure.GcdRingElem;
014import edu.jas.structure.RingFactory;
015
016
017/**
018 * Greatest common divisor algorithms with monic polynomial remainder sequence.
019 * If C is a field, then the monic PRS (on coefficients) is computed otherwise
020 * no simplifications in the reduction are made.
021 * @author Heinz Kredel
022 */
023
024public class GreatestCommonDivisorSimple<C extends GcdRingElem<C>> extends GreatestCommonDivisorAbstract<C> {
025
026
027    private static final Logger logger = LogManager.getLogger(GreatestCommonDivisorSimple.class);
028
029
030    private static final boolean debug = logger.isDebugEnabled();
031
032
033    /**
034     * Univariate GenPolynomial greatest comon divisor. Uses pseudoRemainder for
035     * remainder.
036     * @param P univariate GenPolynomial.
037     * @param S univariate GenPolynomial.
038     * @return gcd(P,S).
039     */
040    @Override
041    public GenPolynomial<C> baseGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
042        if (S == null || S.isZERO()) {
043            return P;
044        }
045        if (P == null || P.isZERO()) {
046            return S;
047        }
048        if (P.ring.nvar > 1) {
049            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
050        }
051        boolean field = P.ring.coFac.isField();
052        long e = P.degree(0);
053        long f = S.degree(0);
054        GenPolynomial<C> q;
055        GenPolynomial<C> r;
056        if (f > e) {
057            r = P;
058            q = S;
059            long g = f;
060            f = e;
061            e = g;
062        } else {
063            q = P;
064            r = S;
065        }
066        if (debug) {
067            logger.debug("degrees: e = " + e + ", f = " + f);
068        }
069        C c;
070        if (field) {
071            r = r.monic();
072            q = q.monic();
073            c = P.ring.getONECoefficient();
074        } else {
075            r = r.abs();
076            q = q.abs();
077            C a = baseContent(r);
078            C b = baseContent(q);
079            c = gcd(a, b); // indirection
080            r = divide(r, a); // indirection
081            q = divide(q, b); // indirection
082        }
083        if (r.isONE()) {
084            return r.multiply(c);
085        }
086        if (q.isONE()) {
087            return q.multiply(c);
088        }
089        GenPolynomial<C> x;
090        //System.out.println("q = " + q);
091        //System.out.println("r = " + r);
092        while (!r.isZERO()) {
093            x = PolyUtil.<C> baseSparsePseudoRemainder(q, r);
094            q = r;
095            if (field) {
096                r = x.monic();
097            } else {
098                r = x;
099            }
100            //System.out.println("q = " + q);
101            //System.out.println("r = " + r);
102        }
103        q = basePrimitivePart(q);
104        return (q.multiply(c)).abs();
105    }
106
107
108    /**
109     * Univariate GenPolynomial recursive greatest comon divisor. Uses
110     * pseudoRemainder for remainder.
111     * @param P univariate recursive GenPolynomial.
112     * @param S univariate recursive GenPolynomial.
113     * @return gcd(P,S).
114     */
115    @Override
116    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateGcd(GenPolynomial<GenPolynomial<C>> P,
117                    GenPolynomial<GenPolynomial<C>> S) {
118        if (S == null || S.isZERO()) {
119            return P;
120        }
121        if (P == null || P.isZERO()) {
122            return S;
123        }
124        if (P.ring.nvar > 1) {
125            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
126        }
127        boolean field = P.leadingBaseCoefficient().ring.coFac.isField();
128        long e = P.degree(0);
129        long f = S.degree(0);
130        GenPolynomial<GenPolynomial<C>> q;
131        GenPolynomial<GenPolynomial<C>> r;
132        if (f > e) {
133            r = P;
134            q = S;
135            long g = f;
136            f = e;
137            e = g;
138        } else {
139            q = P;
140            r = S;
141        }
142        if (debug) {
143            logger.debug("degrees: e = " + e + ", f = " + f);
144        }
145        if (field) {
146            r = PolyUtil.<C> monic(r);
147            q = PolyUtil.<C> monic(q);
148        } else {
149            r = r.abs();
150            q = q.abs();
151        }
152        GenPolynomial<C> a = recursiveContent(r);
153        GenPolynomial<C> b = recursiveContent(q);
154
155        GenPolynomial<C> c = gcd(a, b); // go to recursion
156        //System.out.println("rgcd c = " + c);
157        r = PolyUtil.<C> recursiveDivide(r, a);
158        q = PolyUtil.<C> recursiveDivide(q, b);
159        if (r.isONE()) {
160            return r.multiply(c);
161        }
162        if (q.isONE()) {
163            return q.multiply(c);
164        }
165        GenPolynomial<GenPolynomial<C>> x;
166        while (!r.isZERO()) {
167            x = PolyUtil.<C> recursivePseudoRemainder(q, r);
168            if (logger.isDebugEnabled()) {
169                logger.info("recursivePseudoRemainder.bits = " + x.bitLength());
170            }
171            q = r;
172            if (field) {
173                r = PolyUtil.<C> monic(x);
174            } else {
175                r = x;
176            }
177        }
178        q = recursivePrimitivePart(q);
179        q = q.abs().multiply(c);
180        return q;
181    }
182
183
184    /**
185     * Univariate GenPolynomial resultant.
186     * @param P univariate GenPolynomial.
187     * @param S univariate GenPolynomial.
188     * @return res(P,S).
189     */
190    @Override
191    public GenPolynomial<C> baseResultant(GenPolynomial<C> P, GenPolynomial<C> S) {
192        if (S == null || S.isZERO()) {
193            return S;
194        }
195        if (P == null || P.isZERO()) {
196            return P;
197        }
198        if (P.ring.nvar > 1 || P.ring.nvar == 0) {
199            throw new IllegalArgumentException("no univariate polynomial");
200        }
201        long e = P.degree(0);
202        long f = S.degree(0);
203        if (f == 0 && e == 0) {
204            return P.ring.getONE();
205        }
206        if (e == 0) {
207            return P.power(f);
208        }
209        if (f == 0) {
210            return S.power(e);
211        }
212        GenPolynomial<C> q;
213        GenPolynomial<C> r;
214        int s = 0; // sign is +, 1 for sign is -
215        if (e < f) {
216            r = P;
217            q = S;
218            long t = e;
219            e = f;
220            f = t;
221            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
222                s = 1;
223            }
224        } else {
225            q = P;
226            r = S;
227        }
228        RingFactory<C> cofac = P.ring.coFac;
229        boolean field = cofac.isField();
230        C c = cofac.getONE();
231        GenPolynomial<C> x;
232        long g;
233        do {
234            if (field) {
235                x = q.remainder(r);
236            } else {
237                x = PolyUtil.<C> baseSparsePseudoRemainder(q, r);
238                //System.out.println("x_s = " + x + ", lbcf(r) = " + r.leadingBaseCoefficient());
239            }
240            if (x.isZERO()) {
241                return x;
242            }
243            //System.out.println("x = " + x);
244            e = q.degree(0);
245            f = r.degree(0);
246            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
247                s = 1 - s;
248            }
249            g = x.degree(0);
250            C c2 = r.leadingBaseCoefficient();
251            for (int i = 0; i < (e - g); i++) {
252                c = c.multiply(c2);
253            }
254            q = r;
255            r = x;
256        } while (g != 0);
257        C c2 = r.leadingBaseCoefficient();
258        for (int i = 0; i < f; i++) {
259            c = c.multiply(c2);
260        }
261        if (s == 1) {
262            c = c.negate();
263        }
264        x = P.ring.getONE().multiply(c);
265        return x;
266    }
267
268
269    /**
270     * Univariate GenPolynomial recursive resultant.
271     * @param P univariate recursive GenPolynomial.
272     * @param S univariate recursive GenPolynomial.
273     * @return res(P,S).
274     */
275    @Override
276    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<C>> P,
277                    GenPolynomial<GenPolynomial<C>> S) {
278        if (S == null || S.isZERO()) {
279            return S;
280        }
281        if (P == null || P.isZERO()) {
282            return P;
283        }
284        if (P.ring.nvar > 1 || P.ring.nvar == 0) {
285            throw new IllegalArgumentException("no recursive univariate polynomial");
286        }
287        long e = P.degree(0);
288        long f = S.degree(0);
289        if (f == 0 && e == 0) {
290            // if coeffs are multivariate (and non constant)
291            // otherwise it would be 1
292            GenPolynomial<C> t = resultant(P.leadingBaseCoefficient(), S.leadingBaseCoefficient());
293            return P.ring.getONE().multiply(t);
294        }
295        if (e == 0) {
296            return P.power(f);
297        }
298        if (f == 0) {
299            return S.power(e);
300        }
301        GenPolynomial<GenPolynomial<C>> q;
302        GenPolynomial<GenPolynomial<C>> r;
303        int s = 0; // sign is +, 1 for sign is -
304        if (f > e) {
305            r = P;
306            q = S;
307            long g = f;
308            f = e;
309            e = g;
310            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
311                s = 1;
312            }
313        } else {
314            q = P;
315            r = S;
316        }
317        GenPolynomial<GenPolynomial<C>> x;
318        RingFactory<GenPolynomial<C>> cofac = P.ring.coFac;
319        GenPolynomial<C> c = cofac.getONE();
320        long g;
321        do {
322            x = PolyUtil.<C> recursiveSparsePseudoRemainder(q, r);
323            //x = PolyUtil.<C>recursiveDensePseudoRemainder(q,r);
324            if (x.isZERO()) {
325                return x;
326            }
327            //no: x = recursivePrimitivePart(x);
328            //System.out.println("x = " + x);
329            e = q.degree(0);
330            f = r.degree(0);
331            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
332                s = 1 - s;
333            }
334            g = x.degree(0);
335            GenPolynomial<C> c2 = r.leadingBaseCoefficient();
336            for (int i = 0; i < (e - g); i++) {
337                c = c.multiply(c2);
338            }
339            q = r;
340            r = x;
341        } while (g != 0);
342        GenPolynomial<C> c2 = r.leadingBaseCoefficient();
343        for (int i = 0; i < f; i++) {
344            c = c.multiply(c2);
345        }
346        if (s == 1) {
347            c = c.negate();
348        }
349        x = P.ring.getONE().multiply(c);
350        return x;
351    }
352
353}