001/*
002 * $Id: GreatestCommonDivisorSubres.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.ExpVector;
011import edu.jas.poly.GenPolynomial;
012import edu.jas.poly.GenPolynomialRing;
013import edu.jas.poly.PolyUtil;
014import edu.jas.structure.GcdRingElem;
015import edu.jas.structure.Power;
016import edu.jas.structure.RingElem;
017import edu.jas.structure.RingFactory;
018
019
020/**
021 * Greatest common divisor algorithms with subresultant polynomial remainder
022 * sequence.
023 * @author Heinz Kredel
024 */
025
026public class GreatestCommonDivisorSubres<C extends GcdRingElem<C>> extends GreatestCommonDivisorAbstract<C> {
027
028
029    private static final Logger logger = Logger.getLogger(GreatestCommonDivisorSubres.class);
030
031
032    private boolean debug = logger.isDebugEnabled();
033
034
035    /**
036     * GenPolynomial pseudo remainder. For univariate polynomials.
037     * @param P GenPolynomial.
038     * @param S nonzero GenPolynomial.
039     * @return remainder with ldcf(S)<sup>m</sup> P = quotient * S + remainder.
040     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
041     * @deprecated Use {@link
042     *             edu.jas.poly.PolyUtil#baseDensePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)}
043     *             instead
044     */
045    @Deprecated
046    public GenPolynomial<C> basePseudoRemainder(GenPolynomial<C> P, GenPolynomial<C> S) {
047        return PolyUtil.<C> baseDensePseudoRemainder(P, S);
048    }
049
050
051    /**
052     * GenPolynomial pseudo remainder. For recursive polynomials.
053     * @param P recursive GenPolynomial.
054     * @param S nonzero recursive GenPolynomial.
055     * @return remainder with ldcf(S)<sup>m</sup> P = quotient * S + remainder.
056     * @see edu.jas.poly.GenPolynomial#remainder(edu.jas.poly.GenPolynomial).
057     * @deprecated Use {@link
058     *             edu.jas.poly.PolyUtil#recursiveDensePseudoRemainder(edu.jas.poly.GenPolynomial,edu.jas.poly.GenPolynomial)} instead
059     */
060    @Deprecated
061    public GenPolynomial<GenPolynomial<C>> recursivePseudoRemainder(GenPolynomial<GenPolynomial<C>> P,
062                    GenPolynomial<GenPolynomial<C>> S) {
063        return PolyUtil.<C> recursiveDensePseudoRemainder(P, S);
064    }
065
066
067    /**
068     * Univariate GenPolynomial greatest comon divisor. Uses pseudoRemainder for
069     * remainder.
070     * @param P univariate GenPolynomial.
071     * @param S univariate GenPolynomial.
072     * @return gcd(P,S).
073     */
074    @Override
075    public GenPolynomial<C> baseGcd(GenPolynomial<C> P, GenPolynomial<C> S) {
076        if (S == null || S.isZERO()) {
077            return P;
078        }
079        if (P == null || P.isZERO()) {
080            return S;
081        }
082        if (P.ring.nvar > 1) {
083            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
084        }
085        long e = P.degree(0);
086        long f = S.degree(0);
087        GenPolynomial<C> q;
088        GenPolynomial<C> r;
089        if (f > e) {
090            r = P;
091            q = S;
092            long g = f;
093            f = e;
094            e = g;
095        } else {
096            q = P;
097            r = S;
098        }
099        if (debug) {
100            logger.debug("degrees: e = " + e + ", f = " + f);
101        }
102        r = r.abs();
103        q = q.abs();
104        C a = baseContent(r);
105        C b = baseContent(q);
106        C c = gcd(a, b); // indirection
107        r = divide(r, a); // indirection
108        q = divide(q, b); // indirection
109        if (r.isONE()) {
110            return r.multiply(c);
111        }
112        if (q.isONE()) {
113            return q.multiply(c);
114        }
115        C g = r.ring.getONECoefficient();
116        C h = r.ring.getONECoefficient();
117        GenPolynomial<C> x;
118        C z;
119        while (!r.isZERO()) {
120            long delta = q.degree(0) - r.degree(0);
121            //System.out.println("delta    = " + delta);
122            x = PolyUtil.<C> baseDensePseudoRemainder(q, r);
123            q = r;
124            if (!x.isZERO()) {
125                z = g.multiply(power(P.ring.coFac, h, delta));
126                //System.out.println("z  = " + z);
127                r = x.divide(z);
128                g = q.leadingBaseCoefficient();
129                z = power(P.ring.coFac, g, delta);
130                h = z.divide(power(P.ring.coFac, h, delta - 1));
131                //System.out.println("g  = " + g);
132                //System.out.println("h  = " + h);
133            } else {
134                r = x;
135            }
136        }
137        q = basePrimitivePart(q);
138        return (q.multiply(c)).abs();
139    }
140
141
142    /**
143     * Univariate GenPolynomial recursive greatest comon divisor. Uses
144     * pseudoRemainder for remainder.
145     * @param P univariate recursive GenPolynomial.
146     * @param S univariate recursive GenPolynomial.
147     * @return gcd(P,S).
148     */
149    @Override
150    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateGcd(GenPolynomial<GenPolynomial<C>> P,
151                    GenPolynomial<GenPolynomial<C>> S) {
152        if (S == null || S.isZERO()) {
153            return P;
154        }
155        if (P == null || P.isZERO()) {
156            return S;
157        }
158        if (P.ring.nvar > 1) {
159            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
160        }
161        long e = P.degree(0);
162        long f = S.degree(0);
163        GenPolynomial<GenPolynomial<C>> q;
164        GenPolynomial<GenPolynomial<C>> r;
165        if (f > e) {
166            r = P;
167            q = S;
168            long g = f;
169            f = e;
170            e = g;
171        } else {
172            q = P;
173            r = S;
174        }
175        if (debug) {
176            logger.debug("degrees: e = " + e + ", f = " + f);
177        }
178        r = r.abs();
179        q = q.abs();
180        GenPolynomial<C> a = recursiveContent(r);
181        GenPolynomial<C> b = recursiveContent(q);
182
183        GenPolynomial<C> c = gcd(a, b); // go to recursion
184        //System.out.println("rgcd c = " + c);
185        r = PolyUtil.<C> recursiveDivide(r, a);
186        q = PolyUtil.<C> recursiveDivide(q, b);
187        if (r.isONE()) {
188            return r.multiply(c);
189        }
190        if (q.isONE()) {
191            return q.multiply(c);
192        }
193        GenPolynomial<C> g = r.ring.getONECoefficient();
194        GenPolynomial<C> h = r.ring.getONECoefficient();
195        GenPolynomial<GenPolynomial<C>> x;
196        GenPolynomial<C> z = null;
197        while (!r.isZERO()) {
198            long delta = q.degree(0) - r.degree(0);
199            //System.out.println("rgcd delta = " + delta);
200            x = PolyUtil.<C> recursiveDensePseudoRemainder(q, r);
201            q = r;
202            if (!x.isZERO()) {
203                z = g.multiply(power(P.ring.coFac, h, delta));
204                r = PolyUtil.<C> recursiveDivide(x, z);
205                g = q.leadingBaseCoefficient();
206                z = power(P.ring.coFac, g, delta);
207                h = PolyUtil.<C> basePseudoDivide(z, power(P.ring.coFac, h, delta - 1));
208            } else {
209                r = x;
210            }
211        }
212        q = recursivePrimitivePart(q);
213        return q.abs().multiply(c); //.abs();
214    }
215
216
217    /**
218     * Univariate GenPolynomial resultant. Uses pseudoRemainder for remainder.
219     * @param P univariate GenPolynomial.
220     * @param S univariate GenPolynomial.
221     * @return res(P,S).
222     */
223    @Override
224    public GenPolynomial<C> baseResultant(GenPolynomial<C> P, GenPolynomial<C> S) {
225        if (S == null || S.isZERO()) {
226            return S;
227        }
228        if (P == null || P.isZERO()) {
229            return P;
230        }
231        if (P.ring.nvar > 1) {
232            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
233        }
234        long e = P.degree(0);
235        long f = S.degree(0);
236        GenPolynomial<C> q;
237        GenPolynomial<C> r;
238        if (f > e) {
239            r = P;
240            q = S;
241            long g = f;
242            f = e;
243            e = g;
244        } else {
245            q = P;
246            r = S;
247        }
248        //r = r.abs();
249        //q = q.abs();
250        C a = baseContent(r);
251        C b = baseContent(q);
252        r = divide(r, a); // indirection
253        q = divide(q, b); // indirection
254        RingFactory<C> cofac = P.ring.coFac;
255        C g = cofac.getONE();
256        C h = cofac.getONE();
257        C t = power(cofac, a, e);
258        t = t.multiply(power(cofac, b, f));
259        long s = 1;
260        GenPolynomial<C> x;
261        C z;
262        while (r.degree(0) > 0) {
263            long delta = q.degree(0) - r.degree(0);
264            //System.out.println("delta    = " + delta);
265            if ((q.degree(0) % 2 != 0) && (r.degree(0) % 2 != 0)) {
266                s = -s;
267            }
268            x = PolyUtil.<C> baseDensePseudoRemainder(q, r);
269            //System.out.println("x  = " + x);
270            q = r;
271            if (x.degree(0) > 0) {
272                z = g.multiply(power(cofac, h, delta));
273                //System.out.println("z  = " + z);
274                r = x.divide(z);
275                g = q.leadingBaseCoefficient();
276                z = power(cofac, g, delta);
277                h = z.divide(power(cofac, h, delta - 1));
278            } else {
279                r = x;
280            }
281        }
282        z = power(cofac, r.leadingBaseCoefficient(), q.degree(0));
283        h = z.divide(power(cofac, h, q.degree(0) - 1));
284        z = h.multiply(t);
285        if ( s < 0 ) {
286            z = z.negate();
287        }
288        x = P.ring.getONE().multiply(z);
289        return x;
290    }
291
292
293    /**
294     * Univariate GenPolynomial recursive resultant. Uses pseudoRemainder for
295     * remainder.
296     * @param P univariate recursive GenPolynomial.
297     * @param S univariate recursive GenPolynomial.
298     * @return res(P,S).
299     */
300    @Override
301    public GenPolynomial<GenPolynomial<C>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<C>> P,
302                    GenPolynomial<GenPolynomial<C>> S) {
303        if (S == null || S.isZERO()) {
304            return S;
305        }
306        if (P == null || P.isZERO()) {
307            return P;
308        }
309        if (P.ring.nvar > 1) {
310            throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial");
311        }
312        long e = P.degree(0);
313        long f = S.degree(0);
314        GenPolynomial<GenPolynomial<C>> q;
315        GenPolynomial<GenPolynomial<C>> r;
316        if (f > e) {
317            r = P;
318            q = S;
319            long g = f;
320            f = e;
321            e = g;
322        } else {
323            q = P;
324            r = S;
325        }
326        r = r.abs();
327        q = q.abs();
328        GenPolynomial<C> a = recursiveContent(r);
329        GenPolynomial<C> b = recursiveContent(q);
330        r = PolyUtil.<C> recursiveDivide(r, a);
331        q = PolyUtil.<C> recursiveDivide(q, b);
332        RingFactory<GenPolynomial<C>> cofac = P.ring.coFac;
333        GenPolynomial<C> g = cofac.getONE();
334        GenPolynomial<C> h = cofac.getONE();
335        GenPolynomial<GenPolynomial<C>> x;
336        GenPolynomial<C> t;
337        if (f == 0 && e == 0 && g.ring.nvar > 0) {
338            // if coeffs are multivariate (and non constant)
339            // otherwise it would be 1
340            t = resultant(a, b);
341            x = P.ring.getONE().multiply(t);
342            return x;
343        }
344        t = power(cofac, a, e);
345        t = t.multiply(power(cofac, b, f));
346        long s = 1;
347        GenPolynomial<C> z;
348        while (r.degree(0) > 0) {
349            long delta = q.degree(0) - r.degree(0);
350            //System.out.println("delta    = " + delta);
351            if ((q.degree(0) % 2 != 0) && (r.degree(0) % 2 != 0)) {
352                s = -s;
353            }
354            x = PolyUtil.<C> recursiveDensePseudoRemainder(q, r);
355            //System.out.println("x  = " + x);
356            q = r;
357            if (x.degree(0) > 0) {
358                z = g.multiply(power(P.ring.coFac, h, delta));
359                r = PolyUtil.<C> recursiveDivide(x, z);
360                g = q.leadingBaseCoefficient();
361                z = power(cofac, g, delta);
362                h = PolyUtil.<C> basePseudoDivide(z, power(cofac, h, delta - 1));
363            } else {
364                r = x;
365            }
366        }
367        z = power(cofac, r.leadingBaseCoefficient(), q.degree(0));
368        h = PolyUtil.<C> basePseudoDivide(z, power(cofac, h, q.degree(0) - 1));
369        z = h.multiply(t);
370        if ( s < 0 ) {
371            z = z.negate();
372        }
373        x = P.ring.getONE().multiply(z);
374        return x;
375    }
376
377
378    /**
379     * GenPolynomial base coefficient discriminant.
380     * @param P GenPolynomial.
381     * @return discriminant(P).
382     */
383    public GenPolynomial<C> baseDiscriminant(GenPolynomial<C> P) {
384        if (P == null) {
385            throw new IllegalArgumentException(this.getClass().getName() + " P != null");
386        }
387        if (P.isZERO()) {
388            return P;
389        }
390        GenPolynomialRing<C> pfac = P.ring;
391        if (pfac.nvar > 1) {
392            throw new IllegalArgumentException(this.getClass().getName() + " P not univariate");
393        }
394        C a = P.leadingBaseCoefficient();
395        a = a.inverse();
396        GenPolynomial<C> Pp = PolyUtil.<C> baseDeriviative(P);
397        GenPolynomial<C> res = baseResultant(P, Pp);
398        GenPolynomial<C> disc = res.multiply(a);
399        long n = P.degree(0);
400        n = n * (n - 1);
401        n = n / 2;
402        if (n % 2L != 0L) {
403            disc = disc.negate();
404        }
405        return disc;
406    }
407
408
409    /**
410     * Coefficient power.
411     * @param A coefficient
412     * @param i exponent.
413     * @return A^i.
414     */
415    C power(RingFactory<C> fac, C A, long i) {
416        return Power.<C> power(fac, A, i);
417    }
418
419
420    /**
421     * Polynomial power.
422     * @param A polynomial.
423     * @param i exponent.
424     * @return A^i.
425     */
426    GenPolynomial<C> power(RingFactory<GenPolynomial<C>> fac, GenPolynomial<C> A, long i) {
427        return Power.<GenPolynomial<C>> power(fac, A, i);
428    }
429
430
431}