Example usage for org.apache.commons.math.analysis.polynomials PolynomialFunction getCoefficients

List of usage examples for org.apache.commons.math.analysis.polynomials PolynomialFunction getCoefficients

Introduction

In this page you can find the example usage for org.apache.commons.math.analysis.polynomials PolynomialFunction getCoefficients.

Prototype

public double[] getCoefficients() 

Source Link

Document

Returns a copy of the coefficients array.

Usage

From source file:geogebra.kernel.implicit.PolynomialUtils.java

public static int getDegree(PolynomialFunction p) {
    return getDegree(p.getCoefficients());
}

From source file:geogebra.kernel.implicit.PolynomialUtils.java

public static double getLeadingCoeff(PolynomialFunction p) {
    return getDegree(p.getCoefficients());
}

From source file:geogebra.common.kernel.implicit.PolynomialUtils.java

/**
 * @param p polynomial function/*from   ww  w .j  a  va  2  s  .c o m*/
 * @return degree of the function
 */
public static int getDegree(PolynomialFunction p) {
    return getDegree(p.getCoefficients());
}

From source file:geogebra.common.kernel.implicit.PolynomialUtils.java

/**
 * @param p polynomial function//from w ww.  j a  v a 2 s.  c o m
 * @return leading coefficient
 */
public static double getLeadingCoeff(PolynomialFunction p) {
    return getLeadingCoeff(p.getCoefficients());
}

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/* ww w .j  av  a2s . c  o m*/
 * @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.common.kernel.algos.AlgoSimpleRootsPolynomial.java

protected void doCalc(PolynomialFunction rootsPoly) {
    double roots[] = rootsPoly.getCoefficients();
    int nrRealRoots = 0;
    if (roots.length > 1)
        nrRealRoots = getRoots(roots, eqnSolver);
    makePoints(roots, nrRealRoots);//from  w w w .  j a  va 2  s . c  o m
}

From source file:geogebra.kernel.AlgoSimpleRootsPolynomial.java

protected void doCalc(PolynomialFunction rootsPoly, double min, double max) {
    double roots[] = rootsPoly.getCoefficients();
    int nrRealRoots = 0;
    if (roots.length > 1)
        nrRealRoots = getRoots(roots, eqnSolver);

    for (int i = 0; i < nrRealRoots; ++i) {
        if (Kernel.isGreater(roots[i], max, kernel.getEpsilon())
                || Kernel.isGreater(min, roots[i], kernel.getEpsilon()))
            roots[i] = Double.NaN;
    }//ww  w .ja  v  a 2 s  .c  o m
    makePoints(roots, nrRealRoots);
}

From source file:geogebra.common.kernel.algos.AlgoSimpleRootsPolynomial.java

protected void doCalc(PolynomialFunction rootsPoly, double min, double max) {
    double roots[] = rootsPoly.getCoefficients();
    int nrRealRoots = 0;
    if (roots.length > 1)
        nrRealRoots = getRoots(roots, eqnSolver);

    for (int i = 0; i < nrRealRoots; ++i) {
        if (Kernel.isGreater(roots[i], max, Kernel.STANDARD_PRECISION)
                || Kernel.isGreater(min, roots[i], Kernel.STANDARD_PRECISION))
            roots[i] = Double.NaN;
    }//  w ww.j  av  a  2 s  .  c o  m
    makePoints(roots, nrRealRoots);
}

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  a  va2 s . c  om*/
        }
    }
    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;
}

From source file:geogebra.kernel.implicit.AlgoIntersectImplicitpolys.java

@Override
protected void compute() {
    if (c1 != null) {
        p2 = new GeoImplicitPoly(c1);
    }//from  w  w w . j a  v  a  2s.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;
         */
}