/*
* $Id$
*/
package edu.jas.root;
import java.util.ArrayList;
import java.util.List;
import org.apache.log4j.Logger;
import edu.jas.arith.BigDecimal;
import edu.jas.arith.BigInteger;
import edu.jas.arith.BigRational;
import edu.jas.arith.Rational;
import edu.jas.poly.GenPolynomial;
import edu.jas.poly.GenPolynomialRing;
import edu.jas.poly.PolyUtil;
import edu.jas.structure.RingElem;
import edu.jas.structure.RingFactory;
import edu.jas.structure.UnaryFunctor;
/**
* Real roots abstract class.
* @param coefficient type.
* @author Heinz Kredel
*/
public abstract class RealRootsAbstract & Rational> implements RealRoots {
private static final Logger logger = Logger.getLogger(RealRootsAbstract.class);
//private static boolean debug = logger.isDebugEnabled();
/**
* Real root bound. With f(M) * f(-M) != 0.
* @param f univariate polynomial.
* @return M such that -M < root(f) < M.
*/
public C realRootBound(GenPolynomial f) {
if (f == null) {
return null;
}
RingFactory cfac = f.ring.coFac;
C M = cfac.getONE();
if (f.isZERO()) {
return M;
}
if (f.isConstant()) {
M = f.leadingBaseCoefficient().abs().sum(cfac.getONE());
return M;
}
C a = f.leadingBaseCoefficient().abs();
for (C c : f.getMap().values()) {
C d = c.abs().divide(a);
if (M.compareTo(d) < 0) {
M = d;
}
}
// works also without this case, only for optimization
// to use integral number interval end points
// could be obsolete: can fail if real root is in interval [r,r+1]
// for too low precision or too big r, since r is approximation
/*
if (false&&(Object) M instanceof RealAlgebraicNumber) {
RealAlgebraicNumber Mr = (RealAlgebraicNumber) M;
BigRational r = Mr.magnitude();
//logger.info("rational root bound: " + r);
BigInteger i = new BigInteger(r.numerator().divide(r.denominator()));
i = i.sum(BigInteger.ONE);
//logger.info("integer root bound: " + i);
//M = cfac.fromInteger(r.numerator()).divide(cfac.fromInteger(r.denominator()));
M = cfac.fromInteger(i.getVal());
M = M.sum(f.ring.coFac.getONE());
logger.info("root bound: " + M);
//System.out.println("M = " + M);
return M;
}
*/
BigRational r = M.getRational();
logger.info("rational root bound: " + r);
BigInteger i = new BigInteger(r.numerator().divide(r.denominator()));
i = i.sum(BigInteger.ONE); // ceiling
M = cfac.fromInteger(i.getVal());
M = M.sum(f.ring.coFac.getONE());
logger.info("integer root bound: " + M);
//System.out.println("M = " + M);
return M;
}
/**
* Magnitude bound.
* @param iv interval.
* @param f univariate polynomial.
* @return B such that |f(c)| < B for c in iv.
*/
@SuppressWarnings("unchecked")
public C magnitudeBound(Interval iv, GenPolynomial f) {
if (f == null) {
return null;
}
if (f.isZERO()) {
return f.ring.coFac.getONE();
}
if (f.isConstant()) {
return f.leadingBaseCoefficient().abs();
}
GenPolynomial fa = f.map(new UnaryFunctor() {
public C eval(C a) {
return a.abs();
}
});
//System.out.println("fa = " + fa);
C M = iv.left.abs();
if (M.compareTo(iv.right.abs()) < 0) {
M = iv.right.abs();
}
//System.out.println("M = " + M);
RingFactory cfac = f.ring.coFac;
C B = PolyUtil. evaluateMain(cfac, fa, M);
// works also without this case, only for optimization
// to use rational number interval end points
// can fail if real root is in interval [r,r+1]
// for too low precision or too big r, since r is approximation
if ((Object) B instanceof RealAlgebraicNumber) {
RealAlgebraicNumber Br = (RealAlgebraicNumber) B;
BigRational r = Br.magnitude();
B = cfac.fromInteger(r.numerator()).divide(cfac.fromInteger(r.denominator()));
}
//System.out.println("B = " + B);
return B;
}
/**
* Bi-section point.
* @param iv interval with f(left) * f(right) != 0.
* @param f univariate polynomial, non-zero.
* @return a point c in the interval iv such that f(c) != 0.
*/
public C bisectionPoint(Interval iv, GenPolynomial f) {
if (f == null) {
return null;
}
RingFactory cfac = f.ring.coFac;
C two = cfac.fromInteger(2);
C c = iv.left.sum(iv.right);
c = c.divide(two);
if (f.isZERO() || f.isConstant()) {
return c;
}
C m = PolyUtil. evaluateMain(cfac, f, c);
while (m.isZERO()) {
C d = iv.left.sum(c);
d = d.divide(two);
if (d.equals(c)) {
d = iv.right.sum(c);
d = d.divide(two);
if (d.equals(c)) {
throw new RuntimeException("should not happen " + iv);
}
}
c = d;
m = PolyUtil. evaluateMain(cfac, f, c);
//System.out.println("c = " + c);
}
//System.out.println("c = " + c);
return c;
}
/**
* Isolating intervals for the real roots.
* @param f univariate polynomial.
* @return a list of isolating intervalls for the real roots of f.
*/
public abstract List> realRoots(GenPolynomial f);
/**
* Isolating intervals for the real roots.
* @param f univariate polynomial.
* @param eps requested intervals length.
* @return a list of isolating intervals v such that |v| < eps.
*/
public List> realRoots(GenPolynomial f, C eps) {
List> iv = realRoots(f);
return refineIntervals(iv, f, eps);
}
/**
* Isolating intervals for the real roots.
* @param f univariate polynomial.
* @param eps requested intervals length.
* @return a list of isolating intervals v such that |v| < eps.
*/
public List> realRoots(GenPolynomial f, BigRational eps) {
C e = f.ring.coFac.parse(eps.toString());
return realRoots(f, e);
}
/**
* Sign changes on interval bounds.
* @param iv root isolating interval with f(left) * f(right) != 0.
* @param f univariate polynomial.
* @return true if f(left) * f(right) < 0, else false
*/
public boolean signChange(Interval iv, GenPolynomial f) {
if (f == null) {
return false;
}
RingFactory cfac = f.ring.coFac;
C l = PolyUtil. evaluateMain(cfac, f, iv.left);
C r = PolyUtil. evaluateMain(cfac, f, iv.right);
return l.signum() * r.signum() < 0;
}
/**
* Number of real roots in interval.
* @param iv interval with f(left) * f(right) != 0.
* @param f univariate polynomial.
* @return number of real roots of f in I.
*/
public abstract long realRootCount(Interval iv, GenPolynomial f);
/**
* Half interval.
* @param iv root isolating interval with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @return a new interval v such that |v| < |iv|/2.
*/
public Interval halfInterval(Interval iv, GenPolynomial f) {
if (f == null || f.isZERO()) {
return iv;
}
C len = iv.length();
C two = len.factory().fromInteger(2);
C eps = len.divide(two);
return refineInterval(iv, f, eps);
}
/**
* Refine interval.
* @param iv root isolating interval with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param eps requested interval length.
* @return a new interval v such that |v| < eps.
*/
public Interval refineInterval(Interval iv, GenPolynomial f, C eps) {
if (f == null || f.isZERO() || f.isConstant() || eps == null) {
return iv;
}
if (iv.length().compareTo(eps) < 0) {
return iv;
}
RingFactory cfac = f.ring.coFac;
C two = cfac.fromInteger(2);
Interval v = iv;
while (v.length().compareTo(eps) >= 0) {
C c = v.left.sum(v.right);
c = c.divide(two);
//System.out.println("c = " + c);
//c = RootUtil.bisectionPoint(v,f);
if (PolyUtil. evaluateMain(cfac, f, c).isZERO()) {
v = new Interval(c, c);
break;
}
Interval iv1 = new Interval(v.left, c);
if (signChange(iv1, f)) {
v = iv1;
} else {
v = new Interval(c, v.right);
}
}
return v;
}
/**
* Refine intervals.
* @param V list of isolating intervals with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param eps requested intervals length.
* @return a list of new intervals v such that |v| < eps.
*/
public List> refineIntervals(List> V, GenPolynomial f, C eps) {
if (f == null || f.isZERO() || f.isConstant() || eps == null) {
return V;
}
List> IV = new ArrayList>();
for (Interval v : V) {
Interval iv = refineInterval(v, f, eps);
IV.add(iv);
}
return IV;
}
/**
* Invariant interval for algebraic number sign.
* @param iv root isolating interval for f, with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param g univariate polynomial, gcd(f,g) == 1.
* @return v with v a new interval contained in iv such that g(v) != 0.
*/
public abstract Interval invariantSignInterval(Interval iv, GenPolynomial f, GenPolynomial g);
/**
* Real algebraic number sign.
* @param iv root isolating interval for f, with f(left) * f(right) < 0,
* with iv such that g(iv) != 0.
* @param f univariate polynomial, non-zero.
* @param g univariate polynomial, gcd(f,g) == 1.
* @return sign(g(iv)) .
*/
public int realIntervalSign(Interval iv, GenPolynomial f, GenPolynomial g) {
if (g == null || g.isZERO()) {
return 0;
}
if (f == null || f.isZERO() || f.isConstant()) {
return g.signum();
}
if (g.isConstant()) {
return g.signum();
}
RingFactory cfac = f.ring.coFac;
C c = iv.left.sum(iv.right);
c = c.divide(cfac.fromInteger(2));
C ev = PolyUtil. evaluateMain(cfac, g, c);
//System.out.println("ev = " + ev);
return ev.signum();
}
/**
* Real algebraic number sign.
* @param iv root isolating interval for f, with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param g univariate polynomial, gcd(f,g) == 1.
* @return sign(g(v)), with v a new interval contained in iv such that g(v)
* != 0.
*/
public int realSign(Interval iv, GenPolynomial f, GenPolynomial g) {
if (g == null || g.isZERO()) {
return 0;
}
if (f == null || f.isZERO() || f.isConstant()) {
return g.signum();
}
if (g.isConstant()) {
return g.signum();
}
Interval v = invariantSignInterval(iv, f, g);
return realIntervalSign(v, f, g);
}
/**
* Invariant interval for algebraic number magnitude.
* @param iv root isolating interval for f, with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param g univariate polynomial, gcd(f,g) == 1.
* @param eps length limit for interval length.
* @return v with v a new interval contained in iv such that |g(a) - g(b)|
* < eps for a, b in v in iv.
*/
public Interval invariantMagnitudeInterval(Interval iv, GenPolynomial f, GenPolynomial g,
C eps) {
Interval v = iv;
if (g == null || g.isZERO()) {
return v;
}
if (g.isConstant()) {
return v;
}
if (f == null || f.isZERO() || f.isConstant()) { // ?
return v;
}
GenPolynomial gp = PolyUtil. baseDeriviative(g);
//System.out.println("g = " + g);
//System.out.println("gp = " + gp);
C B = magnitudeBound(iv, gp);
//System.out.println("B = " + B);
RingFactory cfac = f.ring.coFac;
C two = cfac.fromInteger(2);
while (B.multiply(v.length()).compareTo(eps) >= 0) {
C c = v.left.sum(v.right);
c = c.divide(two);
Interval im = new Interval(c, v.right);
if (signChange(im, f)) {
v = im;
} else {
v = new Interval(v.left, c);
}
//System.out.println("v = " + v.toDecimal());
}
return v;
}
/**
* Real algebraic number magnitude.
* @param iv root isolating interval for f, with f(left) * f(right) < 0,
* with iv such that |g(a) - g(b)| < eps for a, b in iv.
* @param f univariate polynomial, non-zero.
* @param g univariate polynomial, gcd(f,g) == 1.
* @return g(iv) .
*/
public C realIntervalMagnitude(Interval iv, GenPolynomial f, GenPolynomial g) {
if (g.isZERO() || g.isConstant()) {
return g.leadingBaseCoefficient();
}
RingFactory cfac = f.ring.coFac;
//if (false) {
// C c = iv.left.sum(iv.right);
// c = c.divide(cfac.fromInteger(2));
// C ev = PolyUtil. evaluateMain(cfac, g, c);
// return ev;
//}
C evl = PolyUtil. evaluateMain(cfac, g, iv.left);
C evr = PolyUtil. evaluateMain(cfac, g, iv.right);
C ev = evl;
if (evl.compareTo(evr) <= 0) {
ev = evr;
}
//System.out.println("ev = " + ev + ", evl = " + evl + ", evr = " + evr + ", iv = " + iv);
return ev;
}
/**
* Real algebraic number magnitude.
* @param iv root isolating interval for f, with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param g univariate polynomial, gcd(f,g) == 1.
* @param eps length limit for interval length.
* @return g(iv) .
*/
public C realMagnitude(Interval iv, GenPolynomial f, GenPolynomial g, C eps) {
if (g.isZERO() || g.isConstant()) {
return g.leadingBaseCoefficient();
}
Interval v = invariantMagnitudeInterval(iv, f, g, eps);
return realIntervalMagnitude(v, f, g);
}
/**
* Approximate real root.
* @param iv real root isolating interval with f(left) * f(right) < 0.
* @param f univariate polynomial, non-zero.
* @param eps requested interval length.
* @return a decimal approximation d such that |d-v| < eps, for f(v) = 0,
* v real.
*/
public BigDecimal approximateRoot(Interval iv, GenPolynomial f, C eps)
throws NoConvergenceException {
if (iv == null) {
throw new IllegalArgumentException("null interval not allowed");
}
BigDecimal d = iv.toDecimal();
if (f == null || f.isZERO() || f.isConstant() || eps == null) {
return d;
}
if (iv.length().compareTo(eps) < 0) {
return d;
}
BigDecimal left = new BigDecimal(iv.left.getRational());
BigDecimal right = new BigDecimal(iv.right.getRational());
BigRational reps = eps.getRational();
BigDecimal e = new BigDecimal(reps);
BigDecimal q = new BigDecimal("0.25");
//System.out.println("left = " + left);
//System.out.println("right = " + right);
e = e.multiply(d); // relative error
//System.out.println("e = " + e);
BigDecimal dc = BigDecimal.ONE;
// polynomials with decimal coefficients
GenPolynomialRing dfac = new GenPolynomialRing(dc, f.ring);
GenPolynomial df = PolyUtil. decimalFromRational(dfac, f);
GenPolynomial fp = PolyUtil. baseDeriviative(f);
GenPolynomial dfp = PolyUtil. decimalFromRational(dfac, fp);
// Newton Raphson iteration: x_{n+1} = x_n - f(x_n)/f'(x_n)
int i = 0;
final int MITER = 50;
int dir = 0;
while (i++ < MITER) {
BigDecimal fx = PolyUtil. evaluateMain(dc, df, d); // f(d)
if (fx.isZERO()) {
return d;
}
BigDecimal fpx = PolyUtil. evaluateMain(dc, dfp, d); // f'(d)
if (fpx.isZERO()) {
throw new NoConvergenceException("zero deriviative should not happen");
}
BigDecimal x = fx.divide(fpx);
BigDecimal dx = d.subtract(x);
//System.out.println("dx = " + dx + ", d = " + d);
if (d.subtract(dx).abs().compareTo(e) <= 0) {
return dx;
}
while (dx.compareTo(left) < 0 || dx.compareTo(right) > 0) { // dx < left: dx - left < 0
// dx > right: dx - right > 0
//System.out.println("trying to leave interval");
if (i++ > MITER) { // dx > right: dx - right > 0
throw new NoConvergenceException("no convergence after " + i + " steps");
}
if (i > MITER / 2 && dir == 0) {
BigDecimal sd = new BigDecimal(iv.randomPoint().getRational());
d = sd;
x = sd.getZERO();
logger.info("trying new random starting point " + d);
i = 0;
dir = 1;
}
if (i > MITER / 2 && dir == 1) {
BigDecimal sd = new BigDecimal(iv.randomPoint().getRational());
d = sd;
x = sd.getZERO();
logger.info("trying new random starting point " + d);
//i = 0;
dir = 2; // end
}
x = x.multiply(q); // x * 1/4
dx = d.subtract(x);
//System.out.println(" x = " + x);
//System.out.println("dx = " + dx);
}
d = dx;
}
throw new NoConvergenceException("no convergence after " + i + " steps");
}
/**
* Approximate real roots.
* @param f univariate polynomial, non-zero.
* @param eps requested interval length.
* @return a list of decimal approximations d such that |d-v| < eps for
* all real v with f(v) = 0.
*/
public List approximateRoots(GenPolynomial f, C eps) {
List> iv = realRoots(f);
List roots = new ArrayList(iv.size());
for (Interval i : iv) {
BigDecimal r = null; //approximateRoot(i, f, eps); roots.add(r);
while (r == null) {
try {
r = approximateRoot(i, f, eps);
roots.add(r);
} catch (NoConvergenceException e) {
// fall back to exact algorithm
//System.out.println("" + e);
C len = i.length();
len = len.divide(f.ring.coFac.fromInteger(1000));
i = refineInterval(i, f, len);
logger.info("fall back rootRefinement = " + i);
}
}
}
return roots;
}
/**
* Test if x is an approximate real root.
* @param x approximate real root.
* @param f univariate polynomial, non-zero.
* @param eps requested interval length.
* @return true if x is a decimal approximation of a real v with f(v) = 0
* with |d-v| < eps, else false.
*/
public boolean isApproximateRoot(BigDecimal x, GenPolynomial f, C eps) {
if (x == null) {
throw new IllegalArgumentException("null root not allowed");
}
if (f == null || f.isZERO() || f.isConstant() || eps == null) {
return true;
}
BigDecimal e = new BigDecimal(eps.getRational());
e = e.multiply(new BigDecimal("1000")); // relax
BigDecimal dc = BigDecimal.ONE;
// polynomials with decimal coefficients
GenPolynomialRing dfac = new GenPolynomialRing(dc, f.ring);
GenPolynomial df = PolyUtil. decimalFromRational(dfac, f);
GenPolynomial fp = PolyUtil. baseDeriviative(f);
GenPolynomial dfp = PolyUtil. decimalFromRational(dfac, fp);
//
return isApproximateRoot(x, df, dfp, e);
}
/**
* Test if x is an approximate real root.
* @param x approximate real root.
* @param f univariate polynomial, non-zero.
* @param fp univariate polynomial, non-zero, deriviative of f.
* @param eps requested interval length.
* @return true if x is a decimal approximation of a real v with f(v) = 0
* with |d-v| < eps, else false.
*/
public boolean isApproximateRoot(BigDecimal x, GenPolynomial f, GenPolynomial fp,
BigDecimal eps) {
if (x == null) {
throw new IllegalArgumentException("null root not allowed");
}
if (f == null || f.isZERO() || f.isConstant() || eps == null) {
return true;
}
BigDecimal dc = BigDecimal.ONE; // only for clarity
// f(x)
BigDecimal fx = PolyUtil. evaluateMain(dc, f, x);
//System.out.println("fx = " + fx);
if (fx.isZERO()) {
return true;
}
// f'(x)
BigDecimal fpx = PolyUtil. evaluateMain(dc, fp, x); // f'(d)
//System.out.println("fpx = " + fpx);
if (fpx.isZERO()) {
return false;
}
BigDecimal d = fx.divide(fpx);
if (d.isZERO()) {
return true;
}
if (d.abs().compareTo(eps) <= 0) {
return true;
}
System.out.println("x = " + x);
System.out.println("d = " + d);
return false;
}
/**
* Test if each x in R is an approximate real root.
* @param R ist of approximate real roots.
* @param f univariate polynomial, non-zero.
* @param eps requested interval length.
* @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.
*/
public boolean isApproximateRoot(List R, GenPolynomial f, C eps) {
if (R == null) {
throw new IllegalArgumentException("null root not allowed");
}
if (f == null || f.isZERO() || f.isConstant() || eps == null) {
return true;
}
BigDecimal e = new BigDecimal(eps.getRational());
e = e.multiply(new BigDecimal("1000")); // relax
BigDecimal dc = BigDecimal.ONE;
// polynomials with decimal coefficients
GenPolynomialRing dfac = new GenPolynomialRing(dc, f.ring);
GenPolynomial df = PolyUtil. decimalFromRational(dfac, f);
GenPolynomial fp = PolyUtil. baseDeriviative(f);
GenPolynomial dfp = PolyUtil. decimalFromRational(dfac, fp);
for (BigDecimal x : R) {
if (!isApproximateRoot(x, df, dfp, e)) {
return false;
}
}
return true;
}
}
Ads help maintain this website.