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