added a PolynomialsUtils class providing factory methods for

Chebyshev, Hermite, Laguerre and Legendre polynomials
the code was extracted from mantissa and modified

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@739840 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-02-01 21:13:55 +00:00
parent 6a965532e6
commit 8ce2128585
16 changed files with 529 additions and 839 deletions

View File

@ -0,0 +1,280 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You 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.apache.commons.math.analysis.polynomials;
import java.util.ArrayList;
import org.apache.commons.math.fraction.Fraction;
/**
* A collection of static methods that operate on or return polynomials.
*
* @version $Revision$ $Date$
* @since 2.0
*/
public class PolynomialsUtils {
/** Coefficients for Chebyshev polynomials. */
private static final ArrayList<Fraction> CHEBYSHEV_COEFFICIENTS;
/** Coefficients for Hermite polynomials. */
private static final ArrayList<Fraction> HERMITE_COEFFICIENTS;
/** Coefficients for Laguerre polynomials. */
private static final ArrayList<Fraction> LAGUERRE_COEFFICIENTS;
/** Coefficients for Legendre polynomials. */
private static final ArrayList<Fraction> 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);
// 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);
// 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);
// 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);
}
/**
* Private constructor, to prevent instantiation.
*/
private PolynomialsUtils() {
}
/**
* Create a Chebyshev polynomial of the first kind.
* <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
* polynomials of the first kind</a> are orthogonal polynomials.
* They can be defined by the following recurrence relations:
* <pre>
* T<sub>0</sub>(X) = 1
* T<sub>1</sub>(X) = X
* T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
* </pre></p>
* @param degree degree of the polynomial
* @return Chebyshev polynomial of specified degree
*/
public static PolynomialFunction createChebyshevPolynomial(final int degree) {
return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
private final Fraction[] coeffs = { Fraction.ZERO, Fraction.TWO, Fraction.ONE};
/** {@inheritDoc} */
public Fraction[] generate(int k) {
return coeffs;
}
});
}
/**
* Create a Hermite polynomial.
* <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
* polynomials</a> are orthogonal polynomials.
* They can be defined by the following recurrence relations:
* <pre>
* H<sub>0</sub>(X) = 1
* H<sub>1</sub>(X) = 2X
* H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
* </pre></p>
* @param degree degree of the polynomial
* @return Hermite polynomial of specified degree
*/
public static PolynomialFunction createHermitePolynomial(final int degree) {
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)};
}
});
}
/**
* Create a Laguerre polynomial.
* <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
* polynomials</a> are orthogonal polynomials.
* They can be defined by the following recurrence relations:
* <pre>
* L<sub>0</sub>(X) = 1
* L<sub>1</sub>(X) = 1 - X
* (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
* </pre></p>
* @param degree degree of the polynomial
* @return Laguerre polynomial of specified degree
*/
public static PolynomialFunction createLaguerrePolynomial(final int degree) {
return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
/** {@inheritDoc} */
public Fraction[] 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)};
}
});
}
/**
* Create a Legendre polynomial.
* <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
* polynomials</a> are orthogonal polynomials.
* They can be defined by the following recurrence relations:
* <pre>
* P<sub>0</sub>(X) = 1
* P<sub>1</sub>(X) = X
* (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
* </pre></p>
* @param degree degree of the polynomial
* @return Legendre polynomial of specified degree
*/
public static PolynomialFunction createLegendrePolynomial(final int degree) {
return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
new RecurrenceCoefficientsGenerator() {
/** {@inheritDoc} */
public Fraction[] generate(int k) {
final int kP1 = k + 1;
return new Fraction[] {
Fraction.ZERO,
new Fraction(k + kP1, kP1),
new Fraction(k, kP1)};
}
});
}
/** Get the coefficients array for a given degree.
* @param degree degree of the polynomial
* @param coefficients list where the computed coefficients are stored
* @param generator recurrence coefficients generator
* @return coefficients array
*/
private static PolynomialFunction buildPolynomial(final int degree,
final ArrayList<Fraction> coefficients,
final RecurrenceCoefficientsGenerator generator) {
final int maxDegree = (int) Math.floor(Math.sqrt(2 * coefficients.size())) - 1;
synchronized (PolynomialsUtils.class) {
if (degree > maxDegree) {
computeUpToDegree(degree, maxDegree, generator, coefficients);
}
}
// coefficient for polynomial 0 is l [0]
// coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
// coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
// coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
// coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
// coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
// coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
// ...
final int start = degree * (degree + 1) / 2;
final double[] a = new double[degree + 1];
for (int i = 0; i <= degree; ++i) {
a[i] = coefficients.get(start + i).doubleValue();
}
// build the polynomial
return new PolynomialFunction(a);
}
/** Compute polynomial coefficients up to a given degree.
* @param degree maximal degree
* @param maxDegree current maximal degree
* @param generator recurrence coefficients generator
* @param coefficients list where the computed coefficients should be appended
*/
private static void computeUpToDegree(final int degree, final int maxDegree,
final RecurrenceCoefficientsGenerator generator,
final ArrayList<Fraction> coefficients) {
int startK = (maxDegree - 1) * maxDegree / 2;
for (int k = maxDegree; k < degree; ++k) {
// start indices of two previous polynomials Pk(X) and Pk-1(X)
int startKm1 = startK;
startK += k;
// Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
Fraction[] ai = generator.generate(k);
Fraction ck = coefficients.get(startK);
Fraction 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;
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;
ck = coefficients.get(startK + k);
coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
// degree k+1 coefficient
coefficients.add(ck.multiply(ai[1]));
}
}
/** Interface for recurrence coefficients generation. */
private static interface RecurrenceCoefficientsGenerator {
/**
* Generate recurrence coefficients.
* @param k highest degree of the polynomials used in the recurrence
* @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);
}
}

View File

@ -29,14 +29,20 @@ import org.apache.commons.math.util.MathUtils;
*/
public class Fraction extends Number implements Comparable<Fraction> {
/** A fraction representing "2 / 1". */
public static final Fraction TWO = new Fraction(2, 1);
/** A fraction representing "1 / 1". */
public static final Fraction ONE = new Fraction(1, 1);
/** A fraction representing "0 / 1". */
public static final Fraction ZERO = new Fraction(0, 1);
/** A fraction representing "-1 / 1". */
public static final Fraction MINUS_ONE = new Fraction(-1, 1);
/** Serializable version identifier */
private static final long serialVersionUID = -5731055832688548463L;
private static final long serialVersionUID = 3071409609509774764L;
/** The denominator. */
private final int denominator;
@ -197,7 +203,7 @@ public class Fraction extends Number implements Comparable<Fraction> {
* reduced to lowest terms.
* @param num the numerator.
* @param den the denominator.
* @throws ArithmeticException if the denomiator is <code>zero</code>
* @throws ArithmeticException if the denominator is <code>zero</code>
*/
public Fraction(int num, int den) {
super();

View File

@ -1,70 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
/**
* This class implements Chebyshev polynomials.
* <p>Chebyshev polynomials can be defined by the following recurrence
* relations:
* <pre>
* T<sub>0</sub>(X) = 1
* T<sub>1</sub>(X) = X
* T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
* </pre></p>
* @version $Id$
* @author L. Maisonobe
*/
public class Chebyshev
extends OrthogonalPolynomial {
/** Generator for the Chebyshev polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(0l),
new RationalNumber(1l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// Tk+1(X) = 2X Tk(X) - Tk-1(X)
setRecurrenceCoefficients(new RationalNumber(0l),
new RationalNumber(2l),
new RationalNumber(1l));
}
};
/** Simple constructor.
* Build a degree 0 Chebyshev polynomial
*/
public Chebyshev() {
super(0, generator);
}
/** Simple constructor.
* Build a degree d Chebyshev polynomial
* @param degree degree of the polynomial
*/
public Chebyshev(int degree) {
super(degree, generator);
}
private static final long serialVersionUID = -893367988717182601L;
}

View File

@ -1,157 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
import java.util.ArrayList;
public abstract class CoefficientsGenerator {
/** Build a generator with coefficients for two polynomials.
* <p>The first polynomial must be a degree 0 polynomial
* P<sub>0</sub>(X)=a<sub>0,0</sub> and the second polynomial
* must be a degree 1 polynomial P<sub>1</sub>(X)=a<sub>0,1</sub>
* +a<sub>1,1</sub>X</p>
* @param a00 constant term for the degree 0 polynomial
* @param a01 constant term for the degree 1 polynomial
* @param a11 X term for the degree 1 polynomial
*/
protected CoefficientsGenerator(RationalNumber a00,
RationalNumber a01, RationalNumber a11) {
l = new ArrayList();
l.add(a00);
l.add(a01);
l.add(a11);
maxDegree = 1;
}
/** Set the recurrence coefficients.
* @param b2k b<sub>2,k</sub> coefficient (b<sub>2,k</sub> = a<sub>2,k</sub> / a<sub>1,k</sub>)
* @param b3k b<sub>3,k</sub> coefficient (b<sub>3,k</sub> = a<sub>3,k</sub> / a<sub>1,k</sub>)
* @param b4k b<sub>4,k</sub> coefficient (b<sub>4,k</sub> = a<sub>4,k</sub> / a<sub>1,k</sub>)
*/
protected void setRecurrenceCoefficients(RationalNumber b2k,
RationalNumber b3k,
RationalNumber b4k) {
this.b2k = b2k;
this.b3k = b3k;
this.b4k = b4k;
}
/** Set the recurrence coefficients.
* The recurrence relation is
* <pre>a<sub>1,k</sub> O<sub>k+1</sub>(X) =(a<sub>2,k</sub> + a<sub>3,k</sub> X) O<sub>k</sub>(X) - a<sub>4,k</sub> O<sub>k-1</sub>(X)</pre>
* the method must call {@link #setRecurrenceCoefficients(RationalNumber,
* RationalNumber, RationalNumber)} to provide the coefficients
* @param k index of the current step
*/
protected abstract void setRecurrenceCoefficients(int k);
/** Compute all the polynomial coefficients up to a given degree.
* @param degree maximal degree
*/
private void computeUpToDegree(int degree) {
int startK = (maxDegree - 1) * maxDegree / 2;
for (int k = maxDegree; k < degree; ++k) {
// start indices of two previous polynomials Ok(X) and Ok-1(X)
int startKm1 = startK;
startK += k;
// a1k Ok+1(X) = (a2k + a3k X) Ok(X) - a4k Ok-1(X)
// we use bik = aik/a1k
setRecurrenceCoefficients(k);
RationalNumber ckPrev = null;
RationalNumber ck = (RationalNumber) l.get(startK);
RationalNumber ckm1 = (RationalNumber) l.get(startKm1);
// degree 0 coefficient
l.add(ck.multiply(b2k).subtract(ckm1.multiply(b4k)));
// degree 1 to degree k-1 coefficients
for (int i = 1; i < k; ++i) {
ckPrev = ck;
ck = (RationalNumber) l.get(startK + i);
ckm1 = (RationalNumber) l.get(startKm1 + i);
l.add(ck.multiply(b2k).add(ckPrev.multiply(b3k)).subtract(ckm1.multiply(b4k)));
}
// degree k coefficient
ckPrev = ck;
ck = (RationalNumber) l.get(startK + k);
l.add(ck.multiply(b2k).add(ckPrev.multiply(b3k)));
// degree k+1 coefficient
l.add(ck.multiply(b3k));
}
maxDegree = degree;
}
/** Get the coefficients array for a given degree.
* @param degree degree of the polynomial
* @return coefficients array
*/
public RationalNumber[] getCoefficients(int degree) {
synchronized (this) {
if (degree > maxDegree) {
computeUpToDegree(degree);
}
}
// coefficient for polynomial 0 is l [0]
// coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
// coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
// coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
// coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
// coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
// coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
// ...
int start = degree * (degree + 1) / 2;
RationalNumber[] a = new RationalNumber[degree + 1];
for (int i = 0; i <= degree; ++i) {
a[i] = (RationalNumber) l.get(start + i);
}
return a;
}
/** List holding the coefficients of the polynomials computed so far. */
private ArrayList l;
/** Maximal degree of the polynomials computed so far. */
private int maxDegree;
/** b<sub>2,k</sub> coefficient to initialize
* (b<sub>2,k</sub> = a<sub>2,k</sub> / a<sub>1,k</sub>). */
private RationalNumber b2k;
/** b<sub>3,k</sub> coefficient to initialize
* (b<sub>3,k</sub> = a<sub>3,k</sub> / a<sub>1,k</sub>). */
private RationalNumber b3k;
/** b<sub>4,k</sub> coefficient to initialize
* (b<sub>4,k</sub> = a<sub>4,k</sub> / a<sub>1,k</sub>). */
private RationalNumber b4k;
}

View File

@ -1,69 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
/**
* This class implements Hermite polynomials.
* <p>Hermite polynomials can be defined by the following recurrence
* relations:
* <pre>
* H<sub>0</sub>(X) = 1
* H<sub>1</sub>(X) = 2X
* H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
* </pre></p>
* @version $Id$
* @author L. Maisonobe
*/
public class Hermite
extends OrthogonalPolynomial {
/** Generator for the Hermite polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(0l),
new RationalNumber(2l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// Hk+1(X) = 2X Hk(X) - 2k Hk-1(X)
setRecurrenceCoefficients(new RationalNumber(0l),
new RationalNumber(2l),
new RationalNumber(k * 2l));
}
};
/** Simple constructor.
* Build a degree 0 Hermite polynomial
*/
public Hermite() {
super(0, generator);
}
/** Simple constructor.
* Build a degree d Hermite polynomial
* @param degree degree of the polynomial
*/
public Hermite(int degree) {
super(degree, generator);
}
private static final long serialVersionUID = 7910082423686662133L;
}

View File

@ -1,70 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
/**
* This class implements Laguerre polynomials.
* <p>Laguerre polynomials can be defined by the following recurrence
* relations:
* <pre>
* L<sub>0</sub>(X) = 1
* L<sub>1</sub>(X) = 1 - X
* (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
* </pre></p>
* @version $Id$
* @author L. Maisonobe
*/
public class Laguerre
extends OrthogonalPolynomial {
/** Generator for the Laguerre polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(1l),
new RationalNumber(-1l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// (k+1) Lk+1(X) = (2k + 1 - X) Lk(X) - k Lk-1(X)
long kP1 = k + 1;
setRecurrenceCoefficients(new RationalNumber(2 * k + 1, kP1),
new RationalNumber(-1l, kP1),
new RationalNumber(k, kP1));
}
};
/** Simple constructor.
* Build a degree 0 Laguerre polynomial
*/
public Laguerre() {
super(0, generator);
}
/** Simple constructor.
* Build a degree d Laguerre polynomial
* @param degree degree of the polynomial
*/
public Laguerre(int degree) {
super(degree, generator);
}
private static final long serialVersionUID = 3213856667479179710L;
}

View File

@ -1,71 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
/**
* This class implements Legendre polynomials.
* <p>Legendre polynomials can be defined by the following recurrence
* relations:
* <pre>
* P<sub>0</sub>(X) = 1
* P<sub>1</sub>(X) = X
* (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
* </pre></p>
* @version $Id$
* @author L. Maisonobe
*/
public class Legendre
extends OrthogonalPolynomial {
/** Generator for the Legendre polynomials. */
private static final CoefficientsGenerator generator =
new CoefficientsGenerator(new RationalNumber(1l),
new RationalNumber(0l),
new RationalNumber(1l)) {
public void setRecurrenceCoefficients(int k) {
// the recurrence relation is
// (k+1) Pk+1(X) = (2k+1) X Pk(X) - k Pk-1(X)
long kP1 = k + 1;
setRecurrenceCoefficients(new RationalNumber(0l),
new RationalNumber(2 * k + 1, kP1),
new RationalNumber(k, kP1));
}
};
/** Simple constructor.
* Build a degree 0 Legendre polynomial
*/
public Legendre() {
super(0, generator);
}
/** Simple constructor.
* Build a degree d Legendre polynomial
* @param degree degree of the polynomial
*/
public Legendre(int degree) {
super(degree, generator);
}
private static final long serialVersionUID = 4014485393845978429L;
}

View File

@ -1,49 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
/**
* This class is the base class for orthogonal polynomials.
* <p>Orthogonal polynomials can be defined by recurrence relations like:
* <pre>
* O<sub>0</sub>(X) = some 0 degree polynomial
* O<sub>1</sub>(X) = some first degree polynomial
* a<sub>1,k</sub> O<sub>k+1</sub>(X) = (a<sub>2,k</sub> + a<sub>3,k</sub> X) O<sub>k</sub>(X) - a<sub>4,k</sub> O<sub>k-1</sub>(X)
* </pre>
* where a<sub>1,k</sub>, a<sub>2,k</sub>, a<sub>3,k</sub> and
* a<sub>4,k</sub> are simple expressions which either are
* constants or depend on k.</p>
* @version $Id$
* @author L. Maisonobe
*/
public abstract class OrthogonalPolynomial
extends Polynomial.Rational {
/** Simple constructor.
* Build a degree d orthogonal polynomial
* @param degree degree of the polynomial
* @param generator coefficients generator for the current type of polynomials
*/
protected OrthogonalPolynomial(int degree, CoefficientsGenerator generator) {
a = generator.getCoefficients(degree);
}
}

View File

@ -28,10 +28,6 @@ public class AllTests {
suite.addTest(RationalNumberTest.suite());
suite.addTest(PolynomialRationalTest.suite());
suite.addTest(PolynomialDoubleTest.suite());
suite.addTest(ChebyshevTest.suite());
suite.addTest(HermiteTest.suite());
suite.addTest(LegendreTest.suite());
suite.addTest(LaguerreTest.suite());
suite.addTest(PolynomialFractionTest.suite());
return suite;

View File

@ -1,85 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
import junit.framework.*;
public class ChebyshevTest
extends TestCase {
public ChebyshevTest(String name) {
super(name);
}
public void testOne() {
assertTrue(new Chebyshev().isOne());
}
public void testFirstPolynomials() {
checkPolynomial(new Chebyshev(3), "-3 x + 4 x^3");
checkPolynomial(new Chebyshev(2), "-1 + 2 x^2");
checkPolynomial(new Chebyshev(1), "x");
checkPolynomial(new Chebyshev(0), "1");
checkPolynomial(new Chebyshev(7), "-7 x + 56 x^3 - 112 x^5 + 64 x^7");
checkPolynomial(new Chebyshev(6), "-1 + 18 x^2 - 48 x^4 + 32 x^6");
checkPolynomial(new Chebyshev(5), "5 x - 20 x^3 + 16 x^5");
checkPolynomial(new Chebyshev(4), "1 - 8 x^2 + 8 x^4");
}
public void testBounds() {
for (int k = 0; k < 12; ++k) {
OrthogonalPolynomial Tk = new Chebyshev(k);
for (double x = -1.0; x <= 1.0; x += 0.02) {
assertTrue(Math.abs(Tk.valueAt(x)) < (1.0 + 1.0e-12));
}
}
}
public void testDifferentials() {
for (int k = 0; k < 12; ++k) {
Polynomial.Rational Tk0 = new Chebyshev(k);
Polynomial.Rational Tk1 = (Polynomial.Rational) Tk0.getDerivative();
Polynomial.Rational Tk2 = (Polynomial.Rational) Tk1.getDerivative();
Polynomial.Rational g0 = new Polynomial.Rational(k * k);
Polynomial.Rational g1 = new Polynomial.Rational(-1l, 0l);
Polynomial.Rational g2 = new Polynomial.Rational(-1l, 0l, 1l);
Polynomial.Rational Tk0g0 = Tk0.multiply(g0);
Polynomial.Rational Tk1g1 = Tk1.multiply(g1);
Polynomial.Rational Tk2g2 = Tk2.multiply(g2);
Polynomial.Rational d = Tk0g0.add(Tk1g1.add(Tk2g2));
assertTrue(d.isZero());
}
}
public void checkPolynomial(Polynomial.Rational p, String reference) {
assertTrue(p.toString().equals(reference));
}
public static Test suite() {
return new TestSuite(ChebyshevTest.class);
}
}

View File

@ -1,76 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
import junit.framework.*;
public class HermiteTest
extends TestCase {
public HermiteTest(String name) {
super(name);
}
public void testOne() {
assertTrue(new Hermite().isOne());
}
public void testFirstPolynomials() {
checkPolynomial(new Hermite(3), "-12 x + 8 x^3");
checkPolynomial(new Hermite(2), "-2 + 4 x^2");
checkPolynomial(new Hermite(1), "2 x");
checkPolynomial(new Hermite(0), "1");
checkPolynomial(new Hermite(7), "-1680 x + 3360 x^3 - 1344 x^5 + 128 x^7");
checkPolynomial(new Hermite(6), "-120 + 720 x^2 - 480 x^4 + 64 x^6");
checkPolynomial(new Hermite(5), "120 x - 160 x^3 + 32 x^5");
checkPolynomial(new Hermite(4), "12 - 48 x^2 + 16 x^4");
}
public void testDifferentials() {
for (int k = 0; k < 12; ++k) {
Polynomial.Rational Hk0 = new Hermite(k);
Polynomial.Rational Hk1 = (Polynomial.Rational) Hk0.getDerivative();
Polynomial.Rational Hk2 = (Polynomial.Rational) Hk1.getDerivative();
Polynomial.Rational g0 = new Polynomial.Rational(2l * k);
Polynomial.Rational g1 = new Polynomial.Rational(-2l, 0l);
Polynomial.Rational g2 = new Polynomial.Rational(1l);
Polynomial.Rational Hk0g0 = Hk0.multiply(g0);
Polynomial.Rational Hk1g1 = Hk1.multiply(g1);
Polynomial.Rational Hk2g2 = Hk2.multiply(g2);
Polynomial.Rational d = Hk0g0.add(Hk1g1.add(Hk2g2));
assertTrue(d.isZero());
}
}
public void checkPolynomial(Polynomial.Rational p, String reference) {
assertTrue(p.toString().equals(reference));
}
public static Test suite() {
return new TestSuite(HermiteTest.class);
}
}

View File

@ -1,82 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
import junit.framework.*;
public class LaguerreTest
extends TestCase {
public LaguerreTest(String name) {
super(name);
}
public void testOne() {
assertTrue(new Laguerre().isOne());
}
public void testFirstPolynomials() {
checkLaguerre(new Laguerre(3), 6l, "6 - 18 x + 9 x^2 - x^3");
checkLaguerre(new Laguerre(2), 2l, "2 - 4 x + x^2");
checkLaguerre(new Laguerre(1), 1l, "1 - x");
checkLaguerre(new Laguerre(0), 1l, "1");
checkLaguerre(new Laguerre(7), 5040l,
"5040 - 35280 x + 52920 x^2 - 29400 x^3"
+ " + 7350 x^4 - 882 x^5 + 49 x^6 - x^7");
checkLaguerre(new Laguerre(6), 720l,
"720 - 4320 x + 5400 x^2 - 2400 x^3 + 450 x^4"
+ " - 36 x^5 + x^6");
checkLaguerre(new Laguerre(5), 120l,
"120 - 600 x + 600 x^2 - 200 x^3 + 25 x^4 - x^5");
checkLaguerre(new Laguerre(4), 24l,
"24 - 96 x + 72 x^2 - 16 x^3 + x^4");
}
public void testDifferentials() {
for (int k = 0; k < 12; ++k) {
Polynomial.Rational Lk0 = new Laguerre(k);
Polynomial.Rational Lk1 = (Polynomial.Rational) Lk0.getDerivative();
Polynomial.Rational Lk2 = (Polynomial.Rational) Lk1.getDerivative();
Polynomial.Rational g0 = new Polynomial.Rational(k);
Polynomial.Rational g1 = new Polynomial.Rational(-1l, 1l);
Polynomial.Rational g2 = new Polynomial.Rational(1l, 0l);
Polynomial.Rational Lk0g0 = Lk0.multiply(g0);
Polynomial.Rational Lk1g1 = Lk1.multiply(g1);
Polynomial.Rational Lk2g2 = Lk2.multiply(g2);
Polynomial.Rational d = Lk0g0.add(Lk1g1.add(Lk2g2));
assertTrue(d.isZero());
}
}
public void checkLaguerre(Laguerre p, long denominator, String reference) {
assertTrue(p.multiply(denominator).toString().equals(reference));
}
public static Test suite() {
return new TestSuite(LaguerreTest.class);
}
}

View File

@ -1,101 +0,0 @@
// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements. See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership. The ASF licenses this file
// to you 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.spaceroots.mantissa.algebra;
import junit.framework.*;
public class LegendreTest
extends TestCase {
public LegendreTest(String name) {
super(name);
}
public void testOne() {
assertTrue(new Legendre().isOne());
}
public void testFirstPolynomials() {
checkLegendre(new Legendre(3), 2l, "-3 x + 5 x^3");
checkLegendre(new Legendre(2), 2l, "-1 + 3 x^2");
checkLegendre(new Legendre(1), 1l, "x");
checkLegendre(new Legendre(0), 1l, "1");
checkLegendre(new Legendre(7), 16l, "-35 x + 315 x^3 - 693 x^5 + 429 x^7");
checkLegendre(new Legendre(6), 16l, "-5 + 105 x^2 - 315 x^4 + 231 x^6");
checkLegendre(new Legendre(5), 8l, "15 x - 70 x^3 + 63 x^5");
checkLegendre(new Legendre(4), 8l, "3 - 30 x^2 + 35 x^4");
}
public void testDifferentials() {
for (int k = 0; k < 12; ++k) {
Polynomial.Rational Pk0 = new Legendre(k);
Polynomial.Rational Pk1 = (Polynomial.Rational) Pk0.getDerivative();
Polynomial.Rational Pk2 = (Polynomial.Rational) Pk1.getDerivative();
Polynomial.Rational g0 = new Polynomial.Rational(k * (k + 1));
Polynomial.Rational g1 = new Polynomial.Rational(-2l, 0l);
Polynomial.Rational g2 = new Polynomial.Rational(-1l, 0l, 1l);
Polynomial.Rational Pk0g0 = Pk0.multiply(g0);
Polynomial.Rational Pk1g1 = Pk1.multiply(g1);
Polynomial.Rational Pk2g2 = Pk2.multiply(g2);
Polynomial.Rational d = Pk0g0.add(Pk1g1.add(Pk2g2));
assertTrue(d.isZero());
}
}
public void testHighDegree() {
checkLegendre(new Legendre(40), 274877906944l,
"34461632205"
+ " - 28258538408100 x^2"
+ " + 3847870979902950 x^4"
+ " - 207785032914759300 x^6"
+ " + 5929294332103310025 x^8"
+ " - 103301483474866556880 x^10"
+ " + 1197358103913226000200 x^12"
+ " - 9763073770369381232400 x^14"
+ " + 58171647881784229843050 x^16"
+ " - 260061484647976556945400 x^18"
+ " + 888315281771246239250340 x^20"
+ " - 2345767627188139419665400 x^22"
+ " + 4819022625419112503443050 x^24"
+ " - 7710436200670580005508880 x^26"
+ " + 9566652323054238154983240 x^28"
+ " - 9104813935044723209570256 x^30"
+ " + 6516550296251767619752905 x^32"
+ " - 3391858621221953912598660 x^34"
+ " + 1211378079007840683070950 x^36"
+ " - 265365894974690562152100 x^38"
+ " + 26876802183334044115405 x^40");
}
public void checkLegendre(Legendre p, long denominator, String reference) {
assertTrue(p.multiply(denominator).toString().equals(reference));
}
public static Test suite() {
return new TestSuite(LegendreTest.class);
}
}

View File

@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="add" >
Added factory methods to create Chebyshev, Hermite, Laguerre and Legendre polynomials.
</action>
<action dev="luc" type="add" >
Added add, subtract, negate, multiply and toString methods to PolynomialFunction.
</action>

View File

@ -319,10 +319,20 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY);</source>
</subsection>
<subsection name="4.6 Polynomials" href="polynomials">
<p>
The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/package.html">
The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/package-summary.html">
org.apache.commons.math.analysis.polynomials</a> package provides real coefficients
polynomials.
</p>
<p>
The <a href="../apidocs/org/apache/commons/math/analysis/polynomials/PolynomialFunction.html">
org.apache.commons.math.analysis.polynomials.PolynomialFunction</a> class is the most general
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.
</p>
</subsection>
</section>
</body>

View File

@ -0,0 +1,225 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You 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.apache.commons.math.analysis.polynomials;
import junit.framework.TestCase;
/**
* Tests the PolynomialsUtils class.
*
* @version $Revision$ $Date$
*/
public class PolynomialsUtilsTest extends TestCase {
public void testFirstChebyshevPolynomials() {
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(3), "-3.0 x + 4.0 x^3");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(2), "-1.0 + 2.0 x^2");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(1), "x");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(0), "1.0");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(7), "-7.0 x + 56.0 x^3 - 112.0 x^5 + 64.0 x^7");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(6), "-1.0 + 18.0 x^2 - 48.0 x^4 + 32.0 x^6");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(5), "5.0 x - 20.0 x^3 + 16.0 x^5");
checkPolynomial(PolynomialsUtils.createChebyshevPolynomial(4), "1.0 - 8.0 x^2 + 8.0 x^4");
}
public void testChebyshevBounds() {
for (int k = 0; k < 12; ++k) {
PolynomialFunction Tk = PolynomialsUtils.createChebyshevPolynomial(k);
for (double x = -1.0; x <= 1.0; x += 0.02) {
assertTrue(k + " " + Tk.value(x), Math.abs(Tk.value(x)) < (1.0 + 1.0e-12));
}
}
}
public void testChebyshevDifferentials() {
for (int k = 0; k < 12; ++k) {
PolynomialFunction Tk0 = PolynomialsUtils.createChebyshevPolynomial(k);
PolynomialFunction Tk1 = Tk0.polynomialDerivative();
PolynomialFunction Tk2 = Tk1.polynomialDerivative();
PolynomialFunction g0 = new PolynomialFunction(new double[] { k * k });
PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -1});
PolynomialFunction g2 = new PolynomialFunction(new double[] { 1, 0, -1 });
PolynomialFunction Tk0g0 = Tk0.multiply(g0);
PolynomialFunction Tk1g1 = Tk1.multiply(g1);
PolynomialFunction Tk2g2 = Tk2.multiply(g2);
checkNullPolynomial(Tk0g0.add(Tk1g1.add(Tk2g2)));
}
}
public void testFirstHermitePolynomials() {
checkPolynomial(PolynomialsUtils.createHermitePolynomial(3), "-12.0 x + 8.0 x^3");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(2), "-2.0 + 4.0 x^2");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(1), "2.0 x");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(0), "1.0");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(7), "-1680.0 x + 3360.0 x^3 - 1344.0 x^5 + 128.0 x^7");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(6), "-120.0 + 720.0 x^2 - 480.0 x^4 + 64.0 x^6");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(5), "120.0 x - 160.0 x^3 + 32.0 x^5");
checkPolynomial(PolynomialsUtils.createHermitePolynomial(4), "12.0 - 48.0 x^2 + 16.0 x^4");
}
public void testHermiteDifferentials() {
for (int k = 0; k < 12; ++k) {
PolynomialFunction Hk0 = PolynomialsUtils.createHermitePolynomial(k);
PolynomialFunction Hk1 = Hk0.polynomialDerivative();
PolynomialFunction Hk2 = Hk1.polynomialDerivative();
PolynomialFunction g0 = new PolynomialFunction(new double[] { 2 * k });
PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -2 });
PolynomialFunction g2 = new PolynomialFunction(new double[] { 1 });
PolynomialFunction Hk0g0 = Hk0.multiply(g0);
PolynomialFunction Hk1g1 = Hk1.multiply(g1);
PolynomialFunction Hk2g2 = Hk2.multiply(g2);
checkNullPolynomial(Hk0g0.add(Hk1g1.add(Hk2g2)));
}
}
public void testFirstLaguerrePolynomials() {
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(3), 6l, "6.0 - 18.0 x + 9.0 x^2 - x^3");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(2), 2l, "2.0 - 4.0 x + x^2");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(1), 1l, "1.0 - x");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(0), 1l, "1.0");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(7), 5040l,
"5040.0 - 35280.0 x + 52920.0 x^2 - 29400.0 x^3"
+ " + 7350.0 x^4 - 882.0 x^5 + 49.0 x^6 - x^7");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(6), 720l,
"720.0 - 4320.0 x + 5400.0 x^2 - 2400.0 x^3 + 450.0 x^4"
+ " - 36.0 x^5 + x^6");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(5), 120l,
"120.0 - 600.0 x + 600.0 x^2 - 200.0 x^3 + 25.0 x^4 - x^5");
checkPolynomial(PolynomialsUtils.createLaguerrePolynomial(4), 24l,
"24.0 - 96.0 x + 72.0 x^2 - 16.0 x^3 + x^4");
}
public void testLaguerreDifferentials() {
for (int k = 0; k < 12; ++k) {
PolynomialFunction Lk0 = PolynomialsUtils.createLaguerrePolynomial(k);
PolynomialFunction Lk1 = Lk0.polynomialDerivative();
PolynomialFunction Lk2 = Lk1.polynomialDerivative();
PolynomialFunction g0 = new PolynomialFunction(new double[] { k });
PolynomialFunction g1 = new PolynomialFunction(new double[] { 1, -1 });
PolynomialFunction g2 = new PolynomialFunction(new double[] { 0, 1 });
PolynomialFunction Lk0g0 = Lk0.multiply(g0);
PolynomialFunction Lk1g1 = Lk1.multiply(g1);
PolynomialFunction Lk2g2 = Lk2.multiply(g2);
checkNullPolynomial(Lk0g0.add(Lk1g1.add(Lk2g2)));
}
}
public void testFirstLegendrePolynomials() {
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(3), 2l, "-3.0 x + 5.0 x^3");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(2), 2l, "-1.0 + 3.0 x^2");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(1), 1l, "x");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(0), 1l, "1.0");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(7), 16l, "-35.0 x + 315.0 x^3 - 693.0 x^5 + 429.0 x^7");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(6), 16l, "-5.0 + 105.0 x^2 - 315.0 x^4 + 231.0 x^6");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(5), 8l, "15.0 x - 70.0 x^3 + 63.0 x^5");
checkPolynomial(PolynomialsUtils.createLegendrePolynomial(4), 8l, "3.0 - 30.0 x^2 + 35.0 x^4");
}
public void testLegendreDifferentials() {
for (int k = 0; k < 12; ++k) {
PolynomialFunction Pk0 = PolynomialsUtils.createLegendrePolynomial(k);
PolynomialFunction Pk1 = Pk0.polynomialDerivative();
PolynomialFunction Pk2 = Pk1.polynomialDerivative();
PolynomialFunction g0 = new PolynomialFunction(new double[] { k * (k + 1) });
PolynomialFunction g1 = new PolynomialFunction(new double[] { 0, -2 });
PolynomialFunction g2 = new PolynomialFunction(new double[] { 1, 0, -1 });
PolynomialFunction Pk0g0 = Pk0.multiply(g0);
PolynomialFunction Pk1g1 = Pk1.multiply(g1);
PolynomialFunction Pk2g2 = Pk2.multiply(g2);
checkNullPolynomial(Pk0g0.add(Pk1g1.add(Pk2g2)));
}
}
public void testHighDegreeLegendre() {
try {
PolynomialsUtils.createLegendrePolynomial(40);
fail("an exception should have been thrown");
} catch (ArithmeticException ae) {
// expected
}
// 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) {
PolynomialFunction q = new PolynomialFunction(new double[] { denominator});
assertEquals(reference, p.multiply(q).toString());
}
private void checkPolynomial(PolynomialFunction p, String reference) {
assertEquals(reference, p.toString());
}
private void checkNullPolynomial(PolynomialFunction p) {
for (double coefficient : p.getCoefficients()) {
assertEquals(0.0, coefficient, 1.0e-13);
}
}
}