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