List of usage examples for org.apache.commons.math.analysis.polynomials PolynomialFunction PolynomialFunction
public PolynomialFunction(double c[])
From source file:cz.cuni.mff.d3s.spl.example.newton.app.Main.java
public static void main(String[] args) { try {/* w w w.j a v a2s .com*/ System.out.printf("Letting things settle down before actual computation.\n"); Thread.sleep(1000 * 1); } catch (InterruptedException ignored) { } /* This seed gives reasonable data, keep it at that ;-). */ Random random = new Random(2); double[] coefficients = generateCoefficients(random); PolynomialFunction function = new PolynomialFunction(coefficients); NewtonSolver solver = new NewtonSolver(); System.out.printf("Will solve polynomial function of degree %d.\n", function.degree()); inspectClass(solver); try { for (int i = 0; i < WARM_UP_LOOPS; i++) { double result = solver.solve(10000, function, -1000, 1000); } long startTimeNanos = System.nanoTime(); for (int i = 0; i < MEASURED_LOOPS; i++) { double result = solver.solve(10000, function, -1000, 1000); } long endTimeNanos = System.nanoTime(); long runTimeNanos = endTimeNanos - startTimeNanos; long runTimeMillis = runTimeNanos / (1000 * 1000); System.out.printf("%d loops of solving took %dns (%dms).\n", MEASURED_LOOPS, runTimeNanos, runTimeMillis); } catch (MaxIterationsExceededException e) { e.printStackTrace(); } catch (FunctionEvaluationException e) { e.printStackTrace(); } inspectClass(solver); }
From source file:geogebra.common.kernel.implicit.PolynomialUtils.java
/** * calculates the quotient of p/d (no calculation of the remainder is done) * @param p dividend/* w w w . j av a 2 s . c om*/ * @param d divisor * @return quotient of p/d */ public static PolynomialFunction polynomialDivision(PolynomialFunction p, PolynomialFunction d) { return new PolynomialFunction(polynomialDivision(p.getCoefficients(), d.getCoefficients())); }
From source file:geogebra.kernel.AlgoIntersectPolynomialConic.java
@Override protected void compute() { double[] A = c.matrix; if (h.isPolynomialFunction(false)) { PolyFunction pf = h.getFunction().getNumericPolynomialDerivative(0); PolynomialFunction y = new PolynomialFunction(pf.getCoeffs()); PolynomialFunction r = new PolynomialFunction(new double[] { A[2], 2 * A[4], A[0] }); r = r.add(y.multiply(new PolynomialFunction(new double[] { 2 * A[5], 2 * A[3] }))); r = r.add(y.multiply(y.multiply(new PolynomialFunction(new double[] { A[1] })))); //Application.debug("r = "+r.toString()); setRootsPolynomial(r);/*w w w .ja v a2 s . c o m*/ } else { Kernel ker = cons.getKernel(); GeoFunction substituteFunctionX = null; ker.setSilentMode(true); /*try { substituteFunctionX = (GeoFunction) ker.getAlgebraProcessor().processValidExpression( ker.getParser().parseGeoGebraExpression("x"))[0]; } catch (MyError e) { // TODO Auto-generated catch block e.printStackTrace(); } catch (ParseException e) { // TODO Auto-generated catch block e.printStackTrace(); } catch (Exception e) { // TODO Auto-generated catch block e.printStackTrace(); }*/ GeoImplicitPoly iPoly = new GeoImplicitPoly(cons); c.toGeoImplicitPoly(iPoly); GeoFunction paramEquation = new GeoFunction(cons, iPoly, null, h); double nroots = 0; double res[] = new double[2]; if (c.getType() == GeoConicND.CONIC_CIRCLE || c.getType() == GeoConicND.CONIC_ELLIPSE) { nroots = kernel.getEquationSolver().solveQuadratic(new double[] { -A[5] * A[5] + A[1] * A[2], 2 * (A[1] * A[4] - A[3] * A[5]), A[0] * A[1] - A[3] * A[3] }, res); } AlgoRoots algo = null; if (nroots == 2) { //assume res[0]>=res[1] if (res[1] > res[0]) { double temp = res[0]; res[0] = res[1]; res[1] = temp; } algo = new AlgoRoots(cons, paramEquation, new GeoNumeric(cons, Math.max(res[1] - Kernel.MIN_PRECISION, h.getMinParameter())), new GeoNumeric(cons, Math.min(res[0] + Kernel.MIN_PRECISION, h.getMaxParameter()))); } else { algo = new AlgoRoots(cons, paramEquation, new GeoNumeric(cons, h.getMinParameter()), new GeoNumeric(cons, h.getMaxParameter())); } GeoPoint[] points = algo.getRootPoints(); List<double[]> valPairs = new ArrayList<double[]>(); for (int i = 0; i < points.length; i++) { double t = points[i].getX(); valPairs.add(new double[] { t, h.evaluate(t) }); } ker.setSilentMode(false); setPoints(valPairs); return; } }
From source file:geogebra.kernel.AlgoIntersectImplicitpolyParametric.java
@Override protected void compute() { if (p == null || !p.isDefined()) { return;/*from w w w .j a v a 2s . c om*/ } if (f != null) { if (!f.isPolynomialFunction(false) || !f.isDefined) { return; } tx = new PolynomialFunction(new double[] { 0, 1 }); //x=t ty = new PolynomialFunction(f.fun.getNumericPolynomialDerivative(0).getCoeffs()); //y=f(t) } else if (l != null) { if (!l.isDefined()) { return; } //get parametrisation of line double startP[] = new double[2]; l.getInhomPointOnLine(startP); tx = new PolynomialFunction(new double[] { startP[0], l.getY() }); //x=p1+t*r1 ty = new PolynomialFunction(new double[] { startP[1], -l.getX() }); //x=p1+t*r1 } else { return; } PolynomialFunction sum = null; PolynomialFunction zs = null; //Insert x and y (univariat)polynomials via the Horner-scheme double[][] coeff = p.getCoeff(); if (coeff != null) for (int i = coeff.length - 1; i >= 0; i--) { zs = new PolynomialFunction(new double[] { coeff[i][coeff[i].length - 1] }); for (int j = coeff[i].length - 2; j >= 0; j--) { zs = zs.multiply(ty).add(new PolynomialFunction(new double[] { coeff[i][j] }));//y*zs+coeff[i][j]; } if (sum == null) sum = zs; else sum = sum.multiply(tx).add(zs);//sum*x+zs; } setRootsPolynomial(sum); }
From source file:geogebra.common.kernel.algos.AlgoIntersectPolynomialConic.java
@Override public void compute() { double[] A = c.getFlatMatrix(); if (h.isPolynomialFunction(false)) { PolyFunction pf = h.getFunction().getNumericPolynomialDerivative(0, false); PolynomialFunction y = new PolynomialFunction(pf.getCoeffs()); PolynomialFunction r = new PolynomialFunction(new double[] { A[2], 2 * A[4], A[0] }); r = r.add(y.multiply(new PolynomialFunction(new double[] { 2 * A[5], 2 * A[3] }))); r = r.add(y.multiply(y.multiply(new PolynomialFunction(new double[] { A[1] })))); // Application.debug("r = "+r.toString()); setRootsPolynomial(r);/* w w w .ja v a 2s. c om*/ } else { Kernel ker = cons.getKernel(); ker.setSilentMode(true); /* * try { substituteFunctionX = (GeoFunction) * ker.getAlgebraProcessor().processValidExpression( * ker.getParser().parseGeoGebraExpression("x"))[0]; } catch * (MyError e) { // TODO Auto-generated catch block * e.printStackTrace(); } catch (ParseException e) { // TODO * Auto-generated catch block e.printStackTrace(); } catch * (Exception e) { // TODO Auto-generated catch block * e.printStackTrace(); } */ GeoImplicitPoly iPoly = new GeoImplicitPoly(cons); c.toGeoImplicitPoly(iPoly); GeoFunction paramEquation = new GeoFunction(cons, iPoly, null, h); double nroots = 0; double res[] = new double[2]; if (c.getType() == GeoConicNDConstants.CONIC_CIRCLE || c.getType() == GeoConicNDConstants.CONIC_ELLIPSE) { nroots = kernel.getEquationSolver().solveQuadratic(new double[] { -A[5] * A[5] + A[1] * A[2], 2 * (A[1] * A[4] - A[3] * A[5]), A[0] * A[1] - A[3] * A[3] }, res, Kernel.STANDARD_PRECISION); } AlgoRoots algo = null; if (nroots == 2) { // assume res[0]>=res[1] if (res[1] > res[0]) { double temp = res[0]; res[0] = res[1]; res[1] = temp; } algo = new AlgoRoots(cons, paramEquation, new GeoNumeric(cons, Math.max(res[1] - Kernel.MIN_PRECISION, h.getMinParameter())), new GeoNumeric(cons, Math.min(res[0] + Kernel.MIN_PRECISION, h.getMaxParameter()))); } else { algo = new AlgoRoots(cons, paramEquation, new GeoNumeric(cons, h.getMinParameter()), new GeoNumeric(cons, h.getMaxParameter())); } GeoPoint[] points = algo.getRootPoints(); List<double[]> valPairs = new ArrayList<double[]>(); for (int i = 0; i < points.length; i++) { double t = points[i].getX(); valPairs.add(new double[] { t, h.evaluate(t) }); } ker.setSilentMode(false); setPoints(valPairs); return; } }
From source file:geogebra.kernel.implicit.AlgoIntersectImplicitpolyParametric.java
@Override protected void compute() { if (!p.isDefined()) { return;/*from ww w . j a v a 2s .com*/ } double maxT; double minT; if (f != null) { if (!f.isDefined()) { return; } if (!f.isPolynomialFunction(false)) { Kernel ker = cons.getKernel(); ker.setSilentMode(true); GeoFunction paramEquation = new GeoFunction(cons, p, null, f); AlgoRoots algo = new AlgoRoots(cons, paramEquation, new GeoNumeric(cons, f.getMinParameter()), new GeoNumeric(cons, f.getMaxParameter())); GeoPoint[] points = algo.getRootPoints(); List<double[]> valPairs = new ArrayList<double[]>(); for (int i = 0; i < points.length; i++) { double t = points[i].getX(); valPairs.add(new double[] { t, f.evaluate(t) }); } ker.setSilentMode(false); setPoints(valPairs); return; } tx = new PolynomialFunction(new double[] { 0, 1 }); //x=t ty = new PolynomialFunction(f.getFunction().getNumericPolynomialDerivative(0).getCoeffs()); //y=f(t) maxT = f.getMaxParameter(); minT = f.getMinParameter(); } else if (l != null) { if (!l.isDefined()) { points.adjustOutputSize(0); return; } //get parametrisation of line double startP[] = new double[2]; l.getInhomPointOnLine(startP); tx = new PolynomialFunction(new double[] { startP[0], l.getY() }); //x=p1+t*r1 ty = new PolynomialFunction(new double[] { startP[1], -l.getX() }); //y=p2+t*r2 maxT = l.getMaxParameter(); minT = l.getMinParameter(); if (l.getParentAlgorithm() instanceof AlgoTangentImplicitpoly) { tangentPoints = ((AlgoTangentImplicitpoly) l.getParentAlgorithm()).getTangentPoints(); } } else { return; } PolynomialFunction sum = null; PolynomialFunction zs = null; //Insert x and y (univariat)polynomials via the Horner-scheme double[][] coeff = p.getCoeff(); if (coeff != null) for (int i = coeff.length - 1; i >= 0; i--) { zs = new PolynomialFunction(new double[] { coeff[i][coeff[i].length - 1] }); for (int j = coeff[i].length - 2; j >= 0; j--) { zs = zs.multiply(ty).add(new PolynomialFunction(new double[] { coeff[i][j] }));//y*zs+coeff[i][j]; } if (sum == null) sum = zs; else sum = sum.multiply(tx).add(zs);//sum*x+zs; } setRootsPolynomialWithinRange(sum, minT, maxT); mergeWithTangentPoints(); }
From source file:geogebra.common.kernel.implicit.AlgoIntersectImplicitpolyParametric.java
@Override public void compute() { if (!p.isDefined()) { return;/* ww w. j ava2 s .co m*/ } double maxT; double minT; if (f != null) { if (!f.isDefined()) { return; } if (!f.isPolynomialFunction(false)) { Kernel ker = cons.getKernel(); ker.setSilentMode(true); GeoFunction paramEquation = new GeoFunction(cons, p, null, f); AlgoRoots algo = new AlgoRoots(cons, paramEquation, new GeoNumeric(cons, f.getMinParameter()), new GeoNumeric(cons, f.getMaxParameter())); GeoPoint[] rootPoints = algo.getRootPoints(); List<double[]> valPairs = new ArrayList<double[]>(); for (int i = 0; i < rootPoints.length; i++) { double t = rootPoints[i].getX(); valPairs.add(new double[] { t, f.evaluate(t) }); } ker.setSilentMode(false); setPoints(valPairs); return; } tx = new PolynomialFunction(new double[] { 0, 1 }); //x=t ty = new PolynomialFunction(f.getFunction().getNumericPolynomialDerivative(0, false).getCoeffs()); //y=f(t) maxT = f.getMaxParameter(); minT = f.getMinParameter(); } else if (l != null) { if (!l.isDefined()) { points.adjustOutputSize(0); return; } //get parametrisation of line double startP[] = new double[2]; l.getInhomPointOnLine(startP); tx = new PolynomialFunction(new double[] { startP[0], l.getY() }); //x=p1+t*r1 ty = new PolynomialFunction(new double[] { startP[1], -l.getX() }); //y=p2+t*r2 maxT = l.getMaxParameter(); minT = l.getMinParameter(); if (l.getParentAlgorithm() instanceof AlgoTangentImplicitpoly) { tangentPoints = ((AlgoTangentImplicitpoly) l.getParentAlgorithm()).getTangentPoints(); } } else { return; } PolynomialFunction sum = null; PolynomialFunction zs = null; //Insert x and y (univariat)polynomials via the Horner-scheme double[][] coeff = p.getCoeff(); if (coeff != null) for (int i = coeff.length - 1; i >= 0; i--) { zs = new PolynomialFunction(new double[] { coeff[i][coeff[i].length - 1] }); for (int j = coeff[i].length - 2; j >= 0; j--) { zs = zs.multiply(ty).add(new PolynomialFunction(new double[] { coeff[i][j] }));//y*zs+coeff[i][j]; } if (sum == null) sum = zs; else sum = sum.multiply(tx).add(zs);//sum*x+zs; } setRootsPolynomialWithinRange(sum, minT, maxT); mergeWithTangentPoints(); }
From source file:geogebra.kernel.implicit.AlgoIntersectImplicitpolys.java
@Override protected void compute() { if (c1 != null) { p2 = new GeoImplicitPoly(c1); }// w w w . j a v a 2 s. co m if (valPairs == null) { valPairs = new LinkedList<double[]>(); } else { valPairs.clear(); } /* * New approach: calculating determinant of Sylvester-matrix to get resolvent * */ // Application.debug("p1="+p1); // Application.debug("p2="+p2); GeoImplicitPoly a = p1, b = p2; if (p1.getDegX() < p2.getDegX()) { a = p2; b = p1; } int m = a.getDegX(); int n = b.getDegX(); //calculate the reduced Sylvester matrix. Complexity will be O(mnpq + m^2nq^2 + n^3pq) //where p=a.getDegY(), q=b.getDegY() //we should minimize m^2 n q^2 by choosing to use polyX or polyY univarType. // int q = a.getDegY(); PolynomialFunction[][] mat = new PolynomialFunction[n][n]; PolynomialFunction[] aNew = new PolynomialFunction[m + n]; PolynomialFunction[] bPolys = new PolynomialFunction[n + 1]; for (int i = 0; i <= n; ++i) bPolys[i] = new PolynomialFunction(b.getCoeff()[i]); for (int i = 0; i < n - 1; ++i) aNew[i] = new PolynomialFunction(new double[] { 0 }); for (int i = n - 1; i < n + m; ++i) aNew[i] = new PolynomialFunction(a.getCoeff()[i - n + 1]); int leadIndex = n + m - 1; //Note: leadIndex of (n+1+t)-th row is equal to X-degree of b, + t. Use //this row to help eliminate aNew[leadIndex]. while (leadIndex >= 2 * n) { if (!(aNew[leadIndex].degree() == 0 && aNew[leadIndex].getCoefficients()[0] == 0)) { for (int j = n - 1; j < leadIndex - n; ++j) aNew[j] = aNew[j].multiply(bPolys[n]); for (int j = leadIndex - n; j < leadIndex; ++j) aNew[j] = aNew[j].multiply(bPolys[n]) .subtract(bPolys[j - leadIndex + n].multiply(aNew[leadIndex])); } --leadIndex; } while (leadIndex >= n) { if (!(aNew[leadIndex].degree() == 0 && aNew[leadIndex].getCoefficients()[0] == 0)) { for (int j = leadIndex - n; j < leadIndex; ++j) aNew[j] = aNew[j].multiply(bPolys[n]) .subtract(bPolys[j - leadIndex + n].multiply(aNew[leadIndex])); } for (int j = 0; j < n; ++j) mat[2 * n - 1 - leadIndex][j] = new PolynomialFunction(aNew[leadIndex - n + j].getCoefficients()); --leadIndex; } //avoid too large coefficients //test case: a: -5 x?+ x+ y = 0m, b: -20 x+2 x+2 x+2 y+4 y = 0 //without reducing coefficients, we get three intersection points: // (0.00000185192649, -0.000000925965389), (0.475635148394481, 0.172245588226639), (2.338809137914722, -12.005665890026151) //after reducing coefficients, we have one more: the tangent point (0.99999997592913, 1.999999891681086) for (int i = 0; i < n; ++i) { double largestCoeff = 0; double reduceFactor = 1; for (int j = 0; j < n; ++j) { for (int k = 0; k < mat[i][j].getCoefficients().length; ++k) { largestCoeff = Math.max(Math.abs(mat[i][j].getCoefficients()[k]), largestCoeff); } } while (largestCoeff > 10) { reduceFactor *= 0.1; largestCoeff *= 0.1; } if (reduceFactor != 1) { for (int j = 0; j < n; ++j) { mat[i][j] = mat[i][j].multiply(new PolynomialFunction(new double[] { reduceFactor })); } } } //Calculate Sylvester matrix by definition. Complexity will be O((m+n)^3 * pq) //where p=a.getDegY(), q=b.getDegY() /* PolynomialFunction[][] mat=new PolynomialFunction[m+n][m+n]; for (int i = 0; i<n; ++i) { for (int j = 0; j<i; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); for (int j = i; j<= i+m; ++j) mat[i][j] = new PolynomialFunction(a.getCoeff()[j-i]); for (int j = i+m+1; j<n+m; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); } for (int i = n; i<m+n; ++i) { for (int j = 0; j<i-n; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); for (int j = i-n; j<= i; ++j) mat[i][j] = new PolynomialFunction(b.getCoeff()[j-i+n]); for (int j = i+1; j<n+m; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); } */ //old code /*PolynomialFunction[][] mat=new PolynomialFunction[n][n]; for (int i=0;i<n;i++){ for (int j=0;j<n;j++){ mat[i][j]=new PolynomialFunction(new double[]{0}); for (int k=Math.max(0, i-j);k<=Math.min(i, m+i-j);k++){ PolynomialFunction p=new PolynomialFunction(b.getCoeff()[k]); mat[i][j]=mat[i][j].add(p.multiply(new PolynomialFunction(a.getCoeff()[m+i-k-j]))); } for (int k=Math.max(0, i+m-j-n);k<=Math.min(i, m+i-j);k++){ PolynomialFunction p=new PolynomialFunction(a.getCoeff()[k]); mat[i][j]=mat[i][j].subtract(p.multiply(new PolynomialFunction(b.getCoeff()[m+i-k-j]))); } } }*/ // Application.debug(Arrays.deepToString(mat)); //Gau-Bareiss for calculating the determinant PolynomialFunction c = new PolynomialFunction(new double[] { 1 }); PolynomialFunction det = null; for (int k = 0; k < n - 1; k++) { int r = 0; double glc = 0; //greatest leading coefficient for (int i = k; i < n; i++) { double lc = PolynomialUtils.getLeadingCoeff(mat[i][k]); if (!Kernel.isZero(lc)) { if (Math.abs(lc) > Math.abs(glc)) { glc = lc; r = i; } } } if (Kernel.isZero(glc)) { det = new PolynomialFunction(new double[] { 0 }); break; } else if (r > k) { for (int j = k; j < n; j++) { //exchange functions PolynomialFunction temp = mat[r][j]; mat[r][j] = mat[k][j]; mat[k][j] = temp; } } for (int i = k + 1; i < n; i++) { for (int j = k + 1; j < n; j++) { PolynomialFunction t1 = mat[i][j].multiply(mat[k][k]); PolynomialFunction t2 = mat[i][k].multiply(mat[k][j]); PolynomialFunction t = t1.subtract(t2); mat[i][j] = PolynomialUtils.polynomialDivision(t, c); } } c = mat[k][k]; } if (det == null) det = mat[n - 1][n - 1]; // Application.debug("resultante = "+det); univarType = PolyY; double roots[] = det.getCoefficients(); int nrRealRoots = 0; if (roots.length > 1) nrRealRoots = getRoots(roots, eqnSolver); double[][] coeff; double[] newCoeff; if (univarType == PolyX) { if (p1.getDegY() < p2.getDegY()) { coeff = p1.getCoeff(); newCoeff = new double[p1.getDegY() + 1]; } else { coeff = p2.getCoeff(); newCoeff = new double[p2.getDegY() + 1]; } } else { if (p1.getDegX() < p2.getDegX()) { coeff = p1.getCoeff(); newCoeff = new double[p1.getDegX() + 1]; } else { coeff = p2.getCoeff(); newCoeff = new double[p2.getDegX() + 1]; } } for (int k = 0; k < nrRealRoots; k++) { double t = roots[k]; if (univarType == PolyX) { for (int j = 0; j < newCoeff.length; j++) { newCoeff[j] = 0; } for (int i = coeff.length - 1; i >= 0; i--) { for (int j = 0; j < coeff[i].length; j++) { newCoeff[j] = newCoeff[j] * t + coeff[i][j]; } for (int j = coeff[i].length; j < newCoeff.length; j++) { newCoeff[j] = newCoeff[j] * t; } } } else { for (int i = 0; i < coeff.length; i++) { newCoeff[i] = 0; for (int j = coeff[i].length - 1; j >= 0; j--) { newCoeff[i] = newCoeff[i] * t + coeff[i][j]; } } } int nr = getRoots(newCoeff, eqnSolver); for (int i = 0; i < nr; i++) { double[] pair = new double[2]; if (univarType == PolyX) { pair[0] = t; pair[1] = newCoeff[i]; } else { pair[0] = newCoeff[i]; pair[1] = t; } if (PolynomialUtils.rootPolishing(pair, p1, p2)) insert(pair); } } if (hints != null) { for (int i = 0; i < hints.size(); i++) { double[] pair = new double[2]; GeoPoint g = hints.get(i); if (g.isDefined() && !Kernel.isZero(g.getZ())) { pair[0] = g.getX() / g.getZ(); pair[1] = g.getY() / g.getZ(); } } } setPoints(valPairs); /* [end new] List<GenPolynomial<BigRational>> polynomials = new ArrayList<GenPolynomial<BigRational>>(); polynomials.add(p1.toGenPolynomial()); polynomials.add(p2.toGenPolynomial()); // Application.debug("dp1: {"+p1.getDegX()+","+p1.getDegY()+"} dp2: {"+p2.getDegX()+","+p2.getDegY()+"}"); // Application.debug("size: "+polynomials.size()); // Application.debug("p: "+polynomials); GroebnerBase<BigRational> gb = GBFactory.getImplementation(BigRational.ONE); List<GenPolynomial<BigRational>> G=gb.GB(polynomials); //G=gb.minimalGB(G); Application.debug("Grbner Basis: "+G); boolean[] var=new boolean[2]; var[0]=var[1]=true; setRootsPolynomial(GeoImplicitPoly.getUnivariatPoly(G,var)); if (var[0]) univarType=0; else univarType=1; */ }
From source file:geogebra.common.kernel.implicit.AlgoIntersectImplicitpolys.java
@Override public void compute() { if (c1 != null) { p2 = new GeoImplicitPoly(c1); }//from ww w. java2s . co m if (valPairs == null) { valPairs = new LinkedList<double[]>(); } else { valPairs.clear(); } /* * New approach: calculating determinant of Sylvester-matrix to get resolvent * */ // Application.debug("p1="+p1); // Application.debug("p2="+p2); GeoImplicitPoly a = p1, b = p2; if (p1.getDegX() < p2.getDegX()) { a = p2; b = p1; } int m = a.getDegX(); int n = b.getDegX(); //calculate the reduced Sylvester matrix. Complexity will be O(mnpq + m^2nq^2 + n^3pq) //where p=a.getDegY(), q=b.getDegY() //we should minimize m^2 n q^2 by choosing to use polyX or polyY univarType. // int q = a.getDegY(); PolynomialFunction[][] mat = new PolynomialFunction[n][n]; PolynomialFunction[] aNew = new PolynomialFunction[m + n]; PolynomialFunction[] bPolys = new PolynomialFunction[n + 1]; for (int i = 0; i <= n; ++i) bPolys[i] = new PolynomialFunction(b.getCoeff()[i]); for (int i = 0; i < n - 1; ++i) aNew[i] = new PolynomialFunction(new double[] { 0 }); for (int i = n - 1; i < n + m; ++i) aNew[i] = new PolynomialFunction(a.getCoeff()[i - n + 1]); int leadIndex = n + m - 1; //Note: leadIndex of (n+1+t)-th row is equal to X-degree of b, + t. Use //this row to help eliminate aNew[leadIndex]. while (leadIndex >= 2 * n) { if (!(aNew[leadIndex].degree() == 0 && aNew[leadIndex].getCoefficients()[0] == 0)) { for (int j = n - 1; j < leadIndex - n; ++j) aNew[j] = aNew[j].multiply(bPolys[n]); for (int j = leadIndex - n; j < leadIndex; ++j) aNew[j] = aNew[j].multiply(bPolys[n]) .subtract(bPolys[j - leadIndex + n].multiply(aNew[leadIndex])); } --leadIndex; } while (leadIndex >= n) { if (!(aNew[leadIndex].degree() == 0 && aNew[leadIndex].getCoefficients()[0] == 0)) { for (int j = leadIndex - n; j < leadIndex; ++j) aNew[j] = aNew[j].multiply(bPolys[n]) .subtract(bPolys[j - leadIndex + n].multiply(aNew[leadIndex])); } for (int j = 0; j < n; ++j) mat[2 * n - 1 - leadIndex][j] = new PolynomialFunction(aNew[leadIndex - n + j].getCoefficients()); --leadIndex; } //avoid too large coefficients //test case: a: -5 x?+ x+ y = 0m, b: -20 x+2 x+2 x+2 y+4 y = 0 //without reducing coefficients, we get three intersection points: // (0.00000185192649, -0.000000925965389), (0.475635148394481, 0.172245588226639), (2.338809137914722, -12.005665890026151) //after reducing coefficients, we have one more: the tangent point (0.99999997592913, 1.999999891681086) for (int i = 0; i < n; ++i) { double largestCoeff = 0; double reduceFactor = 1; for (int j = 0; j < n; ++j) { for (int k = 0; k < mat[i][j].getCoefficients().length; ++k) { largestCoeff = Math.max(Math.abs(mat[i][j].getCoefficients()[k]), largestCoeff); } } while (largestCoeff > 10) { reduceFactor *= 0.1; largestCoeff *= 0.1; } if (reduceFactor != 1) { for (int j = 0; j < n; ++j) { mat[i][j] = mat[i][j].multiply(new PolynomialFunction(new double[] { reduceFactor })); } } } //Calculate Sylvester matrix by definition. Complexity will be O((m+n)^3 * pq) //where p=a.getDegY(), q=b.getDegY() /* PolynomialFunction[][] mat=new PolynomialFunction[m+n][m+n]; for (int i = 0; i<n; ++i) { for (int j = 0; j<i; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); for (int j = i; j<= i+m; ++j) mat[i][j] = new PolynomialFunction(a.getCoeff()[j-i]); for (int j = i+m+1; j<n+m; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); } for (int i = n; i<m+n; ++i) { for (int j = 0; j<i-n; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); for (int j = i-n; j<= i; ++j) mat[i][j] = new PolynomialFunction(b.getCoeff()[j-i+n]); for (int j = i+1; j<n+m; ++j) mat[i][j] = new PolynomialFunction(new double[]{0}); } */ //old code /*PolynomialFunction[][] mat=new PolynomialFunction[n][n]; for (int i=0;i<n;i++){ for (int j=0;j<n;j++){ mat[i][j]=new PolynomialFunction(new double[]{0}); for (int k=Math.max(0, i-j);k<=Math.min(i, m+i-j);k++){ PolynomialFunction p=new PolynomialFunction(b.getCoeff()[k]); mat[i][j]=mat[i][j].add(p.multiply(new PolynomialFunction(a.getCoeff()[m+i-k-j]))); } for (int k=Math.max(0, i+m-j-n);k<=Math.min(i, m+i-j);k++){ PolynomialFunction p=new PolynomialFunction(a.getCoeff()[k]); mat[i][j]=mat[i][j].subtract(p.multiply(new PolynomialFunction(b.getCoeff()[m+i-k-j]))); } } }*/ // Application.debug(Arrays.deepToString(mat)); //Gau-Bareiss for calculating the determinant PolynomialFunction c = new PolynomialFunction(new double[] { 1 }); PolynomialFunction det = null; for (int k = 0; k < n - 1; k++) { int r = 0; double glc = 0; //greatest leading coefficient for (int i = k; i < n; i++) { double lc = PolynomialUtils.getLeadingCoeff(mat[i][k]); if (!Kernel.isZero(lc)) { if (Math.abs(lc) > Math.abs(glc)) { glc = lc; r = i; } } } if (Kernel.isZero(glc)) { det = new PolynomialFunction(new double[] { 0 }); break; } else if (r > k) { for (int j = k; j < n; j++) { //exchange functions PolynomialFunction temp = mat[r][j]; mat[r][j] = mat[k][j]; mat[k][j] = temp; } } for (int i = k + 1; i < n; i++) { for (int j = k + 1; j < n; j++) { PolynomialFunction t1 = mat[i][j].multiply(mat[k][k]); PolynomialFunction t2 = mat[i][k].multiply(mat[k][j]); PolynomialFunction t = t1.subtract(t2); mat[i][j] = PolynomialUtils.polynomialDivision(t, c); } } c = mat[k][k]; } if (det == null) det = mat[n - 1][n - 1]; // Application.debug("resultante = "+det); univarType = PolyY; double roots[] = det.getCoefficients(); // roots[0]-=0.001; int nrRealRoots = 0; if (roots.length > 1) nrRealRoots = getNearRoots(roots, eqnSolver, 1E-1);//getRoots(roots,eqnSolver); double[][] coeff; double[] newCoeff; if (univarType == PolyX) { if (p1.getDegY() < p2.getDegY()) { coeff = p1.getCoeff(); newCoeff = new double[p1.getDegY() + 1]; } else { coeff = p2.getCoeff(); newCoeff = new double[p2.getDegY() + 1]; } } else { if (p1.getDegX() < p2.getDegX()) { coeff = p1.getCoeff(); newCoeff = new double[p1.getDegX() + 1]; } else { coeff = p2.getCoeff(); newCoeff = new double[p2.getDegX() + 1]; } } for (int k = 0; k < nrRealRoots; k++) { double t = roots[k]; if (univarType == PolyX) { for (int j = 0; j < newCoeff.length; j++) { newCoeff[j] = 0; } for (int i = coeff.length - 1; i >= 0; i--) { for (int j = 0; j < coeff[i].length; j++) { newCoeff[j] = newCoeff[j] * t + coeff[i][j]; } for (int j = coeff[i].length; j < newCoeff.length; j++) { newCoeff[j] = newCoeff[j] * t; } } } else { for (int i = 0; i < coeff.length; i++) { newCoeff[i] = 0; for (int j = coeff[i].length - 1; j >= 0; j--) { newCoeff[i] = newCoeff[i] * t + coeff[i][j]; } } } int nr = getNearRoots(newCoeff, eqnSolver, 1E-1);//getRoots(newCoeff,eqnSolver); for (int i = 0; i < nr; i++) { double[] pair = new double[2]; if (univarType == PolyX) { pair[0] = t; pair[1] = newCoeff[i]; } else { pair[0] = newCoeff[i]; pair[1] = t; } if (PolynomialUtils.rootPolishing(pair, p1, p2)) insert(pair); } } if (hints != null) { for (int i = 0; i < hints.size(); i++) { double[] pair = new double[2]; GeoPoint g = hints.get(i); if (g.isDefined() && !Kernel.isZero(g.getZ())) { pair[0] = g.getX() / g.getZ(); pair[1] = g.getY() / g.getZ(); } } } setPoints(valPairs); }
From source file:geogebra.common.kernel.implicit.AlgoIntersectImplicitpolys.java
private static int getNearRoots(double[] roots, EquationSolverInterface solver, double epsilon) { PolynomialFunction poly = new PolynomialFunction(roots); double[] rootsDerivative = poly.polynomialDerivative().getCoefficients(); int nrRoots = getRoots(roots, solver); int nrDeRoots = getRoots(rootsDerivative, solver); for (int i = 0; i < nrDeRoots; i++) { if (Kernel.isEqual(poly.value(rootsDerivative[i]), 0, epsilon)) { if (nrRoots < roots.length) { roots[nrRoots++] = rootsDerivative[i]; }//from w w w . j ava 2 s .c o m } } if (nrRoots == 0) { //a wild guess, test if the root of the n-1 derivative is a root of the original poly as well //works in case of a polynomial with one root of really high multiplicity. double[] c = poly.getCoefficients(); int n = c.length - 1; if (n > 0) { double x = -c[n - 1] / n / c[n]; if (Kernel.isEqual(poly.value(x), 0)) { roots[0] = x; return 1; } } } if (nrRoots == 0) { PolynomialFunction derivative = poly.polynomialDerivative(); double x = 0; double err = Math.abs(poly.value(x)); double lastErr = err * 2; while (err < lastErr && err > Kernel.STANDARD_PRECISION) { double devVal = derivative.value(x); if (!Kernel.isZero(devVal)) x = x - poly.value(x) / devVal; else break; lastErr = err; err = Math.abs(poly.value(x)); } if (Kernel.isEqual(poly.value(x), 0, epsilon)) { roots[0] = x; return 1; } } Arrays.sort(roots, 0, nrRoots); return nrRoots; }