source: src/main/java/agents/anac/y2019/harddealer/math3/analysis/polynomials/PolynomialsUtils.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 16.2 KB
Line 
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 */
17package agents.anac.y2019.harddealer.math3.analysis.polynomials;
18
19import java.util.ArrayList;
20import java.util.HashMap;
21import java.util.List;
22import java.util.Map;
23
24import agents.anac.y2019.harddealer.math3.fraction.BigFraction;
25import agents.anac.y2019.harddealer.math3.util.CombinatoricsUtils;
26import agents.anac.y2019.harddealer.math3.util.FastMath;
27
28/**
29 * A collection of static methods that operate on or return polynomials.
30 *
31 * @since 2.0
32 */
33public class PolynomialsUtils {
34
35 /** Coefficients for Chebyshev polynomials. */
36 private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;
37
38 /** Coefficients for Hermite polynomials. */
39 private static final List<BigFraction> HERMITE_COEFFICIENTS;
40
41 /** Coefficients for Laguerre polynomials. */
42 private static final List<BigFraction> LAGUERRE_COEFFICIENTS;
43
44 /** Coefficients for Legendre polynomials. */
45 private static final List<BigFraction> LEGENDRE_COEFFICIENTS;
46
47 /** Coefficients for Jacobi polynomials. */
48 private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS;
49
50 static {
51
52 // initialize recurrence for Chebyshev polynomials
53 // T0(X) = 1, T1(X) = 0 + 1 * X
54 CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
55 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
56 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
57 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
58
59 // initialize recurrence for Hermite polynomials
60 // H0(X) = 1, H1(X) = 0 + 2 * X
61 HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
62 HERMITE_COEFFICIENTS.add(BigFraction.ONE);
63 HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
64 HERMITE_COEFFICIENTS.add(BigFraction.TWO);
65
66 // initialize recurrence for Laguerre polynomials
67 // L0(X) = 1, L1(X) = 1 - 1 * X
68 LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
69 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
70 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
71 LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
72
73 // initialize recurrence for Legendre polynomials
74 // P0(X) = 1, P1(X) = 0 + 1 * X
75 LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
76 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
77 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
78 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
79
80 // initialize map for Jacobi polynomials
81 JACOBI_COEFFICIENTS = new HashMap<JacobiKey, List<BigFraction>>();
82
83 }
84
85 /**
86 * Private constructor, to prevent instantiation.
87 */
88 private PolynomialsUtils() {
89 }
90
91 /**
92 * Create a Chebyshev polynomial of the first kind.
93 * <p><a href="https://en.wikipedia.org/wiki/Chebyshev_polynomials">Chebyshev
94 * polynomials of the first kind</a> are orthogonal polynomials.
95 * They can be defined by the following recurrence relations:</p><p>
96 * \(
97 * T_0(x) = 1 \\
98 * T_1(x) = x \\
99 * T_{k+1}(x) = 2x T_k(x) - T_{k-1}(x)
100 * \)
101 * </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 /** Fixed recurrence coefficients. */
109 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
110 /** {@inheritDoc} */
111 public BigFraction[] generate(int k) {
112 return coeffs;
113 }
114 });
115 }
116
117 /**
118 * Create a Hermite polynomial.
119 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
120 * polynomials</a> are orthogonal polynomials.
121 * They can be defined by the following recurrence relations:</p><p>
122 * \(
123 * H_0(x) = 1 \\
124 * H_1(x) = 2x \\
125 * H_{k+1}(x) = 2x H_k(X) - 2k H_{k-1}(x)
126 * \)
127 * </p>
128
129 * @param degree degree of the polynomial
130 * @return Hermite polynomial of specified degree
131 */
132 public static PolynomialFunction createHermitePolynomial(final int degree) {
133 return buildPolynomial(degree, HERMITE_COEFFICIENTS,
134 new RecurrenceCoefficientsGenerator() {
135 /** {@inheritDoc} */
136 public BigFraction[] generate(int k) {
137 return new BigFraction[] {
138 BigFraction.ZERO,
139 BigFraction.TWO,
140 new BigFraction(2 * k)};
141 }
142 });
143 }
144
145 /**
146 * Create a Laguerre polynomial.
147 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
148 * polynomials</a> are orthogonal polynomials.
149 * They can be defined by the following recurrence relations:</p><p>
150 * \(
151 * L_0(x) = 1 \\
152 * L_1(x) = 1 - x \\
153 * (k+1) L_{k+1}(x) = (2k + 1 - x) L_k(x) - k L_{k-1}(x)
154 * \)
155 * </p>
156 * @param degree degree of the polynomial
157 * @return Laguerre polynomial of specified degree
158 */
159 public static PolynomialFunction createLaguerrePolynomial(final int degree) {
160 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
161 new RecurrenceCoefficientsGenerator() {
162 /** {@inheritDoc} */
163 public BigFraction[] generate(int k) {
164 final int kP1 = k + 1;
165 return new BigFraction[] {
166 new BigFraction(2 * k + 1, kP1),
167 new BigFraction(-1, kP1),
168 new BigFraction(k, kP1)};
169 }
170 });
171 }
172
173 /**
174 * Create a Legendre polynomial.
175 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
176 * polynomials</a> are orthogonal polynomials.
177 * They can be defined by the following recurrence relations:</p><p>
178 * \(
179 * P_0(x) = 1 \\
180 * P_1(x) = x \\
181 * (k+1) P_{k+1}(x) = (2k+1) x P_k(x) - k P_{k-1}(x)
182 * \)
183 * </p>
184 * @param degree degree of the polynomial
185 * @return Legendre polynomial of specified degree
186 */
187 public static PolynomialFunction createLegendrePolynomial(final int degree) {
188 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
189 new RecurrenceCoefficientsGenerator() {
190 /** {@inheritDoc} */
191 public BigFraction[] generate(int k) {
192 final int kP1 = k + 1;
193 return new BigFraction[] {
194 BigFraction.ZERO,
195 new BigFraction(k + kP1, kP1),
196 new BigFraction(k, kP1)};
197 }
198 });
199 }
200
201 /**
202 * Create a Jacobi polynomial.
203 * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
204 * polynomials</a> are orthogonal polynomials.
205 * They can be defined by the following recurrence relations:</p><p>
206 * \(
207 * P_0^{vw}(x) = 1 \\
208 * P_{-1}^{vw}(x) = 0 \\
209 * 2k(k + v + w)(2k + v + w - 2) P_k^{vw}(x) = \\
210 * (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) x + v^2 - w^2] P_{k-1}^{vw}(x) \\
211 * - 2(k + v - 1)(k + w - 1)(2k + v + w) P_{k-2}^{vw}(x)
212 * \)
213 * </p>
214 * @param degree degree of the polynomial
215 * @param v first exponent
216 * @param w second exponent
217 * @return Jacobi polynomial of specified degree
218 */
219 public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {
220
221 // select the appropriate list
222 final JacobiKey key = new JacobiKey(v, w);
223
224 if (!JACOBI_COEFFICIENTS.containsKey(key)) {
225
226 // allocate a new list for v, w
227 final List<BigFraction> list = new ArrayList<BigFraction>();
228 JACOBI_COEFFICIENTS.put(key, list);
229
230 // Pv,w,0(x) = 1;
231 list.add(BigFraction.ONE);
232
233 // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
234 list.add(new BigFraction(v - w, 2));
235 list.add(new BigFraction(2 + v + w, 2));
236
237 }
238
239 return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
240 new RecurrenceCoefficientsGenerator() {
241 /** {@inheritDoc} */
242 public BigFraction[] generate(int k) {
243 k++;
244 final int kvw = k + v + w;
245 final int twoKvw = kvw + k;
246 final int twoKvwM1 = twoKvw - 1;
247 final int twoKvwM2 = twoKvw - 2;
248 final int den = 2 * k * kvw * twoKvwM2;
249
250 return new BigFraction[] {
251 new BigFraction(twoKvwM1 * (v * v - w * w), den),
252 new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den),
253 new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
254 };
255 }
256 });
257
258 }
259
260 /** Inner class for Jacobi polynomials keys. */
261 private static class JacobiKey {
262
263 /** First exponent. */
264 private final int v;
265
266 /** Second exponent. */
267 private final int w;
268
269 /** Simple constructor.
270 * @param v first exponent
271 * @param w second exponent
272 */
273 JacobiKey(final int v, final int w) {
274 this.v = v;
275 this.w = w;
276 }
277
278 /** Get hash code.
279 * @return hash code
280 */
281 @Override
282 public int hashCode() {
283 return (v << 16) ^ w;
284 }
285
286 /** Check if the instance represent the same key as another instance.
287 * @param key other key
288 * @return true if the instance and the other key refer to the same polynomial
289 */
290 @Override
291 public boolean equals(final Object key) {
292
293 if ((key == null) || !(key instanceof JacobiKey)) {
294 return false;
295 }
296
297 final JacobiKey otherK = (JacobiKey) key;
298 return (v == otherK.v) && (w == otherK.w);
299
300 }
301 }
302
303 /**
304 * Compute the coefficients of the polynomial \(P_s(x)\)
305 * whose values at point {@code x} will be the same as the those from the
306 * original polynomial \(P(x)\) when computed at {@code x + shift}.
307 * <p>
308 * More precisely, let \(\Delta = \) {@code shift} and let
309 * \(P_s(x) = P(x + \Delta)\). The returned array
310 * consists of the coefficients of \(P_s\). So if \(a_0, ..., a_{n-1}\)
311 * are the coefficients of \(P\), then the returned array
312 * \(b_0, ..., b_{n-1}\) satisfies the identity
313 * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
314 *
315 * @param coefficients Coefficients of the original polynomial.
316 * @param shift Shift value.
317 * @return the coefficients \(b_i\) of the shifted
318 * polynomial.
319 */
320 public static double[] shift(final double[] coefficients,
321 final double shift) {
322 final int dp1 = coefficients.length;
323 final double[] newCoefficients = new double[dp1];
324
325 // Pascal triangle.
326 final int[][] coeff = new int[dp1][dp1];
327 for (int i = 0; i < dp1; i++){
328 for(int j = 0; j <= i; j++){
329 coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j);
330 }
331 }
332
333 // First polynomial coefficient.
334 for (int i = 0; i < dp1; i++){
335 newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
336 }
337
338 // Superior order.
339 final int d = dp1 - 1;
340 for (int i = 0; i < d; i++) {
341 for (int j = i; j < d; j++){
342 newCoefficients[i + 1] += coeff[j + 1][j - i] *
343 coefficients[j + 1] * FastMath.pow(shift, j - i);
344 }
345 }
346
347 return newCoefficients;
348 }
349
350
351 /** Get the coefficients array for a given degree.
352 * @param degree degree of the polynomial
353 * @param coefficients list where the computed coefficients are stored
354 * @param generator recurrence coefficients generator
355 * @return coefficients array
356 */
357 private static PolynomialFunction buildPolynomial(final int degree,
358 final List<BigFraction> coefficients,
359 final RecurrenceCoefficientsGenerator generator) {
360 synchronized (coefficients) {
361 final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
362 if (degree > maxDegree) {
363 computeUpToDegree(degree, maxDegree, generator, coefficients);
364 }
365 }
366
367 // coefficient for polynomial 0 is l [0]
368 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
369 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
370 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
371 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
372 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
373 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
374 // ...
375 final int start = degree * (degree + 1) / 2;
376
377 final double[] a = new double[degree + 1];
378 for (int i = 0; i <= degree; ++i) {
379 a[i] = coefficients.get(start + i).doubleValue();
380 }
381
382 // build the polynomial
383 return new PolynomialFunction(a);
384
385 }
386
387 /** Compute polynomial coefficients up to a given degree.
388 * @param degree maximal degree
389 * @param maxDegree current maximal degree
390 * @param generator recurrence coefficients generator
391 * @param coefficients list where the computed coefficients should be appended
392 */
393 private static void computeUpToDegree(final int degree, final int maxDegree,
394 final RecurrenceCoefficientsGenerator generator,
395 final List<BigFraction> coefficients) {
396
397 int startK = (maxDegree - 1) * maxDegree / 2;
398 for (int k = maxDegree; k < degree; ++k) {
399
400 // start indices of two previous polynomials Pk(X) and Pk-1(X)
401 int startKm1 = startK;
402 startK += k;
403
404 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
405 BigFraction[] ai = generator.generate(k);
406
407 BigFraction ck = coefficients.get(startK);
408 BigFraction ckm1 = coefficients.get(startKm1);
409
410 // degree 0 coefficient
411 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
412
413 // degree 1 to degree k-1 coefficients
414 for (int i = 1; i < k; ++i) {
415 final BigFraction ckPrev = ck;
416 ck = coefficients.get(startK + i);
417 ckm1 = coefficients.get(startKm1 + i);
418 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
419 }
420
421 // degree k coefficient
422 final BigFraction ckPrev = ck;
423 ck = coefficients.get(startK + k);
424 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
425
426 // degree k+1 coefficient
427 coefficients.add(ck.multiply(ai[1]));
428
429 }
430
431 }
432
433 /** Interface for recurrence coefficients generation. */
434 private interface RecurrenceCoefficientsGenerator {
435 /**
436 * Generate recurrence coefficients.
437 * @param k highest degree of the polynomials used in the recurrence
438 * @return an array of three coefficients such that
439 * \( P_{k+1}(x) = (a[0] + a[1] x) P_k(x) - a[2] P_{k-1}(x) \)
440 */
441 BigFraction[] generate(int k);
442 }
443
444}
Note: See TracBrowser for help on using the repository browser.