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.math3.analysis.polynomials;
018
019 import java.util.ArrayList;
020 import java.util.HashMap;
021 import java.util.List;
022 import java.util.Map;
023
024 import org.apache.commons.math3.fraction.BigFraction;
025 import org.apache.commons.math3.util.ArithmeticUtils;
026 import org.apache.commons.math3.util.FastMath;
027
028 /**
029 * A collection of static methods that operate on or return polynomials.
030 *
031 * @version $Id: PolynomialsUtils.java 1364387 2012-07-22 18:14:11Z tn $
032 * @since 2.0
033 */
034 public class PolynomialsUtils {
035
036 /** Coefficients for Chebyshev polynomials. */
037 private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;
038
039 /** Coefficients for Hermite polynomials. */
040 private static final List<BigFraction> HERMITE_COEFFICIENTS;
041
042 /** Coefficients for Laguerre polynomials. */
043 private static final List<BigFraction> LAGUERRE_COEFFICIENTS;
044
045 /** Coefficients for Legendre polynomials. */
046 private static final List<BigFraction> LEGENDRE_COEFFICIENTS;
047
048 /** Coefficients for Jacobi polynomials. */
049 private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS;
050
051 static {
052
053 // initialize recurrence for Chebyshev polynomials
054 // T0(X) = 1, T1(X) = 0 + 1 * X
055 CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
056 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
057 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
058 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
059
060 // initialize recurrence for Hermite polynomials
061 // H0(X) = 1, H1(X) = 0 + 2 * X
062 HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
063 HERMITE_COEFFICIENTS.add(BigFraction.ONE);
064 HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
065 HERMITE_COEFFICIENTS.add(BigFraction.TWO);
066
067 // initialize recurrence for Laguerre polynomials
068 // L0(X) = 1, L1(X) = 1 - 1 * X
069 LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
070 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
071 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
072 LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
073
074 // initialize recurrence for Legendre polynomials
075 // P0(X) = 1, P1(X) = 0 + 1 * X
076 LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
077 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
078 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
079 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
080
081 // initialize map for Jacobi polynomials
082 JACOBI_COEFFICIENTS = new HashMap<JacobiKey, List<BigFraction>>();
083
084 }
085
086 /**
087 * Private constructor, to prevent instantiation.
088 */
089 private PolynomialsUtils() {
090 }
091
092 /**
093 * Create a Chebyshev polynomial of the first kind.
094 * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
095 * polynomials of the first kind</a> are orthogonal polynomials.
096 * They can be defined by the following recurrence relations:
097 * <pre>
098 * T<sub>0</sub>(X) = 1
099 * T<sub>1</sub>(X) = X
100 * T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
101 * </pre></p>
102 * @param degree degree of the polynomial
103 * @return Chebyshev polynomial of specified degree
104 */
105 public static PolynomialFunction createChebyshevPolynomial(final int degree) {
106 return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
107 new RecurrenceCoefficientsGenerator() {
108 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
109 /** {@inheritDoc} */
110 public BigFraction[] generate(int k) {
111 return coeffs;
112 }
113 });
114 }
115
116 /**
117 * Create a Hermite polynomial.
118 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
119 * polynomials</a> are orthogonal polynomials.
120 * They can be defined by the following recurrence relations:
121 * <pre>
122 * H<sub>0</sub>(X) = 1
123 * H<sub>1</sub>(X) = 2X
124 * H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
125 * </pre></p>
126
127 * @param degree degree of the polynomial
128 * @return Hermite polynomial of specified degree
129 */
130 public static PolynomialFunction createHermitePolynomial(final int degree) {
131 return buildPolynomial(degree, HERMITE_COEFFICIENTS,
132 new RecurrenceCoefficientsGenerator() {
133 /** {@inheritDoc} */
134 public BigFraction[] generate(int k) {
135 return new BigFraction[] {
136 BigFraction.ZERO,
137 BigFraction.TWO,
138 new BigFraction(2 * k)};
139 }
140 });
141 }
142
143 /**
144 * Create a Laguerre polynomial.
145 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
146 * polynomials</a> are orthogonal polynomials.
147 * They can be defined by the following recurrence relations:
148 * <pre>
149 * L<sub>0</sub>(X) = 1
150 * L<sub>1</sub>(X) = 1 - X
151 * (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
152 * </pre></p>
153 * @param degree degree of the polynomial
154 * @return Laguerre polynomial of specified degree
155 */
156 public static PolynomialFunction createLaguerrePolynomial(final int degree) {
157 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
158 new RecurrenceCoefficientsGenerator() {
159 /** {@inheritDoc} */
160 public BigFraction[] generate(int k) {
161 final int kP1 = k + 1;
162 return new BigFraction[] {
163 new BigFraction(2 * k + 1, kP1),
164 new BigFraction(-1, kP1),
165 new BigFraction(k, kP1)};
166 }
167 });
168 }
169
170 /**
171 * Create a Legendre polynomial.
172 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
173 * polynomials</a> are orthogonal polynomials.
174 * They can be defined by the following recurrence relations:
175 * <pre>
176 * P<sub>0</sub>(X) = 1
177 * P<sub>1</sub>(X) = X
178 * (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
179 * </pre></p>
180 * @param degree degree of the polynomial
181 * @return Legendre polynomial of specified degree
182 */
183 public static PolynomialFunction createLegendrePolynomial(final int degree) {
184 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
185 new RecurrenceCoefficientsGenerator() {
186 /** {@inheritDoc} */
187 public BigFraction[] generate(int k) {
188 final int kP1 = k + 1;
189 return new BigFraction[] {
190 BigFraction.ZERO,
191 new BigFraction(k + kP1, kP1),
192 new BigFraction(k, kP1)};
193 }
194 });
195 }
196
197 /**
198 * Create a Jacobi polynomial.
199 * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
200 * polynomials</a> are orthogonal polynomials.
201 * They can be defined by the following recurrence relations:
202 * <pre>
203 * P<sub>0</sub><sup>vw</sup>(X) = 1
204 * P<sub>-1</sub><sup>vw</sup>(X) = 0
205 * 2k(k + v + w)(2k + v + w - 2) P<sub>k</sub><sup>vw</sup>(X) =
206 * (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) X + v<sup>2</sup> - w<sup>2</sup>] P<sub>k-1</sub><sup>vw</sup>(X)
207 * - 2(k + v - 1)(k + w - 1)(2k + v + w) P<sub>k-2</sub><sup>vw</sup>(X)
208 * </pre></p>
209 * @param degree degree of the polynomial
210 * @param v first exponent
211 * @param w second exponent
212 * @return Jacobi polynomial of specified degree
213 */
214 public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {
215
216 // select the appropriate list
217 final JacobiKey key = new JacobiKey(v, w);
218
219 if (!JACOBI_COEFFICIENTS.containsKey(key)) {
220
221 // allocate a new list for v, w
222 final List<BigFraction> list = new ArrayList<BigFraction>();
223 JACOBI_COEFFICIENTS.put(key, list);
224
225 // Pv,w,0(x) = 1;
226 list.add(BigFraction.ONE);
227
228 // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
229 list.add(new BigFraction(v - w, 2));
230 list.add(new BigFraction(2 + v + w, 2));
231
232 }
233
234 return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
235 new RecurrenceCoefficientsGenerator() {
236 /** {@inheritDoc} */
237 public BigFraction[] generate(int k) {
238 k++;
239 final int kvw = k + v + w;
240 final int twoKvw = kvw + k;
241 final int twoKvwM1 = twoKvw - 1;
242 final int twoKvwM2 = twoKvw - 2;
243 final int den = 2 * k * kvw * twoKvwM2;
244
245 return new BigFraction[] {
246 new BigFraction(twoKvwM1 * (v * v - w * w), den),
247 new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den),
248 new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
249 };
250 }
251 });
252
253 }
254
255 /** Inner class for Jacobi polynomials keys. */
256 private static class JacobiKey {
257
258 /** First exponent. */
259 private final int v;
260
261 /** Second exponent. */
262 private final int w;
263
264 /** Simple constructor.
265 * @param v first exponent
266 * @param w second exponent
267 */
268 public JacobiKey(final int v, final int w) {
269 this.v = v;
270 this.w = w;
271 }
272
273 /** Get hash code.
274 * @return hash code
275 */
276 @Override
277 public int hashCode() {
278 return (v << 16) ^ w;
279 }
280
281 /** Check if the instance represent the same key as another instance.
282 * @param key other key
283 * @return true if the instance and the other key refer to the same polynomial
284 */
285 @Override
286 public boolean equals(final Object key) {
287
288 if ((key == null) || !(key instanceof JacobiKey)) {
289 return false;
290 }
291
292 final JacobiKey otherK = (JacobiKey) key;
293 return (v == otherK.v) && (w == otherK.w);
294
295 }
296 }
297
298 /**
299 * Compute the coefficients of the polynomial <code>P<sub>s</sub>(x)</code>
300 * whose values at point {@code x} will be the same as the those from the
301 * original polynomial <code>P(x)</code> when computed at {@code x + shift}.
302 * Thus, if <code>P(x) = Σ<sub>i</sub> a<sub>i</sub> x<sup>i</sup></code>,
303 * then
304 * <pre>
305 * <table>
306 * <tr>
307 * <td><code>P<sub>s</sub>(x)</td>
308 * <td>= Σ<sub>i</sub> b<sub>i</sub> x<sup>i</sup></code></td>
309 * </tr>
310 * <tr>
311 * <td></td>
312 * <td>= Σ<sub>i</sub> a<sub>i</sub> (x + shift)<sup>i</sup></code></td>
313 * </tr>
314 * </table>
315 * </pre>
316 *
317 * @param coefficients Coefficients of the original polynomial.
318 * @param shift Shift value.
319 * @return the coefficients <code>b<sub>i</sub></code> of the shifted
320 * polynomial.
321 */
322 public static double[] shift(final double[] coefficients,
323 final double shift) {
324 final int dp1 = coefficients.length;
325 final double[] newCoefficients = new double[dp1];
326
327 // Pascal triangle.
328 final int[][] coeff = new int[dp1][dp1];
329 for (int i = 0; i < dp1; i++){
330 for(int j = 0; j <= i; j++){
331 coeff[i][j] = (int) ArithmeticUtils.binomialCoefficient(i, j);
332 }
333 }
334
335 // First polynomial coefficient.
336 for (int i = 0; i < dp1; i++){
337 newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
338 }
339
340 // Superior order.
341 final int d = dp1 - 1;
342 for (int i = 0; i < d; i++) {
343 for (int j = i; j < d; j++){
344 newCoefficients[i + 1] += coeff[j + 1][j - i] *
345 coefficients[j + 1] * FastMath.pow(shift, j - i);
346 }
347 }
348
349 return newCoefficients;
350 }
351
352
353 /** Get the coefficients array for a given degree.
354 * @param degree degree of the polynomial
355 * @param coefficients list where the computed coefficients are stored
356 * @param generator recurrence coefficients generator
357 * @return coefficients array
358 */
359 private static PolynomialFunction buildPolynomial(final int degree,
360 final List<BigFraction> coefficients,
361 final RecurrenceCoefficientsGenerator generator) {
362
363 final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
364 synchronized (PolynomialsUtils.class) {
365 if (degree > maxDegree) {
366 computeUpToDegree(degree, maxDegree, generator, coefficients);
367 }
368 }
369
370 // coefficient for polynomial 0 is l [0]
371 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
372 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
373 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
374 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
375 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
376 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
377 // ...
378 final int start = degree * (degree + 1) / 2;
379
380 final double[] a = new double[degree + 1];
381 for (int i = 0; i <= degree; ++i) {
382 a[i] = coefficients.get(start + i).doubleValue();
383 }
384
385 // build the polynomial
386 return new PolynomialFunction(a);
387
388 }
389
390 /** Compute polynomial coefficients up to a given degree.
391 * @param degree maximal degree
392 * @param maxDegree current maximal degree
393 * @param generator recurrence coefficients generator
394 * @param coefficients list where the computed coefficients should be appended
395 */
396 private static void computeUpToDegree(final int degree, final int maxDegree,
397 final RecurrenceCoefficientsGenerator generator,
398 final List<BigFraction> coefficients) {
399
400 int startK = (maxDegree - 1) * maxDegree / 2;
401 for (int k = maxDegree; k < degree; ++k) {
402
403 // start indices of two previous polynomials Pk(X) and Pk-1(X)
404 int startKm1 = startK;
405 startK += k;
406
407 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
408 BigFraction[] ai = generator.generate(k);
409
410 BigFraction ck = coefficients.get(startK);
411 BigFraction ckm1 = coefficients.get(startKm1);
412
413 // degree 0 coefficient
414 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
415
416 // degree 1 to degree k-1 coefficients
417 for (int i = 1; i < k; ++i) {
418 final BigFraction ckPrev = ck;
419 ck = coefficients.get(startK + i);
420 ckm1 = coefficients.get(startKm1 + i);
421 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
422 }
423
424 // degree k coefficient
425 final BigFraction ckPrev = ck;
426 ck = coefficients.get(startK + k);
427 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
428
429 // degree k+1 coefficient
430 coefficients.add(ck.multiply(ai[1]));
431
432 }
433
434 }
435
436 /** Interface for recurrence coefficients generation. */
437 private interface RecurrenceCoefficientsGenerator {
438 /**
439 * Generate recurrence coefficients.
440 * @param k highest degree of the polynomials used in the recurrence
441 * @return an array of three coefficients such that
442 * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
443 */
444 BigFraction[] generate(int k);
445 }
446
447 }