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
023 /**
024 * A collection of static methods that operate on or return polynomials.
025 *
026 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
027 * @since 2.0
028 */
029 public class PolynomialsUtils {
030
031 /** Coefficients for Chebyshev polynomials. */
032 private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
033
034 /** Coefficients for Hermite polynomials. */
035 private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
036
037 /** Coefficients for Laguerre polynomials. */
038 private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
039
040 /** Coefficients for Legendre polynomials. */
041 private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
042
043 static {
044
045 // initialize recurrence for Chebyshev polynomials
046 // T0(X) = 1, T1(X) = 0 + 1 * X
047 CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
048 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
049 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
050 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
051
052 // initialize recurrence for Hermite polynomials
053 // H0(X) = 1, H1(X) = 0 + 2 * X
054 HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
055 HERMITE_COEFFICIENTS.add(BigFraction.ONE);
056 HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
057 HERMITE_COEFFICIENTS.add(BigFraction.TWO);
058
059 // initialize recurrence for Laguerre polynomials
060 // L0(X) = 1, L1(X) = 1 - 1 * X
061 LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
062 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
063 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
064 LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
065
066 // initialize recurrence for Legendre polynomials
067 // P0(X) = 1, P1(X) = 0 + 1 * X
068 LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
069 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
070 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
071 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
072
073 }
074
075 /**
076 * Private constructor, to prevent instantiation.
077 */
078 private PolynomialsUtils() {
079 }
080
081 /**
082 * Create a Chebyshev polynomial of the first kind.
083 * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
084 * polynomials of the first kind</a> are orthogonal polynomials.
085 * They can be defined by the following recurrence relations:
086 * <pre>
087 * T<sub>0</sub>(X) = 1
088 * T<sub>1</sub>(X) = X
089 * T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
090 * </pre></p>
091 * @param degree degree of the polynomial
092 * @return Chebyshev polynomial of specified degree
093 */
094 public static PolynomialFunction createChebyshevPolynomial(final int degree) {
095 return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
096 new RecurrenceCoefficientsGenerator() {
097 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
098 /** {@inheritDoc} */
099 public BigFraction[] generate(int k) {
100 return coeffs;
101 }
102 });
103 }
104
105 /**
106 * Create a Hermite polynomial.
107 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
108 * polynomials</a> are orthogonal polynomials.
109 * They can be defined by the following recurrence relations:
110 * <pre>
111 * H<sub>0</sub>(X) = 1
112 * H<sub>1</sub>(X) = 2X
113 * H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
114 * </pre></p>
115
116 * @param degree degree of the polynomial
117 * @return Hermite polynomial of specified degree
118 */
119 public static PolynomialFunction createHermitePolynomial(final int degree) {
120 return buildPolynomial(degree, HERMITE_COEFFICIENTS,
121 new RecurrenceCoefficientsGenerator() {
122 /** {@inheritDoc} */
123 public BigFraction[] generate(int k) {
124 return new BigFraction[] {
125 BigFraction.ZERO,
126 BigFraction.TWO,
127 new BigFraction(2 * k)};
128 }
129 });
130 }
131
132 /**
133 * Create a Laguerre polynomial.
134 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
135 * polynomials</a> are orthogonal polynomials.
136 * They can be defined by the following recurrence relations:
137 * <pre>
138 * L<sub>0</sub>(X) = 1
139 * L<sub>1</sub>(X) = 1 - X
140 * (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
141 * </pre></p>
142 * @param degree degree of the polynomial
143 * @return Laguerre polynomial of specified degree
144 */
145 public static PolynomialFunction createLaguerrePolynomial(final int degree) {
146 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
147 new RecurrenceCoefficientsGenerator() {
148 /** {@inheritDoc} */
149 public BigFraction[] generate(int k) {
150 final int kP1 = k + 1;
151 return new BigFraction[] {
152 new BigFraction(2 * k + 1, kP1),
153 new BigFraction(-1, kP1),
154 new BigFraction(k, kP1)};
155 }
156 });
157 }
158
159 /**
160 * Create a Legendre polynomial.
161 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
162 * polynomials</a> are orthogonal polynomials.
163 * They can be defined by the following recurrence relations:
164 * <pre>
165 * P<sub>0</sub>(X) = 1
166 * P<sub>1</sub>(X) = X
167 * (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
168 * </pre></p>
169 * @param degree degree of the polynomial
170 * @return Legendre polynomial of specified degree
171 */
172 public static PolynomialFunction createLegendrePolynomial(final int degree) {
173 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
174 new RecurrenceCoefficientsGenerator() {
175 /** {@inheritDoc} */
176 public BigFraction[] generate(int k) {
177 final int kP1 = k + 1;
178 return new BigFraction[] {
179 BigFraction.ZERO,
180 new BigFraction(k + kP1, kP1),
181 new BigFraction(k, kP1)};
182 }
183 });
184 }
185
186 /** Get the coefficients array for a given degree.
187 * @param degree degree of the polynomial
188 * @param coefficients list where the computed coefficients are stored
189 * @param generator recurrence coefficients generator
190 * @return coefficients array
191 */
192 private static PolynomialFunction buildPolynomial(final int degree,
193 final ArrayList<BigFraction> coefficients,
194 final RecurrenceCoefficientsGenerator generator) {
195
196 final int maxDegree = (int) Math.floor(Math.sqrt(2 * coefficients.size())) - 1;
197 synchronized (PolynomialsUtils.class) {
198 if (degree > maxDegree) {
199 computeUpToDegree(degree, maxDegree, generator, coefficients);
200 }
201 }
202
203 // coefficient for polynomial 0 is l [0]
204 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
205 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
206 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
207 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
208 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
209 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
210 // ...
211 final int start = degree * (degree + 1) / 2;
212
213 final double[] a = new double[degree + 1];
214 for (int i = 0; i <= degree; ++i) {
215 a[i] = coefficients.get(start + i).doubleValue();
216 }
217
218 // build the polynomial
219 return new PolynomialFunction(a);
220
221 }
222
223 /** Compute polynomial coefficients up to a given degree.
224 * @param degree maximal degree
225 * @param maxDegree current maximal degree
226 * @param generator recurrence coefficients generator
227 * @param coefficients list where the computed coefficients should be appended
228 */
229 private static void computeUpToDegree(final int degree, final int maxDegree,
230 final RecurrenceCoefficientsGenerator generator,
231 final ArrayList<BigFraction> coefficients) {
232
233 int startK = (maxDegree - 1) * maxDegree / 2;
234 for (int k = maxDegree; k < degree; ++k) {
235
236 // start indices of two previous polynomials Pk(X) and Pk-1(X)
237 int startKm1 = startK;
238 startK += k;
239
240 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
241 BigFraction[] ai = generator.generate(k);
242
243 BigFraction ck = coefficients.get(startK);
244 BigFraction ckm1 = coefficients.get(startKm1);
245
246 // degree 0 coefficient
247 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
248
249 // degree 1 to degree k-1 coefficients
250 for (int i = 1; i < k; ++i) {
251 final BigFraction ckPrev = ck;
252 ck = coefficients.get(startK + i);
253 ckm1 = coefficients.get(startKm1 + i);
254 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
255 }
256
257 // degree k coefficient
258 final BigFraction ckPrev = ck;
259 ck = coefficients.get(startK + k);
260 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
261
262 // degree k+1 coefficient
263 coefficients.add(ck.multiply(ai[1]));
264
265 }
266
267 }
268
269 /** Interface for recurrence coefficients generation. */
270 private static interface RecurrenceCoefficientsGenerator {
271 /**
272 * Generate recurrence coefficients.
273 * @param k highest degree of the polynomials used in the recurrence
274 * @return an array of three coefficients such that
275 * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
276 */
277 BigFraction[] generate(int k);
278 }
279
280 }