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