001/* 002 * $Id: RealRootsAbstract.java 4063 2012-07-27 13:05:07Z kredel $ 003 */ 004 005package edu.jas.root; 006 007 008import java.util.ArrayList; 009import java.util.List; 010 011import org.apache.log4j.Logger; 012 013import edu.jas.arith.BigDecimal; 014import edu.jas.arith.BigInteger; 015import edu.jas.arith.BigRational; 016import edu.jas.arith.Rational; 017import edu.jas.poly.GenPolynomial; 018import edu.jas.poly.GenPolynomialRing; 019import edu.jas.poly.PolyUtil; 020import edu.jas.structure.RingElem; 021import edu.jas.structure.RingFactory; 022import edu.jas.structure.UnaryFunctor; 023 024 025/** 026 * Real roots abstract class. 027 * @param <C> coefficient type. 028 * @author Heinz Kredel 029 */ 030public abstract class RealRootsAbstract<C extends RingElem<C> & Rational> implements RealRoots<C> { 031 032 033 private static final Logger logger = Logger.getLogger(RealRootsAbstract.class); 034 035 036 //private static boolean debug = logger.isDebugEnabled(); 037 038 039 /** 040 * Real root bound. With f(M) * f(-M) != 0. 041 * @param f univariate polynomial. 042 * @return M such that -M < root(f) < M. 043 */ 044 public C realRootBound(GenPolynomial<C> f) { 045 if (f == null) { 046 return null; 047 } 048 RingFactory<C> cfac = f.ring.coFac; 049 C M = cfac.getONE(); 050 if (f.isZERO()) { 051 return M; 052 } 053 if (f.isConstant()) { 054 M = f.leadingBaseCoefficient().abs().sum(cfac.getONE()); 055 return M; 056 } 057 C a = f.leadingBaseCoefficient().abs(); 058 for (C c : f.getMap().values()) { 059 C d = c.abs().divide(a); 060 if (M.compareTo(d) < 0) { 061 M = d; 062 } 063 } 064 // works also without this case, only for optimization 065 // to use integral number interval end points 066 // could be obsolete: can fail if real root is in interval [r,r+1] 067 // for too low precision or too big r, since r is approximation 068 /* 069 if (false&&(Object) M instanceof RealAlgebraicNumber) { 070 RealAlgebraicNumber Mr = (RealAlgebraicNumber) M; 071 BigRational r = Mr.magnitude(); 072 //logger.info("rational root bound: " + r); 073 BigInteger i = new BigInteger(r.numerator().divide(r.denominator())); 074 i = i.sum(BigInteger.ONE); 075 //logger.info("integer root bound: " + i); 076 //M = cfac.fromInteger(r.numerator()).divide(cfac.fromInteger(r.denominator())); 077 M = cfac.fromInteger(i.getVal()); 078 M = M.sum(f.ring.coFac.getONE()); 079 logger.info("root bound: " + M); 080 //System.out.println("M = " + M); 081 return M; 082 } 083 */ 084 BigRational r = M.getRational(); 085 logger.info("rational root bound: " + r); 086 BigInteger i = new BigInteger(r.numerator().divide(r.denominator())); 087 i = i.sum(BigInteger.ONE); // ceiling 088 M = cfac.fromInteger(i.getVal()); 089 M = M.sum(f.ring.coFac.getONE()); 090 logger.info("integer root bound: " + M); 091 //System.out.println("M = " + M); 092 return M; 093 } 094 095 096 /** 097 * Magnitude bound. 098 * @param iv interval. 099 * @param f univariate polynomial. 100 * @return B such that |f(c)| < B for c in iv. 101 */ 102 public C magnitudeBound(Interval<C> iv, GenPolynomial<C> f) { 103 if (f == null) { 104 return null; 105 } 106 if (f.isZERO()) { 107 return f.ring.coFac.getONE(); 108 } 109 if (f.isConstant()) { 110 return f.leadingBaseCoefficient().abs(); 111 } 112 GenPolynomial<C> fa = f.map(new UnaryFunctor<C, C>() { 113 114 115 public C eval(C a) { 116 return a.abs(); 117 } 118 }); 119 //System.out.println("fa = " + fa); 120 C M = iv.left.abs(); 121 if (M.compareTo(iv.right.abs()) < 0) { 122 M = iv.right.abs(); 123 } 124 //System.out.println("M = " + M); 125 RingFactory<C> cfac = f.ring.coFac; 126 C B = PolyUtil.<C> evaluateMain(cfac, fa, M); 127 // works also without this case, only for optimization 128 // to use rational number interval end points 129 // can fail if real root is in interval [r,r+1] 130 // for too low precision or too big r, since r is approximation 131 if ((Object) B instanceof RealAlgebraicNumber) { 132 RealAlgebraicNumber Br = (RealAlgebraicNumber) B; 133 BigRational r = Br.magnitude(); 134 B = cfac.fromInteger(r.numerator()).divide(cfac.fromInteger(r.denominator())); 135 } 136 //System.out.println("B = " + B); 137 return B; 138 } 139 140 141 /** 142 * Bi-section point. 143 * @param iv interval with f(left) * f(right) != 0. 144 * @param f univariate polynomial, non-zero. 145 * @return a point c in the interval iv such that f(c) != 0. 146 */ 147 public C bisectionPoint(Interval<C> iv, GenPolynomial<C> f) { 148 if (f == null) { 149 return null; 150 } 151 RingFactory<C> cfac = f.ring.coFac; 152 C two = cfac.fromInteger(2); 153 C c = iv.left.sum(iv.right); 154 c = c.divide(two); 155 if (f.isZERO() || f.isConstant()) { 156 return c; 157 } 158 C m = PolyUtil.<C> evaluateMain(cfac, f, c); 159 while (m.isZERO()) { 160 C d = iv.left.sum(c); 161 d = d.divide(two); 162 if (d.equals(c)) { 163 d = iv.right.sum(c); 164 d = d.divide(two); 165 if (d.equals(c)) { 166 throw new RuntimeException("should not happen " + iv); 167 } 168 } 169 c = d; 170 m = PolyUtil.<C> evaluateMain(cfac, f, c); 171 //System.out.println("c = " + c); 172 } 173 //System.out.println("c = " + c); 174 return c; 175 } 176 177 178 /** 179 * Isolating intervals for the real roots. 180 * @param f univariate polynomial. 181 * @return a list of isolating intervalls for the real roots of f. 182 */ 183 public abstract List<Interval<C>> realRoots(GenPolynomial<C> f); 184 185 186 /** 187 * Isolating intervals for the real roots. 188 * @param f univariate polynomial. 189 * @param eps requested intervals length. 190 * @return a list of isolating intervals v such that |v| < eps. 191 */ 192 public List<Interval<C>> realRoots(GenPolynomial<C> f, C eps) { 193 List<Interval<C>> iv = realRoots(f); 194 return refineIntervals(iv, f, eps); 195 } 196 197 198 /** 199 * Isolating intervals for the real roots. 200 * @param f univariate polynomial. 201 * @param eps requested intervals length. 202 * @return a list of isolating intervals v such that |v| < eps. 203 */ 204 public List<Interval<C>> realRoots(GenPolynomial<C> f, BigRational eps) { 205 C e = f.ring.coFac.parse(eps.toString()); 206 return realRoots(f, e); 207 } 208 209 210 /** 211 * Sign changes on interval bounds. 212 * @param iv root isolating interval with f(left) * f(right) != 0. 213 * @param f univariate polynomial. 214 * @return true if f(left) * f(right) < 0, else false 215 */ 216 public boolean signChange(Interval<C> iv, GenPolynomial<C> f) { 217 if (f == null) { 218 return false; 219 } 220 RingFactory<C> cfac = f.ring.coFac; 221 C l = PolyUtil.<C> evaluateMain(cfac, f, iv.left); 222 C r = PolyUtil.<C> evaluateMain(cfac, f, iv.right); 223 return l.signum() * r.signum() < 0; 224 } 225 226 227 /** 228 * Number of real roots in interval. 229 * @param iv interval with f(left) * f(right) != 0. 230 * @param f univariate polynomial. 231 * @return number of real roots of f in I. 232 */ 233 public abstract long realRootCount(Interval<C> iv, GenPolynomial<C> f); 234 235 236 /** 237 * Half interval. 238 * @param iv root isolating interval with f(left) * f(right) < 0. 239 * @param f univariate polynomial, non-zero. 240 * @return a new interval v such that |v| < |iv|/2. 241 */ 242 public Interval<C> halfInterval(Interval<C> iv, GenPolynomial<C> f) { 243 if (f == null || f.isZERO()) { 244 return iv; 245 } 246 C len = iv.length(); 247 C two = len.factory().fromInteger(2); 248 C eps = len.divide(two); 249 return refineInterval(iv, f, eps); 250 } 251 252 253 /** 254 * Refine interval. 255 * @param iv root isolating interval with f(left) * f(right) < 0. 256 * @param f univariate polynomial, non-zero. 257 * @param eps requested interval length. 258 * @return a new interval v such that |v| < eps. 259 */ 260 public Interval<C> refineInterval(Interval<C> iv, GenPolynomial<C> f, C eps) { 261 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 262 return iv; 263 } 264 if (iv.length().compareTo(eps) < 0) { 265 return iv; 266 } 267 RingFactory<C> cfac = f.ring.coFac; 268 C two = cfac.fromInteger(2); 269 Interval<C> v = iv; 270 while (v.length().compareTo(eps) >= 0) { 271 C c = v.left.sum(v.right); 272 c = c.divide(two); 273 //System.out.println("c = " + c); 274 //c = RootUtil.<C>bisectionPoint(v,f); 275 if (PolyUtil.<C> evaluateMain(cfac, f, c).isZERO()) { 276 v = new Interval<C>(c, c); 277 break; 278 } 279 Interval<C> iv1 = new Interval<C>(v.left, c); 280 if (signChange(iv1, f)) { 281 v = iv1; 282 } else { 283 v = new Interval<C>(c, v.right); 284 } 285 } 286 return v; 287 } 288 289 290 /** 291 * Refine intervals. 292 * @param V list of isolating intervals with f(left) * f(right) < 0. 293 * @param f univariate polynomial, non-zero. 294 * @param eps requested intervals length. 295 * @return a list of new intervals v such that |v| < eps. 296 */ 297 public List<Interval<C>> refineIntervals(List<Interval<C>> V, GenPolynomial<C> f, C eps) { 298 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 299 return V; 300 } 301 List<Interval<C>> IV = new ArrayList<Interval<C>>(); 302 for (Interval<C> v : V) { 303 Interval<C> iv = refineInterval(v, f, eps); 304 IV.add(iv); 305 } 306 return IV; 307 } 308 309 310 /** 311 * Invariant interval for algebraic number sign. 312 * @param iv root isolating interval for f, with f(left) * f(right) < 0. 313 * @param f univariate polynomial, non-zero. 314 * @param g univariate polynomial, gcd(f,g) == 1. 315 * @return v with v a new interval contained in iv such that g(v) != 0. 316 */ 317 public abstract Interval<C> invariantSignInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g); 318 319 320 /** 321 * Real algebraic number sign. 322 * @param iv root isolating interval for f, with f(left) * f(right) < 0, 323 * with iv such that g(iv) != 0. 324 * @param f univariate polynomial, non-zero. 325 * @param g univariate polynomial, gcd(f,g) == 1. 326 * @return sign(g(iv)) . 327 */ 328 public int realIntervalSign(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) { 329 if (g == null || g.isZERO()) { 330 return 0; 331 } 332 if (f == null || f.isZERO() || f.isConstant()) { 333 return g.signum(); 334 } 335 if (g.isConstant()) { 336 return g.signum(); 337 } 338 RingFactory<C> cfac = f.ring.coFac; 339 C c = iv.left.sum(iv.right); 340 c = c.divide(cfac.fromInteger(2)); 341 C ev = PolyUtil.<C> evaluateMain(cfac, g, c); 342 //System.out.println("ev = " + ev); 343 return ev.signum(); 344 } 345 346 347 /** 348 * Real algebraic number sign. 349 * @param iv root isolating interval for f, with f(left) * f(right) < 0. 350 * @param f univariate polynomial, non-zero. 351 * @param g univariate polynomial, gcd(f,g) == 1. 352 * @return sign(g(v)), with v a new interval contained in iv such that g(v) 353 * != 0. 354 */ 355 public int realSign(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) { 356 if (g == null || g.isZERO()) { 357 return 0; 358 } 359 if (f == null || f.isZERO() || f.isConstant()) { 360 return g.signum(); 361 } 362 if (g.isConstant()) { 363 return g.signum(); 364 } 365 Interval<C> v = invariantSignInterval(iv, f, g); 366 return realIntervalSign(v, f, g); 367 } 368 369 370 /** 371 * Invariant interval for algebraic number magnitude. 372 * @param iv root isolating interval for f, with f(left) * f(right) < 0. 373 * @param f univariate polynomial, non-zero. 374 * @param g univariate polynomial, gcd(f,g) == 1. 375 * @param eps length limit for interval length. 376 * @return v with v a new interval contained in iv such that |g(a) - g(b)| 377 * < eps for a, b in v in iv. 378 */ 379 public Interval<C> invariantMagnitudeInterval(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g, 380 C eps) { 381 Interval<C> v = iv; 382 if (g == null || g.isZERO()) { 383 return v; 384 } 385 if (g.isConstant()) { 386 return v; 387 } 388 if (f == null || f.isZERO() || f.isConstant()) { // ? 389 return v; 390 } 391 GenPolynomial<C> gp = PolyUtil.<C> baseDeriviative(g); 392 //System.out.println("g = " + g); 393 //System.out.println("gp = " + gp); 394 C B = magnitudeBound(iv, gp); 395 //System.out.println("B = " + B); 396 397 RingFactory<C> cfac = f.ring.coFac; 398 C two = cfac.fromInteger(2); 399 400 while (B.multiply(v.length()).compareTo(eps) >= 0) { 401 C c = v.left.sum(v.right); 402 c = c.divide(two); 403 Interval<C> im = new Interval<C>(c, v.right); 404 if (signChange(im, f)) { 405 v = im; 406 } else { 407 v = new Interval<C>(v.left, c); 408 } 409 //System.out.println("v = " + v.toDecimal()); 410 } 411 return v; 412 } 413 414 415 /** 416 * Real algebraic number magnitude. 417 * @param iv root isolating interval for f, with f(left) * f(right) < 0, 418 * with iv such that |g(a) - g(b)| < eps for a, b in iv. 419 * @param f univariate polynomial, non-zero. 420 * @param g univariate polynomial, gcd(f,g) == 1. 421 * @return g(iv) . 422 */ 423 public C realIntervalMagnitude(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g) { 424 if (g.isZERO() || g.isConstant()) { 425 return g.leadingBaseCoefficient(); 426 } 427 RingFactory<C> cfac = f.ring.coFac; 428 //if (false) { 429 // C c = iv.left.sum(iv.right); 430 // c = c.divide(cfac.fromInteger(2)); 431 // C ev = PolyUtil.<C> evaluateMain(cfac, g, c); 432 // return ev; 433 //} 434 C evl = PolyUtil.<C> evaluateMain(cfac, g, iv.left); 435 C evr = PolyUtil.<C> evaluateMain(cfac, g, iv.right); 436 C ev = evl; 437 if (evl.compareTo(evr) <= 0) { 438 ev = evr; 439 } 440 //System.out.println("ev = " + ev + ", evl = " + evl + ", evr = " + evr + ", iv = " + iv); 441 return ev; 442 } 443 444 445 /** 446 * Real algebraic number magnitude. 447 * @param iv root isolating interval for f, with f(left) * f(right) < 0. 448 * @param f univariate polynomial, non-zero. 449 * @param g univariate polynomial, gcd(f,g) == 1. 450 * @param eps length limit for interval length. 451 * @return g(iv) . 452 */ 453 public C realMagnitude(Interval<C> iv, GenPolynomial<C> f, GenPolynomial<C> g, C eps) { 454 if (g.isZERO() || g.isConstant()) { 455 return g.leadingBaseCoefficient(); 456 } 457 Interval<C> v = invariantMagnitudeInterval(iv, f, g, eps); 458 return realIntervalMagnitude(v, f, g); 459 } 460 461 462 /** 463 * Approximate real root. 464 * @param iv real root isolating interval with f(left) * f(right) < 0. 465 * @param f univariate polynomial, non-zero. 466 * @param eps requested interval length. 467 * @return a decimal approximation d such that |d-v| < eps, for f(v) = 0, 468 * v real. 469 */ 470 public BigDecimal approximateRoot(Interval<C> iv, GenPolynomial<C> f, C eps) 471 throws NoConvergenceException { 472 if (iv == null) { 473 throw new IllegalArgumentException("null interval not allowed"); 474 } 475 BigDecimal d = iv.toDecimal(); 476 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 477 return d; 478 } 479 if (iv.length().compareTo(eps) < 0) { 480 return d; 481 } 482 BigDecimal left = new BigDecimal(iv.left.getRational()); 483 BigDecimal right = new BigDecimal(iv.right.getRational()); 484 BigDecimal e = new BigDecimal(eps.getRational()); 485 BigDecimal q = new BigDecimal("0.25"); 486 //System.out.println("left = " + left); 487 //System.out.println("right = " + right); 488 e = e.multiply(d); // relative error 489 //System.out.println("e = " + e); 490 BigDecimal dc = BigDecimal.ONE; 491 // polynomials with decimal coefficients 492 GenPolynomialRing<BigDecimal> dfac = new GenPolynomialRing<BigDecimal>(dc, f.ring); 493 GenPolynomial<BigDecimal> df = PolyUtil.<C> decimalFromRational(dfac, f); 494 GenPolynomial<C> fp = PolyUtil.<C> baseDeriviative(f); 495 GenPolynomial<BigDecimal> dfp = PolyUtil.<C> decimalFromRational(dfac, fp); 496 497 // Newton Raphson iteration: x_{n+1} = x_n - f(x_n)/f'(x_n) 498 int i = 0; 499 final int MITER = 50; 500 int dir = 0; 501 while (i++ < MITER) { 502 BigDecimal fx = PolyUtil.<BigDecimal> evaluateMain(dc, df, d); // f(d) 503 if (fx.isZERO()) { 504 return d; 505 } 506 BigDecimal fpx = PolyUtil.<BigDecimal> evaluateMain(dc, dfp, d); // f'(d) 507 if (fpx.isZERO()) { 508 throw new NoConvergenceException("zero deriviative should not happen"); 509 } 510 BigDecimal x = fx.divide(fpx); 511 BigDecimal dx = d.subtract(x); 512 //System.out.println("dx = " + dx); 513 if (d.subtract(dx).abs().compareTo(e) <= 0) { 514 return dx; 515 } 516 while (dx.compareTo(left) < 0 || dx.compareTo(right) > 0) { // dx < left: dx - left < 0 517 // dx > right: dx - right > 0 518 //System.out.println("trying to leave interval"); 519 if (i++ > MITER) { // dx > right: dx - right > 0 520 throw new NoConvergenceException("no convergence after " + i + " steps"); 521 } 522 if (i > MITER / 2 && dir == 0) { 523 BigDecimal sd = new BigDecimal(iv.randomPoint().getRational()); 524 d = sd; 525 x = sd.getZERO(); 526 logger.info("trying new random starting point " + d); 527 i = 0; 528 dir = 1; 529 } 530 if (i > MITER / 2 && dir == 1) { 531 BigDecimal sd = new BigDecimal(iv.randomPoint().getRational()); 532 d = sd; 533 x = sd.getZERO(); 534 logger.info("trying new random starting point " + d); 535 //i = 0; 536 dir = 2; // end 537 } 538 x = x.multiply(q); // x * 1/4 539 dx = d.subtract(x); 540 //System.out.println(" x = " + x); 541 //System.out.println("dx = " + dx); 542 } 543 d = dx; 544 } 545 throw new NoConvergenceException("no convergence after " + i + " steps"); 546 } 547 548 549 /** 550 * Approximate real roots. 551 * @param f univariate polynomial, non-zero. 552 * @param eps requested interval length. 553 * @return a list of decimal approximations d such that |d-v| < eps for 554 * all real v with f(v) = 0. 555 */ 556 public List<BigDecimal> approximateRoots(GenPolynomial<C> f, C eps) { 557 List<Interval<C>> iv = realRoots(f); 558 List<BigDecimal> roots = new ArrayList<BigDecimal>(iv.size()); 559 for (Interval<C> i : iv) { 560 BigDecimal r = null; //approximateRoot(i, f, eps); roots.add(r); 561 while (r == null) { 562 try { 563 r = approximateRoot(i, f, eps); 564 roots.add(r); 565 } catch (NoConvergenceException e) { 566 // fall back to exact algorithm 567 //System.out.println("" + e); 568 C len = i.length(); 569 len = len.divide(f.ring.coFac.fromInteger(1000)); 570 i = refineInterval(i, f, len); 571 logger.info("fall back rootRefinement = " + i); 572 } 573 } 574 } 575 return roots; 576 } 577 578 579 /** 580 * Test if x is an approximate real root. 581 * @param x approximate real root. 582 * @param f univariate polynomial, non-zero. 583 * @param eps requested interval length. 584 * @return true if x is a decimal approximation of a real v with f(v) = 0 585 * with |d-v| < eps, else false. 586 */ 587 public boolean isApproximateRoot(BigDecimal x, GenPolynomial<C> f, C eps) { 588 if (x == null) { 589 throw new IllegalArgumentException("null root not allowed"); 590 } 591 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 592 return true; 593 } 594 BigDecimal e = new BigDecimal(eps.getRational()); 595 e = e.multiply(new BigDecimal("1000")); // relax 596 BigDecimal dc = BigDecimal.ONE; 597 // polynomials with decimal coefficients 598 GenPolynomialRing<BigDecimal> dfac = new GenPolynomialRing<BigDecimal>(dc, f.ring); 599 GenPolynomial<BigDecimal> df = PolyUtil.<C> decimalFromRational(dfac, f); 600 GenPolynomial<C> fp = PolyUtil.<C> baseDeriviative(f); 601 GenPolynomial<BigDecimal> dfp = PolyUtil.<C> decimalFromRational(dfac, fp); 602 // 603 return isApproximateRoot(x, df, dfp, e); 604 } 605 606 607 /** 608 * Test if x is an approximate real root. 609 * @param x approximate real root. 610 * @param f univariate polynomial, non-zero. 611 * @param fp univariate polynomial, non-zero, deriviative of f. 612 * @param eps requested interval length. 613 * @return true if x is a decimal approximation of a real v with f(v) = 0 614 * with |d-v| < eps, else false. 615 */ 616 public boolean isApproximateRoot(BigDecimal x, GenPolynomial<BigDecimal> f, GenPolynomial<BigDecimal> fp, 617 BigDecimal eps) { 618 if (x == null) { 619 throw new IllegalArgumentException("null root not allowed"); 620 } 621 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 622 return true; 623 } 624 BigDecimal dc = BigDecimal.ONE; // only for clarity 625 // f(x) 626 BigDecimal fx = PolyUtil.<BigDecimal> evaluateMain(dc, f, x); 627 //System.out.println("fx = " + fx); 628 if (fx.isZERO()) { 629 return true; 630 } 631 // f'(x) 632 BigDecimal fpx = PolyUtil.<BigDecimal> evaluateMain(dc, fp, x); // f'(d) 633 //System.out.println("fpx = " + fpx); 634 if (fpx.isZERO()) { 635 return false; 636 } 637 BigDecimal d = fx.divide(fpx); 638 if (d.isZERO()) { 639 return true; 640 } 641 if (d.abs().compareTo(eps) <= 0) { 642 return true; 643 } 644 System.out.println("x = " + x); 645 System.out.println("d = " + d); 646 return false; 647 } 648 649 650 /** 651 * Test if each x in R is an approximate real root. 652 * @param R ist of approximate real roots. 653 * @param f univariate polynomial, non-zero. 654 * @param eps requested interval length. 655 * @return true if each x in R is a decimal approximation of a real v with 656 * f(v) = 0 with |d-v| < eps, else false. 657 */ 658 public boolean isApproximateRoot(List<BigDecimal> R, GenPolynomial<C> f, C eps) { 659 if (R == null) { 660 throw new IllegalArgumentException("null root not allowed"); 661 } 662 if (f == null || f.isZERO() || f.isConstant() || eps == null) { 663 return true; 664 } 665 BigDecimal e = new BigDecimal(eps.getRational()); 666 e = e.multiply(new BigDecimal("1000")); // relax 667 BigDecimal dc = BigDecimal.ONE; 668 // polynomials with decimal coefficients 669 GenPolynomialRing<BigDecimal> dfac = new GenPolynomialRing<BigDecimal>(dc, f.ring); 670 GenPolynomial<BigDecimal> df = PolyUtil.<C> decimalFromRational(dfac, f); 671 GenPolynomial<C> fp = PolyUtil.<C> baseDeriviative(f); 672 GenPolynomial<BigDecimal> dfp = PolyUtil.<C> decimalFromRational(dfac, fp); 673 for (BigDecimal x : R) { 674 if (!isApproximateRoot(x, df, dfp, e)) { 675 return false; 676 } 677 } 678 return true; 679 } 680 681}