001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.analysis.polynomials;
018
019 import java.util.ArrayList;
020
021 import org.apache.commons.math.fraction.BigFraction;
022 import org.apache.commons.math.util.FastMath;
023
024 /**
025 * A collection of static methods that operate on or return polynomials.
026 *
027 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
028 * @since 2.0
029 */
030 public class PolynomialsUtils {
031
032 /** Coefficients for Chebyshev polynomials. */
033 private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
034
035 /** Coefficients for Hermite polynomials. */
036 private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
037
038 /** Coefficients for Laguerre polynomials. */
039 private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
040
041 /** Coefficients for Legendre polynomials. */
042 private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
043
044 static {
045
046 // initialize recurrence for Chebyshev polynomials
047 // T0(X) = 1, T1(X) = 0 + 1 * X
048 CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
049 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
050 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
051 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
052
053 // initialize recurrence for Hermite polynomials
054 // H0(X) = 1, H1(X) = 0 + 2 * X
055 HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
056 HERMITE_COEFFICIENTS.add(BigFraction.ONE);
057 HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
058 HERMITE_COEFFICIENTS.add(BigFraction.TWO);
059
060 // initialize recurrence for Laguerre polynomials
061 // L0(X) = 1, L1(X) = 1 - 1 * X
062 LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
063 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
064 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
065 LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
066
067 // initialize recurrence for Legendre polynomials
068 // P0(X) = 1, P1(X) = 0 + 1 * X
069 LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
070 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
071 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
072 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
073
074 }
075
076 /**
077 * Private constructor, to prevent instantiation.
078 */
079 private PolynomialsUtils() {
080 }
081
082 /**
083 * Create a Chebyshev polynomial of the first kind.
084 * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
085 * polynomials of the first kind</a> are orthogonal polynomials.
086 * They can be defined by the following recurrence relations:
087 * <pre>
088 * T<sub>0</sub>(X) = 1
089 * T<sub>1</sub>(X) = X
090 * T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
091 * </pre></p>
092 * @param degree degree of the polynomial
093 * @return Chebyshev polynomial of specified degree
094 */
095 public static PolynomialFunction createChebyshevPolynomial(final int degree) {
096 return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
097 new RecurrenceCoefficientsGenerator() {
098 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
099 /** {@inheritDoc} */
100 public BigFraction[] generate(int k) {
101 return coeffs;
102 }
103 });
104 }
105
106 /**
107 * Create a Hermite polynomial.
108 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
109 * polynomials</a> are orthogonal polynomials.
110 * They can be defined by the following recurrence relations:
111 * <pre>
112 * H<sub>0</sub>(X) = 1
113 * H<sub>1</sub>(X) = 2X
114 * H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
115 * </pre></p>
116
117 * @param degree degree of the polynomial
118 * @return Hermite polynomial of specified degree
119 */
120 public static PolynomialFunction createHermitePolynomial(final int degree) {
121 return buildPolynomial(degree, HERMITE_COEFFICIENTS,
122 new RecurrenceCoefficientsGenerator() {
123 /** {@inheritDoc} */
124 public BigFraction[] generate(int k) {
125 return new BigFraction[] {
126 BigFraction.ZERO,
127 BigFraction.TWO,
128 new BigFraction(2 * k)};
129 }
130 });
131 }
132
133 /**
134 * Create a Laguerre polynomial.
135 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
136 * polynomials</a> are orthogonal polynomials.
137 * They can be defined by the following recurrence relations:
138 * <pre>
139 * L<sub>0</sub>(X) = 1
140 * L<sub>1</sub>(X) = 1 - X
141 * (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
142 * </pre></p>
143 * @param degree degree of the polynomial
144 * @return Laguerre polynomial of specified degree
145 */
146 public static PolynomialFunction createLaguerrePolynomial(final int degree) {
147 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
148 new RecurrenceCoefficientsGenerator() {
149 /** {@inheritDoc} */
150 public BigFraction[] generate(int k) {
151 final int kP1 = k + 1;
152 return new BigFraction[] {
153 new BigFraction(2 * k + 1, kP1),
154 new BigFraction(-1, kP1),
155 new BigFraction(k, kP1)};
156 }
157 });
158 }
159
160 /**
161 * Create a Legendre polynomial.
162 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
163 * polynomials</a> are orthogonal polynomials.
164 * They can be defined by the following recurrence relations:
165 * <pre>
166 * P<sub>0</sub>(X) = 1
167 * P<sub>1</sub>(X) = X
168 * (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
169 * </pre></p>
170 * @param degree degree of the polynomial
171 * @return Legendre polynomial of specified degree
172 */
173 public static PolynomialFunction createLegendrePolynomial(final int degree) {
174 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
175 new RecurrenceCoefficientsGenerator() {
176 /** {@inheritDoc} */
177 public BigFraction[] generate(int k) {
178 final int kP1 = k + 1;
179 return new BigFraction[] {
180 BigFraction.ZERO,
181 new BigFraction(k + kP1, kP1),
182 new BigFraction(k, kP1)};
183 }
184 });
185 }
186
187 /** Get the coefficients array for a given degree.
188 * @param degree degree of the polynomial
189 * @param coefficients list where the computed coefficients are stored
190 * @param generator recurrence coefficients generator
191 * @return coefficients array
192 */
193 private static PolynomialFunction buildPolynomial(final int degree,
194 final ArrayList<BigFraction> coefficients,
195 final RecurrenceCoefficientsGenerator generator) {
196
197 final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
198 synchronized (PolynomialsUtils.class) {
199 if (degree > maxDegree) {
200 computeUpToDegree(degree, maxDegree, generator, coefficients);
201 }
202 }
203
204 // coefficient for polynomial 0 is l [0]
205 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
206 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
207 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
208 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
209 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
210 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
211 // ...
212 final int start = degree * (degree + 1) / 2;
213
214 final double[] a = new double[degree + 1];
215 for (int i = 0; i <= degree; ++i) {
216 a[i] = coefficients.get(start + i).doubleValue();
217 }
218
219 // build the polynomial
220 return new PolynomialFunction(a);
221
222 }
223
224 /** Compute polynomial coefficients up to a given degree.
225 * @param degree maximal degree
226 * @param maxDegree current maximal degree
227 * @param generator recurrence coefficients generator
228 * @param coefficients list where the computed coefficients should be appended
229 */
230 private static void computeUpToDegree(final int degree, final int maxDegree,
231 final RecurrenceCoefficientsGenerator generator,
232 final ArrayList<BigFraction> coefficients) {
233
234 int startK = (maxDegree - 1) * maxDegree / 2;
235 for (int k = maxDegree; k < degree; ++k) {
236
237 // start indices of two previous polynomials Pk(X) and Pk-1(X)
238 int startKm1 = startK;
239 startK += k;
240
241 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
242 BigFraction[] ai = generator.generate(k);
243
244 BigFraction ck = coefficients.get(startK);
245 BigFraction ckm1 = coefficients.get(startKm1);
246
247 // degree 0 coefficient
248 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
249
250 // degree 1 to degree k-1 coefficients
251 for (int i = 1; i < k; ++i) {
252 final BigFraction ckPrev = ck;
253 ck = coefficients.get(startK + i);
254 ckm1 = coefficients.get(startKm1 + i);
255 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
256 }
257
258 // degree k coefficient
259 final BigFraction ckPrev = ck;
260 ck = coefficients.get(startK + k);
261 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
262
263 // degree k+1 coefficient
264 coefficients.add(ck.multiply(ai[1]));
265
266 }
267
268 }
269
270 /** Interface for recurrence coefficients generation. */
271 private static interface RecurrenceCoefficientsGenerator {
272 /**
273 * Generate recurrence coefficients.
274 * @param k highest degree of the polynomials used in the recurrence
275 * @return an array of three coefficients such that
276 * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
277 */
278 BigFraction[] generate(int k);
279 }
280
281 }