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: 760901 $ $Date: 2009-04-01 10:29:18 -0400 (Wed, 01 Apr 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 }