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}