Source code

001/*
002 * $Id: GreatestCommonDivisorSimple.java 4025 2012-07-23 16:41:43Z kredel $
003 */
004
005package edu.jas.ufd;
006
007
008import org.apache.log4j.Logger;
009
010import edu.jas.poly.GenPolynomial;
011import edu.jas.poly.PolyUtil;
012import edu.jas.structure.GcdRingElem;
013import edu.jas.structure.RingFactory;
014import edu.jas.structure.Power;
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 = Logger.getLogger(GreatestCommonDivisorSimple.class);
028
029
030    private 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            q = r;
169            if (field) {
170                r = PolyUtil.<C> monic(x);
171            } else {
172                r = x;
173            }
174        }
175        q = recursivePrimitivePart(q);
176        q = q.abs().multiply(c);
177        return q;
178    }
179
180
181    /**
182     * Univariate GenPolynomial resultant. 
183     * @param P univariate GenPolynomial.
184     * @param S univariate GenPolynomial.
185     * @return res(P,S).
186     */
187    @Override
188    public GenPolynomial<C> baseResultant(GenPolynomial<C> P, GenPolynomial<C> S) {
189        if (S == null || S.isZERO()) {
190            return S;
191        }
192        if (P == null || P.isZERO()) {
193            return P;
194        }
195        if (P.ring.nvar > 1 || P.ring.nvar == 0) {
196            throw new IllegalArgumentException("no univariate polynomial");
197        }
198        long e = P.degree(0);
199        long f = S.degree(0);
200        if (f == 0 && e == 0) {
201            return P.ring.getONE();
202        }
203        if ( e == 0 ) {
204            return Power.<GenPolynomial<C>> power(P.ring,P,f);
205        }
206        if ( f == 0 ) {
207            return Power.<GenPolynomial<C>> power(S.ring,S,e);
208        }
209        GenPolynomial<C> q;
210        GenPolynomial<C> r;
211        int s = 0; // sign is +, 1 for sign is -
212        if (e < f) {
213            r = P;
214            q = S;
215            long t = e;
216            e = f;
217            f = t;
218            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
219                s = 1;
220            }
221        } else {
222            q = P;
223            r = S;
224        }
225        RingFactory<C> cofac = P.ring.coFac; 
226        boolean field = cofac.isField();
227        C c = cofac.getONE();
228        GenPolynomial<C> x;
229        long g;
230        do {
231            if (field) {
232                x = q.remainder(r);
233            } else {
234                x = PolyUtil.<C>baseSparsePseudoRemainder(q,r);
235                //System.out.println("x_s = " + x + ", lbcf(r) = " + r.leadingBaseCoefficient());
236            }
237            if ( x.isZERO() ) {
238                return x;
239            }
240            //System.out.println("x = " + x);
241            e = q.degree(0);
242            f = r.degree(0);
243            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
244               s = 1 - s;
245            }
246            g = x.degree(0);
247            C c2 = r.leadingBaseCoefficient();
248            for (int i = 0; i < (e-g); i++ ) {
249                c = c.multiply(c2);
250            }
251            q = r; 
252            r = x;
253        } while (g != 0);
254        C c2 = r.leadingBaseCoefficient();
255        for (int i = 0; i < f; i++ ) {
256            c = c.multiply(c2);
257        }
258        if ( s == 1 ) {
259            c = c.negate();
260        }
261        x = P.ring.getONE().multiply(c);
262        return x;
263    }
264
265
266    /**
267     * Univariate GenPolynomial recursive resultant. 
268     * @param P univariate recursive GenPolynomial.
269     * @param S univariate recursive GenPolynomial.
270     * @return res(P,S).
271     */
272    @Override
273    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<C>> P,
274            GenPolynomial<GenPolynomial<C>> S) {
275        if (S == null || S.isZERO()) {
276            return S;
277        }
278        if (P == null || P.isZERO()) {
279            return P;
280        }
281        if (P.ring.nvar > 1 || P.ring.nvar == 0) {
282            throw new IllegalArgumentException("no recursive univariate polynomial");
283        }
284        long e = P.degree(0);
285        long f = S.degree(0);
286        if ( f == 0 && e == 0 ) {
287            // if coeffs are multivariate (and non constant)
288            // otherwise it would be 1
289            GenPolynomial<C> t = resultant(P.leadingBaseCoefficient(), S.leadingBaseCoefficient());
290            return P.ring.getONE().multiply(t);
291        }
292        if ( e == 0 ) {
293            return Power.<GenPolynomial<GenPolynomial<C>>> power(P.ring,P,f);
294        }
295        if ( f == 0 ) {
296            return Power.<GenPolynomial<GenPolynomial<C>>> power(S.ring,S,e);
297        }
298        GenPolynomial<GenPolynomial<C>> q;
299        GenPolynomial<GenPolynomial<C>> r;
300        int s = 0; // sign is +, 1 for sign is -
301        if (f > e) {
302            r = P;
303            q = S;
304            long g = f;
305            f = e;
306            e = g;
307            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
308                s = 1;
309            }
310        } else {
311            q = P;
312            r = S;
313        }
314        GenPolynomial<GenPolynomial<C>> x;
315        RingFactory<GenPolynomial<C>> cofac = P.ring.coFac; 
316        GenPolynomial<C> c = cofac.getONE();
317        long g;
318        do {
319            x = PolyUtil.<C>recursiveSparsePseudoRemainder(q,r);
320            //x = PolyUtil.<C>recursiveDensePseudoRemainder(q,r);
321            if ( x.isZERO() ) {
322                return x;
323            }
324            //no: x = recursivePrimitivePart(x);
325            //System.out.println("x = " + x);
326            e = q.degree(0);
327            f = r.degree(0);
328            if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f)
329               s = 1 - s;
330            }
331            g = x.degree(0);
332            GenPolynomial<C> c2 = r.leadingBaseCoefficient();
333            for (int i = 0; i < (e-g); i++ ) {
334                c = c.multiply(c2);
335            }
336            q = r; 
337            r = x;
338        } while (g != 0);
339        GenPolynomial<C> c2 = r.leadingBaseCoefficient();
340        for (int i = 0; i < f; i++ ) {
341             c = c.multiply(c2);
342        }
343        if ( s == 1 ) {
344            c = c.negate();
345        }
346        x = P.ring.getONE().multiply(c);
347        return x;
348    }
349
350}