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