removed the constraint on low degree polynomials

when building Chebyshev, Hermite, Laguerre or Legendre polynomials

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@760901 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-04-01 14:29:18 +00:00
parent 3664a2cc7f
commit f5618f7717
3 changed files with 66 additions and 76 deletions

View File

@ -18,7 +18,7 @@ package org.apache.commons.math.analysis.polynomials;
import java.util.ArrayList;
import org.apache.commons.math.fraction.Fraction;
import org.apache.commons.math.fraction.BigFraction;
/**
* A collection of static methods that operate on or return polynomials.
@ -29,46 +29,46 @@ import org.apache.commons.math.fraction.Fraction;
public class PolynomialsUtils {
/** Coefficients for Chebyshev polynomials. */
private static final ArrayList<Fraction> CHEBYSHEV_COEFFICIENTS;
private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
/** Coefficients for Hermite polynomials. */
private static final ArrayList<Fraction> HERMITE_COEFFICIENTS;
private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
/** Coefficients for Laguerre polynomials. */
private static final ArrayList<Fraction> LAGUERRE_COEFFICIENTS;
private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
/** Coefficients for Legendre polynomials. */
private static final ArrayList<Fraction> LEGENDRE_COEFFICIENTS;
private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
static {
// initialize recurrence for Chebyshev polynomials
// T0(X) = 1, T1(X) = 0 + 1 * X
CHEBYSHEV_COEFFICIENTS = new ArrayList<Fraction>();
CHEBYSHEV_COEFFICIENTS.add(Fraction.ONE);
CHEBYSHEV_COEFFICIENTS.add(Fraction.ZERO);
CHEBYSHEV_COEFFICIENTS.add(Fraction.ONE);
CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
// initialize recurrence for Hermite polynomials
// H0(X) = 1, H1(X) = 0 + 2 * X
HERMITE_COEFFICIENTS = new ArrayList<Fraction>();
HERMITE_COEFFICIENTS.add(Fraction.ONE);
HERMITE_COEFFICIENTS.add(Fraction.ZERO);
HERMITE_COEFFICIENTS.add(Fraction.TWO);
HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
HERMITE_COEFFICIENTS.add(BigFraction.ONE);
HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
HERMITE_COEFFICIENTS.add(BigFraction.TWO);
// initialize recurrence for Laguerre polynomials
// L0(X) = 1, L1(X) = 1 - 1 * X
LAGUERRE_COEFFICIENTS = new ArrayList<Fraction>();
LAGUERRE_COEFFICIENTS.add(Fraction.ONE);
LAGUERRE_COEFFICIENTS.add(Fraction.ONE);
LAGUERRE_COEFFICIENTS.add(Fraction.MINUS_ONE);
LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
// initialize recurrence for Legendre polynomials
// P0(X) = 1, P1(X) = 0 + 1 * X
LEGENDRE_COEFFICIENTS = new ArrayList<Fraction>();
LEGENDRE_COEFFICIENTS.add(Fraction.ONE);
LEGENDRE_COEFFICIENTS.add(Fraction.ZERO);
LEGENDRE_COEFFICIENTS.add(Fraction.ONE);
LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
}
@ -94,9 +94,9 @@ public class PolynomialsUtils {
public static PolynomialFunction createChebyshevPolynomial(final int degree) {
return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
private final Fraction[] coeffs = { Fraction.ZERO, Fraction.TWO, Fraction.ONE};
private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
/** {@inheritDoc} */
public Fraction[] generate(int k) {
public BigFraction[] generate(int k) {
return coeffs;
}
});
@ -120,11 +120,11 @@ public class PolynomialsUtils {
return buildPolynomial(degree, HERMITE_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
/** {@inheritDoc} */
public Fraction[] generate(int k) {
return new Fraction[] {
Fraction.ZERO,
Fraction.TWO,
new Fraction(2 * k, 1)};
public BigFraction[] generate(int k) {
return new BigFraction[] {
BigFraction.ZERO,
BigFraction.TWO,
new BigFraction(2 * k)};
}
});
}
@ -146,12 +146,12 @@ public class PolynomialsUtils {
return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
/** {@inheritDoc} */
public Fraction[] generate(int k) {
public BigFraction[] generate(int k) {
final int kP1 = k + 1;
return new Fraction[] {
new Fraction(2 * k + 1, kP1),
new Fraction(-1, kP1),
new Fraction(k, kP1)};
return new BigFraction[] {
new BigFraction(2 * k + 1, kP1),
new BigFraction(-1, kP1),
new BigFraction(k, kP1)};
}
});
}
@ -173,12 +173,12 @@ public class PolynomialsUtils {
return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
/** {@inheritDoc} */
public Fraction[] generate(int k) {
public BigFraction[] generate(int k) {
final int kP1 = k + 1;
return new Fraction[] {
Fraction.ZERO,
new Fraction(k + kP1, kP1),
new Fraction(k, kP1)};
return new BigFraction[] {
BigFraction.ZERO,
new BigFraction(k + kP1, kP1),
new BigFraction(k, kP1)};
}
});
}
@ -190,7 +190,7 @@ public class PolynomialsUtils {
* @return coefficients array
*/
private static PolynomialFunction buildPolynomial(final int degree,
final ArrayList<Fraction> coefficients,
final ArrayList<BigFraction> coefficients,
final RecurrenceCoefficientsGenerator generator) {
final int maxDegree = (int) Math.floor(Math.sqrt(2 * coefficients.size())) - 1;
@ -228,7 +228,7 @@ public class PolynomialsUtils {
*/
private static void computeUpToDegree(final int degree, final int maxDegree,
final RecurrenceCoefficientsGenerator generator,
final ArrayList<Fraction> coefficients) {
final ArrayList<BigFraction> coefficients) {
int startK = (maxDegree - 1) * maxDegree / 2;
for (int k = maxDegree; k < degree; ++k) {
@ -238,24 +238,24 @@ public class PolynomialsUtils {
startK += k;
// Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
Fraction[] ai = generator.generate(k);
BigFraction[] ai = generator.generate(k);
Fraction ck = coefficients.get(startK);
Fraction ckm1 = coefficients.get(startKm1);
BigFraction ck = coefficients.get(startK);
BigFraction ckm1 = coefficients.get(startKm1);
// degree 0 coefficient
coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
// degree 1 to degree k-1 coefficients
for (int i = 1; i < k; ++i) {
final Fraction ckPrev = ck;
final BigFraction ckPrev = ck;
ck = coefficients.get(startK + i);
ckm1 = coefficients.get(startKm1 + i);
coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
}
// degree k coefficient
final Fraction ckPrev = ck;
final BigFraction ckPrev = ck;
ck = coefficients.get(startK + k);
coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
@ -274,7 +274,7 @@ public class PolynomialsUtils {
* @return an array of three coefficients such that
* P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
*/
Fraction[] generate(int k);
BigFraction[] generate(int k);
}
}

View File

@ -324,9 +324,8 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
one, using traditional coefficients arrays. The <a
href="../apidocs/org/apache/commons/math/analysis/polynomials/PolynomialsUtils.html">
org.apache.commons.math.analysis.polynomials.PolynomialsUtils</a> utility class provides static
factory methods to build Chebyshev, Hermite, Lagrange and Legendre polynomials. Beware that due
to overflows in the coefficients computations, these factory methods can only build low degrees
polynomials yet.
factory methods to build Chebyshev, Hermite, Lagrange and Legendre polynomials. Coefficients
are computed using exact fractions so these factory methods can build polynomials up to any degree.
</p>
</subsection>
</section>

View File

@ -177,34 +177,25 @@ public class PolynomialsUtilsTest extends TestCase {
}
public void testHighDegreeLegendre() {
try {
PolynomialsUtils.createLegendrePolynomial(40);
fail("an exception should have been thrown");
} catch (ArithmeticException ae) {
// expected
PolynomialsUtils.createLegendrePolynomial(40);
double[] l40 = PolynomialsUtils.createLegendrePolynomial(40).getCoefficients();
double denominator = 274877906944.0;
double[] numerators = new double[] {
+34461632205.0, -28258538408100.0, +3847870979902950.0, -207785032914759300.0,
+5929294332103310025.0, -103301483474866556880.0, +1197358103913226000200.0, -9763073770369381232400.0,
+58171647881784229843050.0, -260061484647976556945400.0, +888315281771246239250340.0, -2345767627188139419665400.0,
+4819022625419112503443050.0, -7710436200670580005508880.0, +9566652323054238154983240.0, -9104813935044723209570256.0,
+6516550296251767619752905.0, -3391858621221953912598660.0, +1211378079007840683070950.0, -265365894974690562152100.0,
+26876802183334044115405.0
};
for (int i = 0; i < l40.length; ++i) {
if (i % 2 == 0) {
double ci = numerators[i / 2] / denominator;
assertEquals(ci, l40[i], ci * 1.0e-15);
} else {
assertEquals(0.0, l40[i], 0.0);
}
}
// checkPolynomial(PolynomialsUtils.createLegendrePolynomial(40), 274877906944l,
// "34461632205.0"
// + " - 28258538408100.0 x^2"
// + " + 3847870979902950.0 x^4"
// + " - 207785032914759300.0 x^6"
// + " + 5929294332103310025.0 x^8"
// + " - 103301483474866556880.0 x^10"
// + " + 1197358103913226000200.0 x^12"
// + " - 9763073770369381232400.0 x^14"
// + " + 58171647881784229843050.0 x^16"
// + " - 260061484647976556945400.0 x^18"
// + " + 888315281771246239250340.0 x^20"
// + " - 2345767627188139419665400.0 x^22"
// + " + 4819022625419112503443050.0 x^24"
// + " - 7710436200670580005508880.0 x^26"
// + " + 9566652323054238154983240.0 x^28"
// + " - 9104813935044723209570256.0 x^30"
// + " + 6516550296251767619752905.0 x^32"
// + " - 3391858621221953912598660.0 x^34"
// + " + 1211378079007840683070950.0 x^36"
// + " - 265365894974690562152100.0 x^38"
// + " + 26876802183334044115405.0 x^40");
}
private void checkPolynomial(PolynomialFunction p, long denominator, String reference) {