Java tutorial
/* GeoGebra - Dynamic Mathematics for Everyone http://www.geogebra.org This file is part of GeoGebra. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation. */ package geogebra.common.kernel; import geogebra.common.kernel.arithmetic.PolyFunction; import geogebra.common.kernel.roots.RealRootAdapter; import geogebra.common.kernel.roots.RealRootDerivAdapter; import geogebra.common.main.App; import java.util.Arrays; import java.util.Comparator; import java.util.Iterator; import java.util.TreeSet; import org.apache.commons.math.analysis.solvers.LaguerreSolver; import org.apache.commons.math.analysis.solvers.UnivariateRealSolver; import org.apache.commons.math.analysis.solvers.UnivariateRealSolverFactory; import org.apache.commons.math.complex.Complex; /** * Class for solving polynomial equations */ public class EquationSolver implements EquationSolverInterface { private static final double LAGUERRE_EPS = 1E-5; private LaguerreSolver laguerreSolver; private UnivariateRealSolver rootFinderBrent, rootFinderNewton; /** * Createsnew equation solver * * @param kernel * kernel */ public EquationSolver(Kernel kernel) { // we need someone to polish our roots // rootPolisher = new RealRoot(); // extrFinder = kernel.getExtremumFinder(); } final public int polynomialRoots(double[] roots, boolean multiple) { int realRoots; /* * update: if roots[n], roots[n-1], ..., we need to determine "real" * degree Darko Drakulic, 28.6.2011. */ int degree = roots.length - 1; for (int i = degree; i >= 0 && roots[i] == 0; i--) degree--; switch (degree) { // degree of polynomial case 0: realRoots = 0; break; case 1: roots[0] = -roots[0] / roots[1]; realRoots = 1; break; case 2: realRoots = solveQuadratic(roots, roots, Kernel.STANDARD_PRECISION); if (multiple && realRoots == 1) { realRoots = 2; roots[1] = roots[0]; } break; case 3: realRoots = solveCubic(roots, roots, Kernel.STANDARD_PRECISION); break; default: realRoots = laguerreAll(roots); } // solveQuadratic and solveCubic may return -1 return Math.max(0, realRoots); } final public int polynomialComplexRoots(double[] real, double complex[]) { int ret = -1; switch (real.length - 1) { // degree of polynomial case 0: ret = 0; break; case 1: real[0] = -real[0] / real[1]; complex[0] = 0; ret = 1; break; case 2: ret = solveQuadraticComplex(real, complex); break; default: ret = laguerreAllComplex(real, complex); } return ret; } /** * Solves the quadratic whose coefficients are in the <code>eqn</code> array * and places the non-complex roots back into the same array, returning the * number of roots. The quadratic solved is represented by the equation: * * <pre> * eqn = {C, B, A}; * ax^2 + bx + c = 0 * </pre> * * A return value of <code>-1</code> is used to distinguish a constant * equation, which might be always 0 or never 0, from an equation that has * no zeroes. * * @param eqn * the array that contains the quadratic coefficients * @return the number of roots, or <code>-1</code> if the equation is a * constant */ final public int solveQuadratic(double eqn[]) { return solveQuadratic(eqn, eqn, Kernel.STANDARD_PRECISION); } final public int solveQuadratic(double eqn[], double res[], double eps) { double a = eqn[2]; double b = eqn[1]; double c = eqn[0]; int roots = 0; if (Math.abs(a) < eps) { // The quadratic parabola has degenerated to a line. if (Math.abs(b) < eps) // The line has degenerated to a constant. return -1; res[roots++] = -c / b; } else if (Math.abs(b) < eps * a) { // a*x^2 + c = 0 double x2 = -c / a; if (Kernel.isZero(x2, eps)) { res[roots++] = 0; } else if (x2 < 0) { // no roots return 0; } else { res[roots++] = Math.sqrt(x2); res[roots++] = -Math.sqrt(x2); } } else { // From Numerical Recipes, 5.6, Quadratic and Cubic Equations double d = b * b - 4.0 * a * c; if (Math.abs(d) < eps * b * b) { res[roots++] = -b / (2.0 * a); } else { if (d < 0.0) // If d < 0.0, then there are no roots return 0; d = Math.sqrt(d); // For accuracy, calculate one root using: // (-b +/- d) / 2a // and the other using: // 2c / (-b +/- d) // Choose the sign of the +/- so that b+d gets larger in // magnitude if (b < 0.0) { d = -d; } double q = (b + d) / -2.0; // We already tested a for being 0 above res[roots++] = q / a; res[roots++] = c / q; } } return roots; } /** * @param real * put coefficients here on input, receive real parts on output * @param complex * imaginary parts on output * @return number of roots */ final public static int solveQuadraticComplex(double real[], double complex[]) { double a = real[2]; double b = real[1]; double c = real[0]; int roots = 0; if (Math.abs(a) < Kernel.STANDARD_PRECISION) { // The quadratic parabola has degenerated to a line. if (Math.abs(b) < Kernel.STANDARD_PRECISION) // The line has degenerated to a constant. return -1; complex[roots] = 0; real[roots++] = -c / b; } else { // From Numerical Recipes, 5.6, Quadratic and Cubic Equations double d = b * b - 4.0 * a * c; if (Math.abs(d) < Kernel.STANDARD_PRECISION) { complex[roots] = 0; real[roots++] = -b / (2.0 * a); } else { if (d < 0.0) { // If d < 0.0, then there are 2 complex roots d = Math.sqrt(-d); complex[0] = d / (2.0 * a); real[0] = -b / (2.0 * a); complex[1] = -complex[0]; real[1] = real[0]; roots = 2; } else { d = Math.sqrt(d); // For accuracy, calculate one root using: // (-b +/- d) / 2a // and the other using: // 2c / (-b +/- d) // Choose the sign of the +/- so that b+d gets larger in // magnitude if (b < 0.0) { d = -d; } double q = (b + d) / -2.0; // We already tested a for being 0 above complex[roots] = 0; real[roots++] = q / a; complex[roots] = 0; real[roots++] = c / q; } } } return roots; } /** * Solves the cubic whose coefficients are in the <code>eqn</code> array and * places the non-complex roots back into the same array, returning the * number of roots. The solved cubic is represented by the equation: * * <pre> * eqn = {c, b, a, d} * dx^3 + ax^2 + bx + c = 0 * </pre> * * A return value of -1 is used to distinguish a constant equation that * might be always 0 or never 0 from an equation that has no zeroes. * * @param eqn * an array containing coefficients for a cubic * @return the number of roots, or -1 if the equation is a constant. */ final public int solveCubic(double eqn[]) { return solveCubic(eqn, eqn, Kernel.STANDARD_PRECISION); } /* * adapted from gsl/poly/solve_cubic.c * * added solveQuadratic() added fixRoots() removed sorting of roots * * solve_cubic.c - finds the real roots of x^3 + a x^2 + b x + c = 0 */ final public int solveCubic(double eqn[], double res[], double eps) { int roots = 0; double d = eqn[3]; if (Math.abs(d) < eps) { // The cubic has degenerated to quadratic (or line or ...). return solveQuadratic(eqn, res, eps); } double a = eqn[2] / d; double b = eqn[1] / d; double c = eqn[0] / d; double q = (a * a - 3 * b); // D(q) = 2aD(a) + 3D(b) double r = (2 * a * a * a - 9 * a * b + 27 * c); // D(r) = (3aa-9b)D(a) // - 9aD(b) + 27D(c) double Q = q / 9; // D(Q) = D(q)/9 double R = r / 54; double Q3 = Q * Q * Q; double R2 = R * R; double CR2 = 729 * r * r; // D(CR2) = 729*2rD(r) ( |1458r(3aa-9b) - // 8748q*2a| + |-9a*1458r -8748q*3| + // |27*1458r| )*kernel.STANDARD_PRECISION double CQ3 = 2916 * q * q * q; // D(CQ3) = 2916*3qD(q) (D(root) ~ // D(2sqrt(Q))= -1/sqrt(Q) D(Q), // |D(root)| |2a+3|/sqrt(9q) // *kernel.STANDARD_PRECISION // Application.debug("Q = "+Q+" R = "+R+" Q3 = "+Q3+" R2 = "+R2+" CR2 = "+CR2+" CQ3 = "+CQ3); // Application.debug(Math.abs(CR2 - CQ3)+""); if (Math.abs(R) < Kernel.STANDARD_PRECISION && Math.abs(Q) < Kernel.STANDARD_PRECISION) // if // (R // == // 0 // && // Q // == // 0) { res[roots++] = -a / 3; res[roots++] = -a / 3; res[roots++] = -a / 3; return 3; } // Michael Borcherds changed to check CR2 equal to CQ3 to first 8 // significant digits // important for eg y=(ax-b)^2(cx-d) // |D(CR2-CQ3)|< (|r(aa-3b) - 4qa| + |-3ar -6q| + |9r|)*13122*sqrt(q) / // |2a+3| *kernel.STANDARD_PRECISION // for simplicity, it (may be) about 10* max(CR2,CR3)/|2a+3| * // kernel.STANDARD_PRECISION // else if (Math.abs(CR2 - CQ3) < Math.max(CR2, CQ3) * // kernel.STANDARD_PRECISION) // // else if (CR2 == CQ3) else if (Math.abs(CR2 - CQ3) < Math.max(CR2, CQ3) * 10 / Math.max(1, Math.abs(2 * a + 3)) * Kernel.STANDARD_PRECISION) // else // if // (CR2 // == // CQ3) { // this test is actually R2 == Q3, written in a form suitable // for exact computation with integers // Due to finite precision some double roots may be missed, and // considered to be a pair of complex roots z = x +/- // kernel.STANDARD_PRECISION // i // close to the real axis. double sqrtQ = Math.sqrt(Q); if (R > 0) { res[roots++] = -2 * sqrtQ - a / 3; res[roots++] = sqrtQ - a / 3; res[roots++] = sqrtQ - a / 3; } else { res[roots++] = -sqrtQ - a / 3; res[roots++] = -sqrtQ - a / 3; res[roots++] = 2 * sqrtQ - a / 3; } return 3; } else if (CR2 < CQ3) // equivalent to R2 < Q3 { double sqrtQ = Math.sqrt(Q); double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ; double theta = Math.acos(R / sqrtQ3); double norm = -2 * sqrtQ; res[roots++] = norm * Math.cos(theta / 3) - a / 3; res[roots++] = norm * Math.cos((theta + 2.0 * Math.PI) / 3) - a / 3; res[roots++] = norm * Math.cos((theta - 2.0 * Math.PI) / 3) - a / 3; // GeoGebra addition // TODO: find a better way to deal with this if (res != eqn) // if res has the same reference as eqn, it's not // possible to fix roots fixRoots(res, eqn); return 3; } else { double sgnR = (R >= 0 ? 1 : -1); double A = -sgnR * Math.pow(Math.abs(R) + Math.sqrt(R2 - Q3), 1.0 / 3.0); double B = Q / A; res[roots++] = A + B - a / 3; return 1; } }// */ /* * This pruning step is necessary since solveCubic uses the cosine function * to calculate the roots when there are 3 of them. Since the cosine method * can have an error of +/- 1E-14 we need to make sure that we don't make * any bad decisions due to an error. * * If the root is not near one of the endpoints, then we will only have a * slight inaccuracy in calculating the x intercept which will only cause a * slightly wrong answer for some points very close to the curve. While the * results in that case are not as accurate as they could be, they are not * disastrously inaccurate either. * * On the other hand, if the error happens near one end of the curve, then * our processing to reject values outside of the t=[0,1] range will fail * and the results of that failure will be disastrous since for an entire * horizontal range of test points, we will either overcount or undercount * the crossings and get a wrong answer for all of them, even when they are * clearly and obviously inside or outside the curve. * * To work around this problem, we try a couple of Newton-Raphson iterations * to see if the true root is closer to the endpoint or further away. If it * is further away, then we can stop since we know we are on the right side * of the endpoint. If we change direction, then either we are now being * dragged away from the endpoint in which case the first condition will * cause us to stop, or we have passed the endpoint and are headed back. In * the second case, we simply evaluate the slope at the endpoint itself and * place ourselves on the appropriate side of it or on it depending on that * result. */ private static void fixRoots(double res[], double eqn[]) { double myepsilon = 1E-5; for (int i = 0; i < 3; i++) { double t = res[i]; if (Math.abs(t) < myepsilon) { res[i] = findZero(t, 0, eqn); } else if (Math.abs(t - 1) < myepsilon) { res[i] = findZero(t, 1, eqn); } } } private static double solveEqn(double eqn[], int order, double t) { double v = eqn[order]; int counter = order; while (--counter >= 0) { v = v * t + eqn[counter]; } return v; } private static double findZero(double init, double target, double eqn[]) { double t = init; double slopeqn[] = { eqn[1], 2 * eqn[2], 3 * eqn[3] }; double slope; double origdelta = 0; double origt = t; while (true) { slope = solveEqn(slopeqn, 2, t); if (slope == 0.0) // At a local minima - must return return t; double y = solveEqn(eqn, 3, t); if (y == 0.0) // Found it! - return it return t; // assert(slope != 0 && y != 0); double delta = -(y / slope); // assert(delta != 0); if (origdelta == 0.0) { origdelta = delta; } if (t < target) { if (delta < 0) return t; } else if (t > target) { if (delta > 0) return t; } else return (delta > 0 ? (target + java.lang.Double.MIN_VALUE) : (target - java.lang.Double.MIN_VALUE)); double newt = t + delta; if (t == newt) // The deltas are so small that we aren't moving... return t; if (delta * origdelta < 0) { // We have reversed our path. int tag = (origt < t ? getTag(target, origt, t) : getTag(target, t, origt)); if (tag != INSIDE) // Local minima found away from target - return the middle return (origt + t) / 2; // Local minima somewhere near target - move to target // and let the slope determine the resulting t. t = target; } else { t = newt; } } } private static final int BELOW = -2; private static final int LOWEDGE = -1; private static final int INSIDE = 0; private static final int HIGHEDGE = 1; private static final int ABOVE = 2; /* * Determine where coord lies with respect to the range from low to high. It * is assumed that low <= high. The return value is one of the 5 values * BELOW, LOWEDGE, INSIDE, HIGHEDGE, or ABOVE. */ private static int getTag(double coord, double low, double high) { if (coord <= low) return (coord < low ? BELOW : LOWEDGE); if (coord >= high) return (coord > high ? ABOVE : HIGHEDGE); return INSIDE; } /* ************************************************* */ private static final double LAGUERRE_START = -1; /** * Calculates all roots of a polynomial given by eqn using Laguerres method. * Polishes roots found. The roots are stored in eqn again. * * @param eqn * coefficients of polynomial */ @SuppressWarnings("deprecation") private int laguerreAll(double[] eqn) { // for fast evaluation of polynomial (used for root polishing) PolyFunction polyFunc = new PolyFunction(eqn); PolyFunction derivFunc = polyFunc.getDerivative(); /* * double estx = 0; try { estx = rootPolisher.newtonRaphson(polyFunc, * LAGUERRE_START); Application.debug("newton estx: " + estx); if * (Double.isNaN(estx)) { estx = LAGUERRE_START; * Application.debug("corrected estx: " + estx); } } catch (Exception e) * {} */ // calc roots with Laguerre method /* * old code using Flanagan library //ComplexPoly poly = new * ComplexPoly(eqn); //Complex [] complexRoots = poly.roots(false, new * Complex(LAGUERRE_START, 0)); // don't polish here */ Complex[] complexRoots = null; try { if (laguerreSolver == null) { laguerreSolver = new LaguerreSolver(); } complexRoots = laguerreSolver.solveAll(eqn, LAGUERRE_START); } catch (Exception e) { System.err.println("EquationSolver.LaguerreSolver: " + e.getLocalizedMessage()); } if (complexRoots == null) complexRoots = new Complex[0]; // sort complexRoots by real part into laguerreRoots double[] laguerreRoots = new double[complexRoots.length]; for (int i = 0; i < laguerreRoots.length; i++) { laguerreRoots[i] = complexRoots[i].getReal(); } Arrays.sort(laguerreRoots); // get the roots from Laguerre method int realRoots = 0; double root; for (int i = 0; i < laguerreRoots.length; i++) { // System.out.println("laguerreRoots[" + i + "] = " + // laguerreRoots[i]); // let's polish all complex roots to get all real roots root = laguerreRoots[i]; // check if root is bounded in intervall [root-eps, root+eps] double left = i == 0 ? root - 1 : (root + laguerreRoots[i - 1]) / 2; double right = i == laguerreRoots.length - 1 ? root + 1 : (root + laguerreRoots[i + 1]) / 2; double f_left = polyFunc.evaluate(left); double f_right = polyFunc.evaluate(right); boolean bounded = f_left * f_right < 0.0; try { if (bounded) { if (rootFinderBrent == null) { UnivariateRealSolverFactory fact = UnivariateRealSolverFactory.newInstance(); rootFinderBrent = fact.newBrentSolver(); } // small f'(root): don't go too far from our laguerre root ! double brentRoot = rootFinderBrent.solve(new RealRootAdapter(polyFunc), left, right, root); if (Math.abs(polyFunc.evaluate(brentRoot)) < Math.abs(polyFunc.evaluate(root))) { root = brentRoot; } // System.out.println("Polish bisectNewtonRaphson: " + // root); } else { if (rootFinderNewton == null) { UnivariateRealSolverFactory fact = UnivariateRealSolverFactory.newInstance(); rootFinderNewton = fact.newNewtonSolver(); } // the root is not bounded: give Mr. Newton a chance double newtonRoot = rootFinderNewton.solve(new RealRootDerivAdapter(polyFunc), left, right, root); if (Math.abs(polyFunc.evaluate(newtonRoot)) < Math.abs(polyFunc.evaluate(root))) { root = newtonRoot; } // System.out.println("Polish newtonRaphson: " + root); } } catch (Exception e) { // Application.debug("Polish FAILED: "); // polishing failed: maybe we have an extremum here // try to find a local extremum try { if (rootFinderBrent == null) { UnivariateRealSolverFactory fact = UnivariateRealSolverFactory.newInstance(); rootFinderBrent = fact.newBrentSolver(); } double brentRoot = rootFinderBrent.solve(new RealRootAdapter(derivFunc), left, right); if (Math.abs(polyFunc.evaluate(brentRoot)) < Math.abs(polyFunc.evaluate(root))) { root = brentRoot; } // System.out.println(" find extremum successfull: " + // root); } catch (Exception ex) { App.debug(ex.getMessage()); } } // check if the found root is really ok double[] val = polyFunc.evaluateDerivFunc(root); // get f(root) and // f'(root) double error = Math.abs(val[0]); // | f(root) | double slope = Math.abs(val[1]); boolean success; if (slope < 1) success = error < LAGUERRE_EPS; else success = error < LAGUERRE_EPS * slope; if (success) { // Application.debug("FOUND ROOT: " + root); eqn[realRoots] = root; realRoots++; } else { // Application.debug("no root: " + root + ", error " + error); } } return realRoots; } /** * Calculates all roots of a polynomial given by eqn using Laguerres method. * Polishes roots found. The roots are stored in eqn again. */ @SuppressWarnings("deprecation") private int laguerreAllComplex(double[] real, double[] complex) { Complex[] complexRoots = null; try { if (laguerreSolver == null) { laguerreSolver = new LaguerreSolver(); } complexRoots = laguerreSolver.solveAll(real, LAGUERRE_START); } catch (Exception e) { App.debug("Problem solving with LaguerreSolver" + e.getLocalizedMessage()); return 0; } // sort by real part & remove duplicates TreeSet<Complex> sortedSet = new TreeSet<Complex>(getComparatorReal()); for (int i = 0; i < complexRoots.length; i++) { sortedSet.add(complexRoots[i]); } int roots = 0; Complex temp; Iterator<Complex> iterator = sortedSet.iterator(); while (iterator.hasNext()) { temp = iterator.next(); real[roots] = temp.getReal(); complex[roots] = temp.getImaginary(); roots++; } return roots; } // avoid too big polynomial coefficients in laguerreAll // private static final int MAX_POLYNOMIAL_COEFFICIENT = 1000; // // // avoid huge coefficients // private void checkCoefficients(double [] eqn) { // double max = Math.abs(eqn[0]); // double temp; // for (int i=1; i < eqn.length; i++) { // temp = Math.abs(eqn[i]); // if (temp > max) max = temp; // } // if (max > MAX_POLYNOMIAL_COEFFICIENT) { // Application.debug("changed coefficients"); // double factor = MAX_POLYNOMIAL_COEFFICIENT/max; // for (int i=0; i < eqn.length; i++) { // eqn[i] *= factor; // } // } // } /* * solve_quartic.c - finds the real roots of x^4 + a x^3 + b x^2 + c x + d = * 0 * * modified from GSL Quartic extension * http://www.network-theory.co.uk/download/gslextras/Quartic/ */ public int solveQuartic(double eqn[], double res[], double eps) { if (Math.abs(eqn[4]) < 0) return solveCubic(eqn, res, Kernel.STANDARD_PRECISION); double a = eqn[3] / eqn[4], b = eqn[2] / eqn[4], c = eqn[1] / eqn[4], d = eqn[0] / eqn[4]; /* * This code is based on a simplification of the algorithm from * zsolve_quartic.c for real roots */ double u[] = new double[3], v[] = new double[3], zarr[] = new double[4]; double aa, pp, qq, rr, rc, sc, tc, mt; double w1r, w1i, w2r, w2i, w3r; double v1, v2, arg, theta; double disc, h; int k1 = 0, k2 = 0; int roots = 0; /* * Deal easily with the cases where the quartic is degenerate. The * ordering of solutions is done explicitly. */ if (0 == b && 0 == c) { if (0 == d) { if (a > 0) { res[roots++] = -a; res[roots++] = 0.0; res[roots++] = 0.0; res[roots++] = 0.0; } else { res[roots++] = 0.0; res[roots++] = 0.0; res[roots++] = 0.0; res[roots++] = -a; } return 4; } else if (0 == a) { if (d > 0) { return 0; } res[roots++] = Math.sqrt(Math.sqrt(-d)); res[roots] = -res[roots - 1]; roots++; return 2; } } if (0.0 == c && 0.0 == d) { res[roots++] = 0.0; res[roots++] = 0.0; double[] res2 = new double[3]; res2[2] = 1.0; res2[1] = a; res2[0] = b; int n = solveQuadratic(res2, res2, eps); res[roots++] = res2[0]; res[roots++] = res2[1]; // if (gsl_poly_solve_quadratic(1.0,a,b,x2,x3)==0) { if (n == 0) { mt = 3; } else { mt = 1; } } else { /* * For non-degenerate solutions, proceed by constructing and solving * the resolvent cubic */ aa = a * a; pp = b - (3.0 / 8.0) * aa; qq = c - (1.0 / 2.0) * a * (b - (1.0 / 4.0) * aa); rr = d - (1.0 / 4.0) * (a * c - (1.0 / 4.0) * aa * (b - (3.0 / 16.0) * aa)); rc = (1.0 / 2.0) * pp; sc = (1.0 / 4.0) * ((1.0 / 4.0) * pp * pp - rr); tc = -((1.0 / 8.0) * qq * (1.0 / 8.0) * qq); /* * This code solves the resolvent cubic in a convenient fashion for * this implementation of the quartic. If there are three real * roots, then they are placed directly into u[]. If two are * complex, then the real root is put into u[0] and the real and * imaginary part of the complex roots are placed into u[1] and * u[2], respectively. Additionally, this calculates the * discriminant of the cubic and puts it into the variable disc. */ { double qcub = (rc * rc - 3 * sc); double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc); double Q = qcub / 9; double R = rcub / 54; double Q3 = Q * Q * Q; double R2 = R * R; double CR2 = 729 * rcub * rcub; double CQ3 = 2916 * qcub * qcub * qcub; disc = (CR2 - CQ3) / 2125764.0; if (0 == R && 0 == Q) { u[0] = -rc / 3; u[1] = -rc / 3; u[2] = -rc / 3; } else if (CR2 == CQ3) { double sqrtQ = Math.sqrt(Q); if (R > 0) { u[0] = -2 * sqrtQ - rc / 3; u[1] = sqrtQ - rc / 3; u[2] = sqrtQ - rc / 3; } else { u[0] = -sqrtQ - rc / 3; u[1] = -sqrtQ - rc / 3; u[2] = 2 * sqrtQ - rc / 3; } } else if (CR2 < CQ3) { double sqrtQ = Math.sqrt(Q); double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ; theta = Math.acos(R / sqrtQ3); if (R / sqrtQ3 >= 1.0) theta = 0.0; { double norm = -2 * sqrtQ; u[0] = norm * Math.cos(theta / 3) - rc / 3; u[1] = norm * Math.cos((theta + 2.0 * Math.PI) / 3) - rc / 3; u[2] = norm * Math.cos((theta - 2.0 * Math.PI) / 3) - rc / 3; } } else { double sgnR = (R >= 0 ? 1 : -1); double modR = Math.abs(R); double sqrt_disc = Math.sqrt(R2 - Q3); double A = -sgnR * Math.pow(modR + sqrt_disc, 1.0 / 3.0); double B = Q / A; double mod_diffAB = Math.abs(A - B); u[0] = A + B - rc / 3; u[1] = -0.5 * (A + B) - rc / 3; u[2] = -(Math.sqrt(3.0) / 2.0) * mod_diffAB; } } /* End of solution to resolvent cubic */ /* * Combine the square roots of the roots of the cubic resolvent * appropriately. Also, calculate 'mt' which designates the nature * of the roots: mt=1 : 4 real roots (disc == 0) mt=2 : 0 real roots * (disc < 0) mt=3 : 2 real roots (disc > 0) */ if (0.0 == disc) u[2] = u[1]; if (0 >= disc) { mt = 2; /* * One would think that we could return 0 here and exit, since * mt=2. However, this assignment is temporary and changes to * mt=1 under certain conditions below. */ v[0] = Math.abs(u[0]); v[1] = Math.abs(u[1]); v[2] = Math.abs(u[2]); v1 = Math.max(Math.max(v[0], v[1]), v[2]); /* Work out which two roots have the largest moduli */ k1 = 0; k2 = 0; if (v1 == v[0]) { k1 = 0; v2 = Math.max(v[1], v[2]); } else if (v1 == v[1]) { k1 = 1; v2 = Math.max(v[0], v[2]); } else { k1 = 2; v2 = Math.max(v[0], v[1]); } if (v2 == v[0]) { k2 = 0; } else if (v2 == v[1]) { k2 = 1; } else { k2 = 2; } if (0.0 <= u[k1]) { w1r = Math.sqrt(u[k1]); w1i = 0.0; } else { w1r = 0.0; w1i = Math.sqrt(-u[k1]); } if (0.0 <= u[k2]) { w2r = Math.sqrt(u[k2]); w2i = 0.0; } else { w2r = 0.0; w2i = Math.sqrt(-u[k2]); } } else { mt = 3; if (0.0 == u[1] && 0.0 == u[2]) { arg = 0.0; } else { arg = Math.sqrt(Math.sqrt(u[1] * u[1] + u[2] * u[2])); } theta = Math.atan2(u[2], u[1]); w1r = arg * Math.cos(theta / 2.0); w1i = arg * Math.sin(theta / 2.0); w2r = w1r; w2i = -w1i; } /* Solve the quadratic to obtain the roots to the quartic */ w3r = qq / 8.0 * (w1i * w2i - w1r * w2r) / (w1i * w1i + w1r * w1r) / (w2i * w2i + w2r * w2r); h = a / 4.0; zarr[0] = w1r + w2r + w3r - h; zarr[1] = -w1r - w2r + w3r - h; zarr[2] = -w1r + w2r - w3r - h; zarr[3] = w1r - w2r - w3r - h; /* Arrange the roots into the variables z0, z1, z2, z3 */ if (2 == mt) { if (u[k1] >= 0 && u[k2] >= 0) { mt = 1; res[roots++] = zarr[0]; res[roots++] = zarr[1]; res[roots++] = zarr[2]; res[roots++] = zarr[3]; } else { return 0; } } else { res[roots++] = zarr[0]; res[roots++] = zarr[1]; } } /* Sort the roots as usual */ if (1 == mt) { /* * Roots are all real, sort them by the real part if (*x0 > *x1) * SWAPD (*x0, *x1); if (*x0 > *x2) SWAPD (*x0, *x2); if (*x0 > *x3) * SWAPD (*x0, *x3); * * if (*x1 > *x2) SWAPD (*x1, *x2); if (*x2 > *x3) { SWAPD (*x2, * *x3); if (*x1 > *x2) SWAPD (*x1, *x2); } */ return 4; } /* * 2 real roots if (*x0 > *x1) SWAPD (*x0, *x1); */ return 2; } /** * Returns a comparator for Complex objects. (sorts on real part coordinate) * If equal does return zero (deletes duplicates!) * * @return comparator for complex numbers */ public static Comparator<Complex> getComparatorReal() { if (comparatorReal == null) { comparatorReal = new Comparator<Complex>() { public int compare(Complex itemA, Complex itemB) { double compReal = itemA.getReal() - itemB.getReal(); if (Kernel.isZero(compReal)) { double compImaginary = itemA.getImaginary() - itemB.getImaginary(); // if real parts equal, sort on imaginary if (!Kernel.isZero(compImaginary)) return compImaginary < 0 ? -1 : +1; // return 0 -> remove duplicates! return 0; } return compReal < 0 ? -1 : +1; } }; } return comparatorReal; } private static Comparator<Complex> comparatorReal; }