org.ddogleg.solver.PolynomialOps Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ddogleg Show documentation
Show all versions of ddogleg Show documentation
DDogleg Numerics is a high performance Java library for non-linear optimization, robust model fitting, polynomial root finding, sorting, and more.
/*
* Copyright (c) 2012-2015, Peter Abeles. All Rights Reserved.
*
* This file is part of DDogleg (http://ddogleg.org).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.ddogleg.solver;
import org.ddogleg.solver.impl.FindRealRootsSturm;
import org.ddogleg.solver.impl.RootFinderCompanion;
import org.ddogleg.solver.impl.SturmSequence;
import org.ddogleg.solver.impl.WrapRealRootsSturm;
/**
* @author Peter Abeles
*/
public class PolynomialOps {
/**
* Given the coefficients compute the vertex (minimum/maximum) of the quadratic.
*
* y = a*x2 + b*x + c
*
* @param a quadratic coefficient.
* @param b quadratic coefficient.
* @return The quadratic's vertex.
*/
public static double quadraticVertex( double a, double b ) {
return -b/(2.0*a);
}
public static Polynomial derivative( Polynomial poly , Polynomial deriv ) {
if( deriv == null ) {
deriv = new Polynomial(poly.size-1);
} else {
deriv.size = poly.size - 1;
}
for( int i = 1; i < poly.size; i++ ) {
deriv.c[i-1] = poly.c[i]*i;
}
return deriv;
}
// TODO try using a linear search alg here
public static double refineRoot( Polynomial poly, double root , int maxIterations ) {
// for( int i = 0; i < maxIterations; i++ ) {
//
// double v = poly.c[poly.size-1];
// double d = v*(poly.size-1);
//
// for( int j = poly.size-1; j > 0; j-- ) {
// v = poly.c[j] + v*root;
// d = poly.c[j]*j + d*root;
// }
// v = poly.c[0] + v*root;
//
// if( d == 0 )
// return root;
//
// root -= v/d;
// }
//
// return root;
Polynomial deriv = new Polynomial(poly.size());
derivative(poly,deriv);
for( int i = 0; i < maxIterations; i++ ) {
double v = poly.evaluate(root);
double d = deriv.evaluate(root);
if( d == 0 )
return root;
root -= v/d;
}
return root;
}
/**
*
* Polynomial division. Computes both the quotient and the remainder.
*
* quotient = numerator/denominator
* remainder = numerator % denominator
*
*
* @param numerator Numerator in the division. Not modified.
* @param denominator Denominator in the division. Not modified.
* @param quotient Output quotient, Modified.
* @param remainder Output remainder. Modified.
*/
public static void divide( Polynomial numerator , Polynomial denominator , Polynomial quotient , Polynomial remainder ) {
if( denominator.size <= 0 )
throw new IllegalArgumentException("Trying to device by a polynomial of size 0");
int nn = numerator.size-1; int nd = denominator.size-1;
while( nd >= 0 && denominator.c[nd] == 0 )
nd -= 1;
// divide by zero error
if( nd < 0 ) {
throw new IllegalArgumentException("Divide by zero error");
}
quotient.size = nn-nd+1;
remainder.setTo(numerator);
for( int k = nn-nd; k >= 0; k-- ) {
if( nd < 0 || k < 0 )
System.out.println("EGADS");
double c = quotient.c[k] = remainder.c[nd+k]/denominator.c[nd];
for( int j = k+nd; j >= k; j-- ) {
remainder.c[j] -= c*denominator.c[j-k];
}
}
// The remainder can't be larger than the denominator
remainder.size = nd;
}
/**
* Multiplies the two polynomials together.
* @param a Polynomial
* @param b Polynomial
* @param result Optional storage parameter for the results. Must be have enough coefficients to store the results.
* If null a new instance is declared.
* @return Results of the multiplication
*/
public static Polynomial multiply( Polynomial a , Polynomial b , Polynomial result ) {
int N = Math.max(0,a.size() + b.size() - 1);
if( result == null ) {
result = new Polynomial(N);
} else {
if( result.size < N )
throw new IllegalArgumentException("Unexpected length of 'result'");
result.zero();
}
for( int i = 0; i < a.size; i++ ) {
double coef = a.c[i];
int index = i;
for( int j = 0; j < b.size; j++ ) {
result.c[index++] += coef*b.c[j];
}
}
return result;
}
/**
* Adds two polynomials together. The lengths of the polynomials do not need to be the zero, but the 'missing'
* coefficients are assumed to be zero.
*
* @param a Polynomial
* @param b Polynomial
* @param results Optional storage for resulting polynomial. If null a new instance is declared. If not null
* its length must be the same as the largest polynomial 'a' or 'b'.
* @return Polynomial 'a' and 'b' added together.
*/
public static Polynomial add( Polynomial a , Polynomial b , Polynomial results ) {
int N = Math.max(a.size,b.size);
if( results == null ) {
results = new Polynomial(N);
} else if( results.size < N ) {
throw new IllegalArgumentException("storage for results must be at least as large as the the largest polynomial");
} else {
for( int i = N; i < results.size; i++ ) {
results.c[i] = 0;
}
}
int M = Math.min(a.size,b.size);
for( int i = 0; i < M; i++ ) {
results.c[i] = a.c[i] + b.c[i];
}
// copy the non-overlapping coefficients for the larger polynomial
if( a.size > b.size ) {
for( int i = b.size; i < N; i++ )
results.c[i] = a.c[i];
} else {
for( int i = a.size; i < N; i++ )
results.c[i] = b.c[i];
}
return results;
}
public static int countRealRoots( Polynomial poly ) {
SturmSequence sturm = new SturmSequence(poly.size);
sturm.initialize(poly);
return sturm.countRealRoots(Double.NEGATIVE_INFINITY,Double.POSITIVE_INFINITY);
}
/**
* Creates different polynomial root finders.
*
* @param maxCoefficients The maximum number of coefficients that will be processed. This is the order + 1
* @param which 0 = Sturm and 1 = companion matrix.
* @return PolynomialRoots
*/
public static PolynomialRoots createRootFinder( int maxCoefficients , RootFinderType which ) {
switch( which ) {
case STURM:
FindRealRootsSturm sturm = new FindRealRootsSturm(maxCoefficients,-1,1e-10,200,200);
return new WrapRealRootsSturm(sturm);
case EVD:
return new RootFinderCompanion();
default:
throw new IllegalArgumentException("Unknown algorithm: "+which);
}
}
/**
* Returns a real root to the cubic polynomial: 0 = c0 + c1*x + c2*x^2 + c3*c^3. There can be other
* real roots.
*
* WARNING: This technique is much less stable than using one of the RootFinder algorithms
*
* @param c0 Polynomial coefficient for power 0
* @param c1 Polynomial coefficient for power 1
* @param c2 Polynomial coefficient for power 2
* @param c3 Polynomial coefficient for power 3
* @return A real root of the polynomial
*/
public static double cubicRealRoot( double c0 , double c1 , double c2 , double c3 ) {
// convert it into this format: r + q*x + p*x^2 + x^3 = 0
double p = c2/c3;
double q = c1/c3;
double r = c0/c3;
// reduce by substitution x = y - p/3 which results in y^3 + a*y + b = 0
double a = (3*q - p*p)/3.0;
double b = (2*p*p*p - 9*p*q + 27*r)/27.0;
double left = -b/2.0;
double right = Math.sqrt(b*b/4.0 + a*a*a/27.0);
double inner1 = left+right;
double inner2 = left-right;
double A,B;
if( inner1 < 0 )
A = -Math.pow(-inner1,1.0/3.0);
else
A = Math.pow(inner1,1.0/3.0);
if( inner2 < 0 )
B = -Math.pow(-inner2,1.0/3.0);
else
B = Math.pow(inner2,1.0/3.0);
return (A + B) - p/3.0;
}
}