List of usage examples for org.apache.commons.math.analysis.polynomials PolynomialFunction PolynomialFunction
public PolynomialFunction(double c[])
From source file:geogebra.kernel.GeoImplicitPoly.java
/** * // w w w . j ava 2 s.c o m * @param p GenPolynomial * @param varN connected to p.ring.nvar; if varN[i] set, p univariat in x_i is fine * @return p univariat in x_i, with varN[i]==true converted to PolynomialFunction, else null<br> * varN[i] remains true, all other set to false */ public static PolynomialFunction getUnivariatPoly(GenPolynomial<BigRational> p, boolean[] varN) { ExpVector e = p.degreeVector(); double[] coeff; if (varN.length != p.ring.nvar) { Application.debug("Length of varN doesn't match ring.nvar"); return null; } int k = -1; for (int i = 0; i < varN.length; i++) { if (e.getVal(i) != 0) { if (k >= 0) { // Application.debug("not univariat: "+p); return null; } else { k = i; } } } for (int i = 0; i < varN.length && k < 0; i++) //if k<0 => constant => 'univariat' in every var, take first one allowed if (varN[i]) k = i; if (k < 0) { return null; } // Application.debug("univar-poly = "+p); if (!varN[k]) { // Application.debug("univariat in wrong var: "+p); return null; } for (int i = 0; i < varN.length; i++) { varN[i] = i == k; } ExpVector dir = ExpVector.create(varN.length, k, 1); coeff = new double[(int) e.getVal(k) + 1]; for (int i = coeff.length - 1; i >= 0; i--) { BigRational b = p.coefficient(e); coeff[i] = b.numerator().doubleValue() / b.denominator().doubleValue(); e = e.subtract(dir); } return new PolynomialFunction(coeff); }
From source file:org.geogebra.common.kernel.implicit.AlgoIntersectImplicitpolys.java
@Override public void compute() { if (c1 != null) { p2 = new GeoImplicitPoly(c1); }//w w w . j a v a 2s . c o 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: Intersect[-5 x^4+ x^2+ y^2 = 0, -20 x^3+2 x^2+2 x+2 y^2+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)); //Gauss-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); }