1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math.analysis.polynomials;
18
19 import java.util.ArrayList;
20
21 import org.apache.commons.math.fraction.BigFraction;
22
23 /**
24 * A collection of static methods that operate on or return polynomials.
25 *
26 * @version $Revision: 760901 $ $Date: 2009-04-01 10:29:18 -0400 (Wed, 01 Apr 2009) $
27 * @since 2.0
28 */
29 public class PolynomialsUtils {
30
31 /** Coefficients for Chebyshev polynomials. */
32 private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
33
34 /** Coefficients for Hermite polynomials. */
35 private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
36
37 /** Coefficients for Laguerre polynomials. */
38 private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
39
40 /** Coefficients for Legendre polynomials. */
41 private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
42
43 static {
44
45 // initialize recurrence for Chebyshev polynomials
46 // T0(X) = 1, T1(X) = 0 + 1 * X
47 CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
48 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
49 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
50 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
51
52 // initialize recurrence for Hermite polynomials
53 // H0(X) = 1, H1(X) = 0 + 2 * X
54 HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
55 HERMITE_COEFFICIENTS.add(BigFraction.ONE);
56 HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
57 HERMITE_COEFFICIENTS.add(BigFraction.TWO);
58
59 // initialize recurrence for Laguerre polynomials
60 // L0(X) = 1, L1(X) = 1 - 1 * X
61 LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
62 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
63 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
64 LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
65
66 // initialize recurrence for Legendre polynomials
67 // P0(X) = 1, P1(X) = 0 + 1 * X
68 LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
69 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
70 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
71 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
72
73 }
74
75 /**
76 * Private constructor, to prevent instantiation.
77 */
78 private PolynomialsUtils() {
79 }
80
81 /**
82 * Create a Chebyshev polynomial of the first kind.
83 * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
84 * polynomials of the first kind</a> are orthogonal polynomials.
85 * They can be defined by the following recurrence relations:
86 * <pre>
87 * T<sub>0</sub>(X) = 1
88 * T<sub>1</sub>(X) = X
89 * T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
90 * </pre></p>
91 * @param degree degree of the polynomial
92 * @return Chebyshev polynomial of specified degree
93 */
94 public static PolynomialFunction createChebyshevPolynomial(final int degree) {
95 return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
96 new RecurrenceCoefficientsGenerator() {
97 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
98 /** {@inheritDoc} */
99 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 }