source: src/main/java/agents/anac/y2019/harddealer/math3/special/BesselJ.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: 27.8 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 */
17
18package agents.anac.y2019.harddealer.math3.special;
19
20import agents.anac.y2019.harddealer.math3.analysis.UnivariateFunction;
21import agents.anac.y2019.harddealer.math3.exception.ConvergenceException;
22import agents.anac.y2019.harddealer.math3.exception.MathIllegalArgumentException;
23import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
24import agents.anac.y2019.harddealer.math3.util.FastMath;
25import agents.anac.y2019.harddealer.math3.util.MathArrays;
26
27/**
28 * This class provides computation methods related to Bessel
29 * functions of the first kind. Detailed descriptions of these functions are
30 * available in <a
31 * href="http://en.wikipedia.org/wiki/Bessel_function">Wikipedia</a>, <a
32 * href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Abrabowitz and
33 * Stegun</a> (Ch. 9-11), and <a href="http://dlmf.nist.gov/">DLMF</a> (Ch. 10).
34 * <p>
35 * This implementation is based on the rjbesl Fortran routine at
36 * <a href="http://www.netlib.org/specfun/rjbesl">Netlib</a>.</p>
37 * <p>
38 * From the Fortran code: </p>
39 * <p>
40 * This program is based on a program written by David J. Sookne (2) that
41 * computes values of the Bessel functions J or I of real argument and integer
42 * order. Modifications include the restriction of the computation to the J
43 * Bessel function of non-negative real argument, the extension of the
44 * computation to arbitrary positive order, and the elimination of most
45 * underflow.</p>
46 * <p>
47 * References:</p>
48 * <ul>
49 * <li>"A Note on Backward Recurrence Algorithms," Olver, F. W. J., and Sookne,
50 * D. J., Math. Comp. 26, 1972, pp 941-947.</li>
51 * <li>"Bessel Functions of Real Argument and Integer Order," Sookne, D. J., NBS
52 * Jour. of Res. B. 77B, 1973, pp 125-132.</li>
53 * </ul> </p>
54 * @since 3.4
55 */
56public class BesselJ
57 implements UnivariateFunction {
58
59 // ---------------------------------------------------------------------
60 // Mathematical constants
61 // ---------------------------------------------------------------------
62
63 /** -2 / pi */
64 private static final double PI2 = 0.636619772367581343075535;
65
66 /** first few significant digits of 2pi */
67 private static final double TOWPI1 = 6.28125;
68
69 /** 2pi - TWOPI1 to working precision */
70 private static final double TWOPI2 = 1.935307179586476925286767e-3;
71
72 /** TOWPI1 + TWOPI2 */
73 private static final double TWOPI = TOWPI1 + TWOPI2;
74
75 // ---------------------------------------------------------------------
76 // Machine-dependent parameters
77 // ---------------------------------------------------------------------
78
79 /**
80 * 10.0^K, where K is the largest integer such that ENTEN is
81 * machine-representable in working precision
82 */
83 private static final double ENTEN = 1.0e308;
84
85 /**
86 * Decimal significance desired. Should be set to (INT(log_{10}(2) * (it)+1)).
87 * Setting NSIG lower will result in decreased accuracy while setting
88 * NSIG higher will increase CPU time without increasing accuracy.
89 * The truncation error is limited to a relative error of
90 * T=.5(10^(-NSIG)).
91 */
92 private static final double ENSIG = 1.0e16;
93
94 /** 10.0 ** (-K) for the smallest integer K such that K >= NSIG/4 */
95 private static final double RTNSIG = 1.0e-4;
96
97 /** Smallest ABS(X) such that X/4 does not underflow */
98 private static final double ENMTEN = 8.90e-308;
99
100 /** Minimum acceptable value for x */
101 private static final double X_MIN = 0.0;
102
103 /**
104 * Upper limit on the magnitude of x. If abs(x) = n, then at least
105 * n iterations of the backward recursion will be executed. The value of
106 * 10.0 ** 4 is used on every machine.
107 */
108 private static final double X_MAX = 1.0e4;
109
110 /** First 25 factorials as doubles */
111 private static final double[] FACT = {
112 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
113 3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
114 1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15,
115 1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19,
116 1.12400072777760768e21, 2.585201673888497664e22,
117 6.2044840173323943936e23
118 };
119
120 /** Order of the function computed when {@link #value(double)} is used */
121 private final double order;
122
123 /**
124 * Create a new BesselJ with the given order.
125 *
126 * @param order order of the function computed when using {@link #value(double)}.
127 */
128 public BesselJ(double order) {
129 this.order = order;
130 }
131
132 /**
133 * Returns the value of the constructed Bessel function of the first kind,
134 * for the passed argument.
135 *
136 * @param x Argument
137 * @return Value of the Bessel function at x
138 * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
139 * @throws ConvergenceException if the algorithm fails to converge
140 */
141 public double value(double x)
142 throws MathIllegalArgumentException, ConvergenceException {
143 return BesselJ.value(order, x);
144 }
145
146 /**
147 * Returns the first Bessel function, \(J_{order}(x)\).
148 *
149 * @param order Order of the Bessel function
150 * @param x Argument
151 * @return Value of the Bessel function of the first kind, \(J_{order}(x)\)
152 * @throws MathIllegalArgumentException if {@code x} is too large relative to {@code order}
153 * @throws ConvergenceException if the algorithm fails to converge
154 */
155 public static double value(double order, double x)
156 throws MathIllegalArgumentException, ConvergenceException {
157 final int n = (int) order;
158 final double alpha = order - n;
159 final int nb = n + 1;
160 final BesselJResult res = rjBesl(x, alpha, nb);
161
162 if (res.nVals >= nb) {
163 return res.vals[n];
164 } else if (res.nVals < 0) {
165 throw new MathIllegalArgumentException(LocalizedFormats.BESSEL_FUNCTION_BAD_ARGUMENT,order, x);
166 } else if (FastMath.abs(res.vals[res.nVals - 1]) < 1e-100) {
167 return res.vals[n]; // underflow; return value (will be zero)
168 }
169 throw new ConvergenceException(LocalizedFormats.BESSEL_FUNCTION_FAILED_CONVERGENCE, order, x);
170 }
171
172 /**
173 * Encapsulates the results returned by {@link BesselJ#rjBesl(double, double, int)}.
174 * <p>
175 * {@link #getVals()} returns the computed function values.
176 * {@link #getnVals()} is the number of values among those returned by {@link #getnVals()}
177 * that can be considered accurate.
178 * </p><p>
179 * <ul>
180 * <li>nVals < 0: An argument is out of range. For example, nb <= 0, alpha
181 * < 0 or > 1, or x is too large. In this case, b(0) is set to zero, the
182 * remainder of the b-vector is not calculated, and nVals is set to
183 * MIN(nb,0) - 1 so that nVals != nb.</li>
184 * <li>nb > nVals > 0: Not all requested function values could be calculated
185 * accurately. This usually occurs because nb is much larger than abs(x). In
186 * this case, b(n) is calculated to the desired accuracy for n < nVals, but
187 * precision is lost for nVals < n <= nb. If b(n) does not vanish for n >
188 * nVals (because it is too small to be represented), and b(n)/b(nVals) =
189 * \(10^{-k}\), then only the first NSIG-k significant figures of b(n) can be
190 * trusted.</li></ul></p>
191 */
192 public static class BesselJResult {
193
194 /** Bessel function values */
195 private final double[] vals;
196
197 /** Valid value count */
198 private final int nVals;
199
200 /**
201 * Create a new BesselJResult with the given values and valid value count.
202 *
203 * @param b values
204 * @param n count of valid values
205 */
206 public BesselJResult(double[] b, int n) {
207 vals = MathArrays.copyOf(b, b.length);
208 nVals = n;
209 }
210
211 /**
212 * @return the computed function values
213 */
214 public double[] getVals() {
215 return MathArrays.copyOf(vals, vals.length);
216 }
217
218 /**
219 * @return the number of valid function values (normally the same as the
220 * length of the array returned by {@link #getnVals()})
221 */
222 public int getnVals() {
223 return nVals;
224 }
225 }
226
227 /**
228 * Calculates Bessel functions \(J_{n+alpha}(x)\) for
229 * non-negative argument x, and non-negative order n + alpha.
230 * <p>
231 * Before using the output vector, the user should check that
232 * nVals = nb, i.e., all orders have been calculated to the desired accuracy.
233 * See BesselResult class javadoc for details on return values.
234 * </p>
235 * @param x non-negative real argument for which J's are to be calculated
236 * @param alpha fractional part of order for which J's or exponentially
237 * scaled J's (\(J\cdot e^{x}\)) are to be calculated. 0 <= alpha < 1.0.
238 * @param nb integer number of functions to be calculated, nb > 0. The first
239 * function calculated is of order alpha, and the last is of order
240 * nb - 1 + alpha.
241 * @return BesselJResult a vector of the functions
242 * \(J_{alpha}(x)\) through \(J_{nb-1+alpha}(x)\), or the corresponding exponentially
243 * scaled functions and an integer output variable indicating possible errors
244 */
245 public static BesselJResult rjBesl(double x, double alpha, int nb) {
246 final double[] b = new double[nb];
247
248 int ncalc = 0;
249 double alpem = 0;
250 double alp2em = 0;
251
252 // ---------------------------------------------------------------------
253 // Check for out of range arguments.
254 // ---------------------------------------------------------------------
255 final int magx = (int) x;
256 if ((nb > 0) && (x >= X_MIN) && (x <= X_MAX) && (alpha >= 0) &&
257 (alpha < 1)) {
258 // ---------------------------------------------------------------------
259 // Initialize result array to zero.
260 // ---------------------------------------------------------------------
261 ncalc = nb;
262 for (int i = 0; i < nb; ++i) {
263 b[i] = 0;
264 }
265
266 // ---------------------------------------------------------------------
267 // Branch to use 2-term ascending series for small X and asymptotic
268 // form for large X when NB is not too large.
269 // ---------------------------------------------------------------------
270 double tempa;
271 double tempb;
272 double tempc;
273 double tover;
274 if (x < RTNSIG) {
275 // ---------------------------------------------------------------------
276 // Two-term ascending series for small X.
277 // ---------------------------------------------------------------------
278 tempa = 1;
279 alpem = 1 + alpha;
280 double halfx = 0;
281 if (x > ENMTEN) {
282 halfx = 0.5 * x;
283 }
284 if (alpha != 0) {
285 tempa = FastMath.pow(halfx, alpha) /
286 (alpha * Gamma.gamma(alpha));
287 }
288 tempb = 0;
289 if (x + 1 > 1) {
290 tempb = -halfx * halfx;
291 }
292 b[0] = tempa + (tempa * tempb / alpem);
293 if ((x != 0) && (b[0] == 0)) {
294 ncalc = 0;
295 }
296 if (nb != 1) {
297 if (x <= 0) {
298 for (int n = 1; n < nb; ++n) {
299 b[n] = 0;
300 }
301 } else {
302 // ---------------------------------------------------------------------
303 // Calculate higher order functions.
304 // ---------------------------------------------------------------------
305 tempc = halfx;
306 tover = tempb != 0 ? ENMTEN / tempb : 2 * ENMTEN / x;
307 for (int n = 1; n < nb; ++n) {
308 tempa /= alpem;
309 alpem += 1;
310 tempa *= tempc;
311 if (tempa <= tover * alpem) {
312 tempa = 0;
313 }
314 b[n] = tempa + (tempa * tempb / alpem);
315 if ((b[n] == 0) && (ncalc > n)) {
316 ncalc = n;
317 }
318 }
319 }
320 }
321 } else if ((x > 25.0) && (nb <= magx + 1)) {
322 // ---------------------------------------------------------------------
323 // Asymptotic series for X > 25
324 // ---------------------------------------------------------------------
325 final double xc = FastMath.sqrt(PI2 / x);
326 final double mul = 0.125 / x;
327 final double xin = mul * mul;
328 int m = 0;
329 if (x >= 130.0) {
330 m = 4;
331 } else if (x >= 35.0) {
332 m = 8;
333 } else {
334 m = 11;
335 }
336
337 final double xm = 4.0 * m;
338 // ---------------------------------------------------------------------
339 // Argument reduction for SIN and COS routines.
340 // ---------------------------------------------------------------------
341 double t = (double) ((int) ((x / TWOPI) + 0.5));
342 final double z = x - t * TOWPI1 - t * TWOPI2 - (alpha + 0.5) / PI2;
343 double vsin = FastMath.sin(z);
344 double vcos = FastMath.cos(z);
345 double gnu = 2 * alpha;
346 double capq;
347 double capp;
348 double s;
349 double t1;
350 double xk;
351 for (int i = 1; i <= 2; i++) {
352 s = (xm - 1 - gnu) * (xm - 1 + gnu) * xin * 0.5;
353 t = (gnu - (xm - 3.0)) * (gnu + (xm - 3.0));
354 capp = (s * t) / FACT[2 * m];
355 t1 = (gnu - (xm + 1)) * (gnu + (xm + 1));
356 capq = (s * t1) / FACT[2 * m + 1];
357 xk = xm;
358 int k = 2 * m;
359 t1 = t;
360
361 for (int j = 2; j <= m; j++) {
362 xk -= 4.0;
363 s = (xk - 1 - gnu) * (xk - 1 + gnu);
364 t = (gnu - (xk - 3.0)) * (gnu + (xk - 3.0));
365 capp = (capp + 1 / FACT[k - 2]) * s * t * xin;
366 capq = (capq + 1 / FACT[k - 1]) * s * t1 * xin;
367 k -= 2;
368 t1 = t;
369 }
370
371 capp += 1;
372 capq = (capq + 1) * ((gnu * gnu) - 1) * (0.125 / x);
373 b[i - 1] = xc * (capp * vcos - capq * vsin);
374 if (nb == 1) {
375 return new BesselJResult(MathArrays.copyOf(b, b.length),
376 ncalc);
377 }
378 t = vsin;
379 vsin = -vcos;
380 vcos = t;
381 gnu += 2.0;
382 }
383
384 // ---------------------------------------------------------------------
385 // If NB > 2, compute J(X,ORDER+I) I = 2, NB-1
386 // ---------------------------------------------------------------------
387 if (nb > 2) {
388 gnu = 2 * alpha + 2.0;
389 for (int j = 2; j < nb; ++j) {
390 b[j] = gnu * b[j - 1] / x - b[j - 2];
391 gnu += 2.0;
392 }
393 }
394 } else {
395 // ---------------------------------------------------------------------
396 // Use recurrence to generate results. First initialize the
397 // calculation of P*S.
398 // ---------------------------------------------------------------------
399 final int nbmx = nb - magx;
400 int n = magx + 1;
401 int nstart = 0;
402 int nend = 0;
403 double en = 2 * (n + alpha);
404 double plast = 1;
405 double p = en / x;
406 double pold;
407 // ---------------------------------------------------------------------
408 // Calculate general significance test.
409 // ---------------------------------------------------------------------
410 double test = 2 * ENSIG;
411 boolean readyToInitialize = false;
412 if (nbmx >= 3) {
413 // ---------------------------------------------------------------------
414 // Calculate P*S until N = NB-1. Check for possible
415 // overflow.
416 // ---------------------------------------------------------------------
417 tover = ENTEN / ENSIG;
418 nstart = magx + 2;
419 nend = nb - 1;
420 en = 2 * (nstart - 1 + alpha);
421 double psave;
422 double psavel;
423 for (int k = nstart; k <= nend; k++) {
424 n = k;
425 en += 2.0;
426 pold = plast;
427 plast = p;
428 p = (en * plast / x) - pold;
429 if (p > tover) {
430 // ---------------------------------------------------------------------
431 // To avoid overflow, divide P*S by TOVER. Calculate
432 // P*S until
433 // ABS(P) > 1.
434 // ---------------------------------------------------------------------
435 tover = ENTEN;
436 p /= tover;
437 plast /= tover;
438 psave = p;
439 psavel = plast;
440 nstart = n + 1;
441 do {
442 n += 1;
443 en += 2.0;
444 pold = plast;
445 plast = p;
446 p = (en * plast / x) - pold;
447 } while (p <= 1);
448 tempb = en / x;
449 // ---------------------------------------------------------------------
450 // Calculate backward test and find NCALC, the
451 // highest N such that
452 // the test is passed.
453 // ---------------------------------------------------------------------
454 test = pold * plast * (0.5 - 0.5 / (tempb * tempb));
455 test /= ENSIG;
456 p = plast * tover;
457 n -= 1;
458 en -= 2.0;
459 nend = FastMath.min(nb, n);
460 for (int l = nstart; l <= nend; l++) {
461 pold = psavel;
462 psavel = psave;
463 psave = (en * psavel / x) - pold;
464 if (psave * psavel > test) {
465 ncalc = l - 1;
466 readyToInitialize = true;
467 break;
468 }
469 }
470 ncalc = nend;
471 readyToInitialize = true;
472 break;
473 }
474 }
475 if (!readyToInitialize) {
476 n = nend;
477 en = 2 * (n + alpha);
478 // ---------------------------------------------------------------------
479 // Calculate special significance test for NBMX > 2.
480 // ---------------------------------------------------------------------
481 test = FastMath.max(test, FastMath.sqrt(plast * ENSIG) *
482 FastMath.sqrt(2 * p));
483 }
484 }
485 // ---------------------------------------------------------------------
486 // Calculate P*S until significance test passes.
487 // ---------------------------------------------------------------------
488 if (!readyToInitialize) {
489 do {
490 n += 1;
491 en += 2.0;
492 pold = plast;
493 plast = p;
494 p = (en * plast / x) - pold;
495 } while (p < test);
496 }
497 // ---------------------------------------------------------------------
498 // Initialize the backward recursion and the normalization sum.
499 // ---------------------------------------------------------------------
500 n += 1;
501 en += 2.0;
502 tempb = 0;
503 tempa = 1 / p;
504 int m = (2 * n) - 4 * (n / 2);
505 double sum = 0;
506 double em = (double) (n / 2);
507 alpem = em - 1 + alpha;
508 alp2em = 2 * em + alpha;
509 if (m != 0) {
510 sum = tempa * alpem * alp2em / em;
511 }
512 nend = n - nb;
513
514 boolean readyToNormalize = false;
515 boolean calculatedB0 = false;
516
517 // ---------------------------------------------------------------------
518 // Recur backward via difference equation, calculating (but not
519 // storing) B(N), until N = NB.
520 // ---------------------------------------------------------------------
521 for (int l = 1; l <= nend; l++) {
522 n -= 1;
523 en -= 2.0;
524 tempc = tempb;
525 tempb = tempa;
526 tempa = (en * tempb / x) - tempc;
527 m = 2 - m;
528 if (m != 0) {
529 em -= 1;
530 alp2em = 2 * em + alpha;
531 if (n == 1) {
532 break;
533 }
534 alpem = em - 1 + alpha;
535 if (alpem == 0) {
536 alpem = 1;
537 }
538 sum = (sum + tempa * alp2em) * alpem / em;
539 }
540 }
541
542 // ---------------------------------------------------------------------
543 // Store B(NB).
544 // ---------------------------------------------------------------------
545 b[n - 1] = tempa;
546 if (nend >= 0) {
547 if (nb <= 1) {
548 alp2em = alpha;
549 if (alpha + 1 == 1) {
550 alp2em = 1;
551 }
552 sum += b[0] * alp2em;
553 readyToNormalize = true;
554 } else {
555 // ---------------------------------------------------------------------
556 // Calculate and store B(NB-1).
557 // ---------------------------------------------------------------------
558 n -= 1;
559 en -= 2.0;
560 b[n - 1] = (en * tempa / x) - tempb;
561 if (n == 1) {
562 calculatedB0 = true;
563 } else {
564 m = 2 - m;
565 if (m != 0) {
566 em -= 1;
567 alp2em = 2 * em + alpha;
568 alpem = em - 1 + alpha;
569 if (alpem == 0) {
570 alpem = 1;
571 }
572
573 sum = (sum + (b[n - 1] * alp2em)) * alpem / em;
574 }
575 }
576 }
577 }
578 if (!readyToNormalize && !calculatedB0) {
579 nend = n - 2;
580 if (nend != 0) {
581 // ---------------------------------------------------------------------
582 // Calculate via difference equation and store B(N),
583 // until N = 2.
584 // ---------------------------------------------------------------------
585
586 for (int l = 1; l <= nend; l++) {
587 n -= 1;
588 en -= 2.0;
589 b[n - 1] = (en * b[n] / x) - b[n + 1];
590 m = 2 - m;
591 if (m != 0) {
592 em -= 1;
593 alp2em = 2 * em + alpha;
594 alpem = em - 1 + alpha;
595 if (alpem == 0) {
596 alpem = 1;
597 }
598
599 sum = (sum + b[n - 1] * alp2em) * alpem / em;
600 }
601 }
602 }
603 }
604 // ---------------------------------------------------------------------
605 // Calculate b[0]
606 // ---------------------------------------------------------------------
607 if (!readyToNormalize) {
608 if (!calculatedB0) {
609 b[0] = 2.0 * (alpha + 1) * b[1] / x - b[2];
610 }
611 em -= 1;
612 alp2em = 2 * em + alpha;
613 if (alp2em == 0) {
614 alp2em = 1;
615 }
616 sum += b[0] * alp2em;
617 }
618 // ---------------------------------------------------------------------
619 // Normalize. Divide all B(N) by sum.
620 // ---------------------------------------------------------------------
621
622 if (FastMath.abs(alpha) > 1e-16) {
623 sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha);
624 }
625 tempa = ENMTEN;
626 if (sum > 1) {
627 tempa *= sum;
628 }
629
630 for (n = 0; n < nb; n++) {
631 if (FastMath.abs(b[n]) < tempa) {
632 b[n] = 0;
633 }
634 b[n] /= sum;
635 }
636 }
637 // ---------------------------------------------------------------------
638 // Error return -- X, NB, or ALPHA is out of range.
639 // ---------------------------------------------------------------------
640 } else {
641 if (b.length > 0) {
642 b[0] = 0;
643 }
644 ncalc = FastMath.min(nb, 0) - 1;
645 }
646 return new BesselJResult(MathArrays.copyOf(b, b.length), ncalc);
647 }
648}
Note: See TracBrowser for help on using the repository browser.