?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Compilation: kill warning.
[gnumeric.git] / src / sf-bessel.c
blob2494a3edef41f1499157b99e23da448d4a2556da
1 #include <gnumeric-config.h>
2 #include <sf-bessel.h>
3 #include <sf-gamma.h>
4 #include <sf-trig.h>
5 #include <mathfunc.h>
6
7 #define MATHLIB_STANDALONE
8 #define ML_ERROR(cause) do { } while(0)
9 #define MATHLIB_ERROR(_a,_b) return gnm_nan;
10 #define M_SQRT_2dPI GNM_const(0.797884560802865355879892119869) /* sqrt(2/pi) */
11 #define MATHLIB_WARNING2 g_warning
12 #define MATHLIB_WARNING4 g_warning
13
14 static gnm_float bessel_k(gnm_float x, gnm_float alpha, gnm_float expo);
15
16 static inline int imin2 (int x, int y) { return MIN (x, y); }
17 static inline int imax2 (int x, int y) { return MAX (x, y); }
18 static inline gnm_float fmax2 (gnm_float x, gnm_float y) { return MAX (x, y); }
19
20 static gnm_float
21 gnm_cbrt (gnm_float x)
22 {
23 gnm_float ar = gnm_pow (gnm_abs (x), 1.0 / 3.0);
24 return x < 0 ? 0 - ar : ar;
25 }
26
27
28 /* ------------------------------------------------------------------------- */
29
30 /* Imported src/nmath/bessel.h from R. */
31
32 /* Constants und Documentation that apply to several of the
33 * ./bessel_[ijky].c files */
34
35 /* *******************************************************************
36
37 Explanation of machine-dependent constants
38
39 beta = Radix for the floating-point system
40 minexp = Smallest representable power of beta
41 maxexp = Smallest power of beta that overflows
42 it = p = Number of bits (base-beta digits) in the mantissa
43 (significand) of a working precision (floating-point) variable
44 NSIG = Decimal significance desired. Should be set to
45 INT(LOG10(2)*it+1). Setting NSIG lower will result
46 in decreased accuracy while setting NSIG higher will
47 increase CPU time without increasing accuracy. The
48 truncation error is limited to a relative error of
49 T=.5*10^(-NSIG).
50 ENTEN = 10 ^ K, where K is the largest long such that
51 ENTEN is machine-representable in working precision
52 ENSIG = 10 ^ NSIG
53 RTNSIG = 10 ^ (-K) for the smallest long K such that
54 K >= NSIG/4
55 ENMTEN = Smallest ABS(X) such that X/4 does not underflow
56 XINF = Largest positive machine number; approximately beta ^ maxexp
57 == GNM_MAX (defined in #include <float.h>)
58 SQXMIN = Square root of beta ^ minexp = gnm_sqrt(GNM_MIN)
59
60 EPS = The smallest positive floating-point number such that 1.0+EPS > 1.0
61 = beta ^ (-p) == GNM_EPSILON
62
63
64 For I :
65
66 EXPARG = Largest working precision argument that the library
67 EXP routine can handle and upper limit on the
68 magnitude of X when IZE=1; approximately LOG(beta ^ maxexp)
69
70 For I and J :
71
72 xlrg_IJ = (was = XLARGE). Upper limit on the magnitude of X (when
73 IZE=2 for I()). Bear in mind that if ABS(X)=N, then at least
74 N iterations of the backward recursion will be executed.
75 The value of 10 ^ 4 is used on every machine.
76
77 For j :
78 XMIN_J = Smallest acceptable argument for RBESY; approximately
79 max(2*beta ^ minexp, 2/XINF), rounded up
80
81 For Y :
82
83 xlrg_Y = (was = XLARGE). Upper bound on X;
84 approximately 1/DEL, because the sine and cosine functions
85 have lost about half of their precision at that point.
86
87 EPS_SINC = Machine number below which gnm_sin(x)/x = 1; approximately SQRT(EPS).
88 THRESH = Lower bound for use of the asymptotic form;
89 approximately AINT(-LOG10(EPS/2.0))+1.0
90
91
92 For K :
93
94 xmax_k = (was = XMAX). Upper limit on the magnitude of X when ize = 1;
95 i.e. maximal x for UNscaled answer.
96
97 Solution to equation:
98 W(X) * (1 -1/8 X + 9/128 X^2) = beta ^ minexp
99 where W(X) = EXP(-X)*SQRT(PI/2X)
100
101 --------------------------------------------------------------------
102
103 Approximate values for some important machines are:
104
105 beta minexp maxexp it NSIG ENTEN ENSIG RTNSIG ENMTEN EXPARG
106 IEEE (IBM/XT,
107 SUN, etc.) (S.P.) 2 -126 128 24 8 1e38 1e8 1e-2 4.70e-38 88
108 IEEE (...) (D.P.) 2 -1022 1024 53 16 1e308 1e16 1e-4 8.90e-308 709
109 CRAY-1 (S.P.) 2 -8193 8191 48 15 1e2465 1e15 1e-4 1.84e-2466 5677
110 Cyber 180/855
111 under NOS (S.P.) 2 -975 1070 48 15 1e322 1e15 1e-4 1.25e-293 741
112 IBM 3033 (D.P.) 16 -65 63 14 5 1e75 1e5 1e-2 2.16e-78 174
113 VAX (S.P.) 2 -128 127 24 8 1e38 1e8 1e-2 1.17e-38 88
114 VAX D-Format (D.P.) 2 -128 127 56 17 1e38 1e17 1e-5 1.17e-38 88
115 VAX G-Format (D.P.) 2 -1024 1023 53 16 1e307 1e16 1e-4 2.22e-308 709
116
117
118 And routine specific :
119
120 xlrg_IJ xlrg_Y xmax_k EPS_SINC XMIN_J XINF THRESH
121 IEEE (IBM/XT,
122 SUN, etc.) (S.P.) 1e4 1e4 85.337 1e-4 2.36e-38 3.40e38 8.
123 IEEE (...) (D.P.) 1e4 1e8 705.342 1e-8 4.46e-308 1.79e308 16.
124 CRAY-1 (S.P.) 1e4 2e7 5674.858 5e-8 3.67e-2466 5.45e2465 15.
125 Cyber 180/855
126 under NOS (S.P.) 1e4 2e7 672.788 5e-8 6.28e-294 1.26e322 15.
127 IBM 3033 (D.P.) 1e4 1e8 177.852 1e-8 2.77e-76 7.23e75 17.
128 VAX (S.P.) 1e4 1e4 86.715 1e-4 1.18e-38 1.70e38 8.
129 VAX e-Format (D.P.) 1e4 1e9 86.715 1e-9 1.18e-38 1.70e38 17.
130 VAX G-Format (D.P.) 1e4 1e8 706.728 1e-8 2.23e-308 8.98e307 16.
131
132 */
133 #define nsig_BESS 16
134 #define ensig_BESS 1e16
135 #define rtnsig_BESS 1e-4
136 #define enmten_BESS 8.9e-308
137 #define enten_BESS 1e308
138
139 #define exparg_BESS 709.
140 #define xlrg_BESS_IJ 1e5
141 #define xlrg_BESS_Y 1e8
142 #define thresh_BESS_Y 16.
143
144 #define xmax_BESS_K 705.342/* maximal x for UNscaled answer */
145
146
147 /* gnm_sqrt(GNM_MIN) = GNM_const(1.491668e-154) */
148 #define sqxmin_BESS_K 1.49e-154
149
150 /* x < eps_sinc <==> gnm_sin(x)/x == 1 (particularly "==>");
151 Linux (around 2001-02) gives GNM_const(2.14946906753213e-08)
152 Solaris 2.5.1 gives GNM_const(2.14911933289084e-08)
153 */
154 #define M_eps_sinc 2.149e-8
155
156 /* ------------------------------------------------------------------------ */
157 /* Imported src/nmath/bessel_i.c from R. */
158 /*
159 * Mathlib : A C Library of Special Functions
160 * Copyright (C) 1998-2001 Ross Ihaka and the R Development Core team.
161 *
162 * This program is free software; you can redistribute it and/or modify
163 * it under the terms of the GNU General Public License as published by
164 * the Free Software Foundation; either version 2 of the License, or
165 * (at your option) any later version.
166 *
167 * This program is distributed in the hope that it will be useful,
168 * but WITHOUT ANY WARRANTY; without even the implied warranty of
169 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
170 * GNU General Public License for more details.
171 *
172 * You should have received a copy of the GNU General Public License
173 * along with this program; if not, write to the Free Software
174 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
175 * USA.
176 */
177
178 /* DESCRIPTION --> see below */
179
180
181 /* From http://www.netlib.org/specfun/ribesl Fortran translated by f2c,...
182 * ------------------------------=#---- Martin Maechler, ETH Zurich
183 */
184
185 #ifndef MATHLIB_STANDALONE
186 #endif
187
188 static void I_bessel(gnm_float *x, gnm_float *alpha, long *nb,
189 long *ize, gnm_float *bi, long *ncalc);
190
191 static gnm_float bessel_i(gnm_float x, gnm_float alpha, gnm_float expo)
192 {
193 long nb, ncalc, ize;
194 gnm_float *bi;
195 #ifndef MATHLIB_STANDALONE
196 char *vmax;
197 #endif
198
199 #ifdef IEEE_754
200 /* NaNs propagated correctly */
201 if (gnm_isnan(x) || gnm_isnan(alpha)) return x + alpha;
202 #endif
203 if (x < 0) {
204 ML_ERROR(ME_RANGE);
205 return gnm_nan;
206 }
207 ize = (long)expo;
208 if (alpha < 0) {
209 /* Using Abramowitz & Stegun 9.6.2
210 * this may not be quite optimal (CPU and accuracy wise) */
211 return(bessel_i(x, -alpha, expo) +
212 bessel_k(x, -alpha, expo) * ((ize == 1)? 2. : 2.*gnm_exp(-x))/M_PIgnum
213 * gnm_sinpi(-alpha));
214 }
215 nb = 1+ (long)gnm_floor(alpha);/* nb-1 <= alpha < nb */
216 alpha -= (nb-1);
217 #ifdef MATHLIB_STANDALONE
218 bi = (gnm_float *) calloc(nb, sizeof(gnm_float));
219 if (!bi) MATHLIB_ERROR("%s", ("bessel_i allocation error"));
220 #else
221 vmax = vmaxget();
222 bi = (gnm_float *) R_alloc(nb, sizeof(gnm_float));
223 #endif
224 I_bessel(&x, &alpha, &nb, &ize, bi, &ncalc);
225 if(ncalc != nb) {/* error input */
226 if(ncalc < 0)
227 MATHLIB_WARNING4(("bessel_i(%" GNM_FORMAT_g "): ncalc (=%ld) != nb (=%ld); alpha=%" GNM_FORMAT_g ". Arg. out of range?\n"),
228 x, ncalc, nb, alpha);
229 else
230 MATHLIB_WARNING2(("bessel_i(%" GNM_FORMAT_g ",nu=%" GNM_FORMAT_g "): precision lost in result\n"),
231 x, alpha+nb-1);
232 }
233 x = bi[nb-1];
234 #ifdef MATHLIB_STANDALONE
235 free(bi);
236 #else
237 vmaxset(vmax);
238 #endif
239 return x;
240 }
241
242 static void I_bessel(gnm_float *x, gnm_float *alpha, long *nb,
243 long *ize, gnm_float *bi, long *ncalc)
244 {
245 /* -------------------------------------------------------------------
246
247 This routine calculates Bessel functions I_(N+ALPHA) (X)
248 for non-negative argument X, and non-negative order N+ALPHA,
249 with or without exponential scaling.
250
251
252 Explanation of variables in the calling sequence
253
254 X - Non-negative argument for which
255 I's or exponentially scaled I's (I*EXP(-X))
256 are to be calculated. If I's are to be calculated,
257 X must be less than EXPARG_BESS (see bessel.h).
258 ALPHA - Fractional part of order for which
259 I's or exponentially scaled I's (I*EXP(-X)) are
260 to be calculated. 0 <= ALPHA < 1.0.
261 NB - Number of functions to be calculated, NB > 0.
262 The first function calculated is of order ALPHA, and the
263 last is of order (NB - 1 + ALPHA).
264 IZE - Type. IZE = 1 if unscaled I's are to be calculated,
265 = 2 if exponentially scaled I's are to be calculated.
266 BI - Output vector of length NB. If the routine
267 terminates normally (NCALC=NB), the vector BI contains the
268 functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the
269 corresponding exponentially scaled functions.
270 NCALC - Output variable indicating possible errors.
271 Before using the vector BI, the user should check that
272 NCALC=NB, i.e., all orders have been calculated to
273 the desired accuracy. See error returns below.
274
275
276 *******************************************************************
277 *******************************************************************
278
279 Error returns
280
281 In case of an error, NCALC != NB, and not all I's are
282 calculated to the desired accuracy.
283
284 NCALC < 0: An argument is out of range. For example,
285 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= EXPARG_BESS.
286 In this case, the BI-vector is not calculated, and NCALC is
287 set to MIN0(NB,0)-1 so that NCALC != NB.
288
289 NB > NCALC > 0: Not all requested function values could
290 be calculated accurately. This usually occurs because NB is
291 much larger than ABS(X). In this case, BI[N] is calculated
292 to the desired accuracy for N <= NCALC, but precision
293 is lost for NCALC < N <= NB. If BI[N] does not vanish
294 for N > NCALC (because it is too small to be represented),
295 and BI[N]/BI[NCALC] = 10**(-K), then only the first NSIG-K
296 significant figures of BI[N] can be trusted.
297
298
299 Intrinsic functions required are:
300
301 DBLE, EXP, gamma_cody, INT, MAX, MIN, REAL, SQRT
302
303
304 Acknowledgement
305
306 This program is based on a program written by David J.
307 Sookne (2) that computes values of the Bessel functions J or
308 I of float argument and long order. Modifications include
309 the restriction of the computation to the I Bessel function
310 of non-negative float argument, the extension of the computation
311 to arbitrary positive order, the inclusion of optional
312 exponential scaling, and the elimination of most underflow.
313 An earlier version was published in (3).
314
315 References: "A Note on Backward Recurrence Algorithms," Olver,
316 F. W. J., and Sookne, D. J., Math. Comp. 26, 1972,
317 pp 941-947.
318
319 "Bessel Functions of Real Argument and Integer Order,"
320 Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp
321 125-132.
322
323 "ALGORITHM 597, Sequence of Modified Bessel Functions
324 of the First Kind," Cody, W. J., Trans. Math. Soft.,
325 1983, pp. 242-245.
326
327 Latest modification: May 30, 1989
328
329 Modified by: W. J. Cody and L. Stoltz
330 Applied Mathematics Division
331 Argonne National Laboratory
332 Argonne, IL 60439
333 */
334
335 /*-------------------------------------------------------------------
336 Mathematical constants
337 -------------------------------------------------------------------*/
338 const gnm_float const__ = 1.585;
339
340 /* Local variables */
341 long nend, intx, nbmx, k, l, n, nstart;
342 gnm_float pold, test, p, em, en, empal, emp2al, halfx,
343 aa, bb, cc, psave, plast, tover, psavel, sum, nu, twonu;
344
345 /*Parameter adjustments */
346 --bi;
347 nu = *alpha;
348 twonu = nu + nu;
349
350 /*-------------------------------------------------------------------
351 Check for X, NB, OR IZE out of range.
352 ------------------------------------------------------------------- */
353 if (*nb > 0 && *x >= 0. && (0. <= nu && nu < 1.) &&
354 (1 <= *ize && *ize <= 2) ) {
355
356 *ncalc = *nb;
357 if((*ize == 1 && *x > exparg_BESS) ||
358 (*ize == 2 && *x > xlrg_BESS_IJ)) {
359 ML_ERROR(ME_RANGE);
360 for(k=1; k <= *nb; k++)
361 bi[k]=gnm_pinf;
362 return;
363 }
364 intx = (long) (*x);/* --> we will probably fail when *x > LONG_MAX */
365 if (*x >= rtnsig_BESS) { /* "non-small" x */
366 /* -------------------------------------------------------------------
367 Initialize the forward sweep, the P-sequence of Olver
368 ------------------------------------------------------------------- */
369 nbmx = *nb - intx;
370 n = intx + 1;
371 en = (gnm_float) (n + n) + twonu;
372 plast = 1.;
373 p = en / *x;
374 /* ------------------------------------------------
375 Calculate general significance test
376 ------------------------------------------------ */
377 test = ensig_BESS + ensig_BESS;
378 if (intx << 1 > nsig_BESS * 5) {
379 test = gnm_sqrt(test * p);
380 } else {
381 test /= gnm_pow(const__, (gnm_float)intx);
382 }
383 if (nbmx >= 3) {
384 /* --------------------------------------------------
385 Calculate P-sequence until N = NB-1
386 Check for possible overflow.
387 ------------------------------------------------ */
388 tover = enten_BESS / ensig_BESS;
389 nstart = intx + 2;
390 nend = *nb - 1;
391 for (k = nstart; k <= nend; ++k) {
392 n = k;
393 en += 2.;
394 pold = plast;
395 plast = p;
396 p = en * plast / *x + pold;
397 if (p > tover) {
398 /* ------------------------------------------------
399 To avoid overflow, divide P-sequence by TOVER.
400 Calculate P-sequence until ABS(P) > 1.
401 ---------------------------------------------- */
402 tover = enten_BESS;
403 p /= tover;
404 plast /= tover;
405 psave = p;
406 psavel = plast;
407 nstart = n + 1;
408 do {
409 ++n;
410 en += 2.;
411 pold = plast;
412 plast = p;
413 p = en * plast / *x + pold;
414 }
415 while (p <= 1.);
416
417 bb = en / *x;
418 /* ------------------------------------------------
419 Calculate backward test, and find NCALC,
420 the highest N such that the test is passed.
421 ------------------------------------------------ */
422 test = pold * plast / ensig_BESS;
423 test *= .5 - .5 / (bb * bb);
424 p = plast * tover;
425 --n;
426 en -= 2.;
427 nend = imin2(*nb,n);
428 for (l = nstart; l <= nend; ++l) {
429 *ncalc = l;
430 pold = psavel;
431 psavel = psave;
432 psave = en * psavel / *x + pold;
433 if (psave * psavel > test) {
434 goto L90;
435 }
436 }
437 *ncalc = nend + 1;
438 L90:
439 --(*ncalc);
440 goto L120;
441 }
442 }
443 n = nend;
444 en = (gnm_float)(n + n) + twonu;
445 /*---------------------------------------------------
446 Calculate special significance test for NBMX > 2.
447 --------------------------------------------------- */
448 test = fmax2(test,gnm_sqrt(plast * ensig_BESS) * gnm_sqrt(p + p));
449 }
450 /* --------------------------------------------------------
451 Calculate P-sequence until significance test passed.
452 -------------------------------------------------------- */
453 do {
454 ++n;
455 en += 2.;
456 pold = plast;
457 plast = p;
458 p = en * plast / *x + pold;
459 } while (p < test);
460
461 L120:
462 /* -------------------------------------------------------------------
463 Initialize the backward recursion and the normalization sum.
464 ------------------------------------------------------------------- */
465 ++n;
466 en += 2.;
467 bb = 0.;
468 aa = 1. / p;
469 em = (gnm_float) n - 1.;
470 empal = em + nu;
471 emp2al = em - 1. + twonu;
472 sum = aa * empal * emp2al / em;
473 nend = n - *nb;
474 if (nend < 0) {
475 /* -----------------------------------------------------
476 N < NB, so store BI[N] and set higher orders to 0..
477 ----------------------------------------------------- */
478 bi[n] = aa;
479 nend = -nend;
480 for (l = 1; l <= nend; ++l) {
481 bi[n + l] = 0.;
482 }
483 } else {
484 if (nend > 0) {
485 /* -----------------------------------------------------
486 Recur backward via difference equation,
487 calculating (but not storing) BI[N], until N = NB.
488 --------------------------------------------------- */
489 for (l = 1; l <= nend; ++l) {
490 --n;
491 en -= 2.;
492 cc = bb;
493 bb = aa;
494 aa = en * bb / *x + cc;
495 em -= 1.;
496 emp2al -= 1.;
497 if (n == 1) {
498 break;
499 }
500 if (n == 2) {
501 emp2al = 1.;
502 }
503 empal -= 1.;
504 sum = (sum + aa * empal) * emp2al / em;
505 }
506 }
507 /* ---------------------------------------------------
508 Store BI[NB]
509 --------------------------------------------------- */
510 bi[n] = aa;
511 if (*nb <= 1) {
512 sum = sum + sum + aa;
513 goto L230;
514 }
515 /* -------------------------------------------------
516 Calculate and Store BI[NB-1]
517 ------------------------------------------------- */
518 --n;
519 en -= 2.;
520 bi[n] = en * aa / *x + bb;
521 if (n == 1) {
522 goto L220;
523 }
524 em -= 1.;
525 if (n == 2)
526 emp2al = 1.;
527 else
528 emp2al -= 1.;
529
530 empal -= 1.;
531 sum = (sum + bi[n] * empal) * emp2al / em;
532 }
533 nend = n - 2;
534 if (nend > 0) {
535 /* --------------------------------------------
536 Calculate via difference equation
537 and store BI[N], until N = 2.
538 ------------------------------------------ */
539 for (l = 1; l <= nend; ++l) {
540 --n;
541 en -= 2.;
542 bi[n] = en * bi[n + 1] / *x + bi[n + 2];
543 em -= 1.;
544 if (n == 2)
545 emp2al = 1.;
546 else
547 emp2al -= 1.;
548 empal -= 1.;
549 sum = (sum + bi[n] * empal) * emp2al / em;
550 }
551 }
552 /* ----------------------------------------------
553 Calculate BI[1]
554 -------------------------------------------- */
555 bi[1] = 2. * empal * bi[2] / *x + bi[3];
556 L220:
557 sum = sum + sum + bi[1];
558
559 L230:
560 /* ---------------------------------------------------------
561 Normalize. Divide all BI[N] by sum.
562 --------------------------------------------------------- */
563 if (nu != 0.)
564 sum *= (gnm_exp(lgamma1p (nu)) * gnm_pow(*x * .5, -nu));
565 if (*ize == 1)
566 sum *= gnm_exp(-(*x));
567 aa = enmten_BESS;
568 if (sum > 1.)
569 aa *= sum;
570 for (n = 1; n <= *nb; ++n) {
571 if (bi[n] < aa)
572 bi[n] = 0.;
573 else
574 bi[n] /= sum;
575 }
576 return;
577 } else {
578 /* -----------------------------------------------------------
579 Two-term ascending series for small X.
580 -----------------------------------------------------------*/
581 aa = 1.;
582 empal = 1. + nu;
583 if (*x > enmten_BESS)
584 halfx = .5 * *x;
585 else
586 halfx = 0.;
587 if (nu != 0.)
588 aa = gnm_pow(halfx, nu) / gnm_exp(lgamma1p(nu));
589 if (*ize == 2)
590 aa *= gnm_exp(-(*x));
591 if (*x + 1. > 1.)
592 bb = halfx * halfx;
593 else
594 bb = 0.;
595
596 bi[1] = aa + aa * bb / empal;
597 if (*x != 0. && bi[1] == 0.)
598 *ncalc = 0;
599 if (*nb > 1) {
600 if (*x == 0.) {
601 for (n = 2; n <= *nb; ++n) {
602 bi[n] = 0.;
603 }
604 } else {
605 /* -------------------------------------------------
606 Calculate higher-order functions.
607 ------------------------------------------------- */
608 cc = halfx;
609 tover = (enmten_BESS + enmten_BESS) / *x;
610 if (bb != 0.)
611 tover = enmten_BESS / bb;
612 for (n = 2; n <= *nb; ++n) {
613 aa /= empal;
614 empal += 1.;
615 aa *= cc;
616 if (aa <= tover * empal)
617 bi[n] = aa = 0.;
618 else
619 bi[n] = aa + aa * bb / empal;
620 if (bi[n] == 0. && *ncalc > n)
621 *ncalc = n - 1;
622 }
623 }
624 }
625 }
626 } else {
627 *ncalc = imin2(*nb,0) - 1;
628 }
629 }
630
631 /* ------------------------------------------------------------------------ */
632 /* Imported src/nmath/bessel_k.c from R. */
633 /*
634 * Mathlib : A C Library of Special Functions
635 * Copyright (C) 1998-2001 Ross Ihaka and the R Development Core team.
636 * Copyright (C) 2002-3 The R Foundation
637 *
638 * This program is free software; you can redistribute it and/or modify
639 * it under the terms of the GNU General Public License as published by
640 * the Free Software Foundation; either version 2 of the License, or
641 * (at your option) any later version.
642 *
643 * This program is distributed in the hope that it will be useful,
644 * but WITHOUT ANY WARRANTY; without even the implied warranty of
645 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
646 * GNU General Public License for more details.
647 *
648 * You should have received a copy of the GNU General Public License
649 * along with this program; if not, write to the Free Software
650 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
651 * USA.
652 */
653
654 /* DESCRIPTION --> see below */
655
656
657 /* From http://www.netlib.org/specfun/rkbesl Fortran translated by f2c,...
658 * ------------------------------=#---- Martin Maechler, ETH Zurich
659 */
660
661 #ifndef MATHLIB_STANDALONE
662 #endif
663
664 static void K_bessel(gnm_float *x, gnm_float *alpha, long *nb,
665 long *ize, gnm_float *bk, long *ncalc);
666
667 static gnm_float bessel_k(gnm_float x, gnm_float alpha, gnm_float expo)
668 {
669 long nb, ncalc, ize;
670 gnm_float *bk;
671 #ifndef MATHLIB_STANDALONE
672 char *vmax;
673 #endif
674
675 #ifdef IEEE_754
676 /* NaNs propagated correctly */
677 if (gnm_isnan(x) || gnm_isnan(alpha)) return x + alpha;
678 #endif
679 if (x < 0) {
680 ML_ERROR(ME_RANGE);
681 return gnm_nan;
682 }
683 ize = (long)expo;
684 if(alpha < 0)
685 alpha = -alpha;
686 nb = 1+ (long)gnm_floor(alpha);/* nb-1 <= |alpha| < nb */
687 alpha -= (nb-1);
688 #ifdef MATHLIB_STANDALONE
689 bk = (gnm_float *) calloc(nb, sizeof(gnm_float));
690 if (!bk) MATHLIB_ERROR("%s", ("bessel_k allocation error"));
691 #else
692 vmax = vmaxget();
693 bk = (gnm_float *) R_alloc(nb, sizeof(gnm_float));
694 #endif
695 K_bessel(&x, &alpha, &nb, &ize, bk, &ncalc);
696 if(ncalc != nb) {/* error input */
697 if(ncalc < 0)
698 MATHLIB_WARNING4(("bessel_k(%" GNM_FORMAT_g "): ncalc (=%ld) != nb (=%ld); alpha=%" GNM_FORMAT_g ". Arg. out of range?\n"),
699 x, ncalc, nb, alpha);
700 else
701 MATHLIB_WARNING2(("bessel_k(%" GNM_FORMAT_g ",nu=%" GNM_FORMAT_g "): precision lost in result\n"),
702 x, alpha+nb-1);
703 }
704 x = bk[nb-1];
705 #ifdef MATHLIB_STANDALONE
706 free(bk);
707 #else
708 vmaxset(vmax);
709 #endif
710 return x;
711 }
712
713 static void K_bessel(gnm_float *x, gnm_float *alpha, long *nb,
714 long *ize, gnm_float *bk, long *ncalc)
715 {
716 /*-------------------------------------------------------------------
717
718 This routine calculates modified Bessel functions
719 of the third kind, K_(N+ALPHA) (X), for non-negative
720 argument X, and non-negative order N+ALPHA, with or without
721 exponential scaling.
722
723 Explanation of variables in the calling sequence
724
725 X - Non-negative argument for which
726 K's or exponentially scaled K's (K*EXP(X))
727 are to be calculated. If K's are to be calculated,
728 X must not be greater than XMAX_BESS_K.
729 ALPHA - Fractional part of order for which
730 K's or exponentially scaled K's (K*EXP(X)) are
731 to be calculated. 0 <= ALPHA < 1.0.
732 NB - Number of functions to be calculated, NB > 0.
733 The first function calculated is of order ALPHA, and the
734 last is of order (NB - 1 + ALPHA).
735 IZE - Type. IZE = 1 if unscaled K's are to be calculated,
736 = 2 if exponentially scaled K's are to be calculated.
737 BK - Output vector of length NB. If the
738 routine terminates normally (NCALC=NB), the vector BK
739 contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X),
740 or the corresponding exponentially scaled functions.
741 If (0 < NCALC < NB), BK(I) contains correct function
742 values for I <= NCALC, and contains the ratios
743 K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
744 NCALC - Output variable indicating possible errors.
745 Before using the vector BK, the user should check that
746 NCALC=NB, i.e., all orders have been calculated to
747 the desired accuracy. See error returns below.
748
749
750 *******************************************************************
751
752 Error returns
753
754 In case of an error, NCALC != NB, and not all K's are
755 calculated to the desired accuracy.
756
757 NCALC < -1: An argument is out of range. For example,
758 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K.
759 In this case, the B-vector is not calculated,
760 and NCALC is set to MIN0(NB,0)-2 so that NCALC != NB.
761 NCALC = -1: Either K(ALPHA,X) >= XINF or
762 K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF. In this case,
763 the B-vector is not calculated. Note that again
764 NCALC != NB.
765
766 0 < NCALC < NB: Not all requested function values could
767 be calculated accurately. BK(I) contains correct function
768 values for I <= NCALC, and contains the ratios
769 K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
770
771
772 Intrinsic functions required are:
773
774 ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT
775
776
777 Acknowledgement
778
779 This program is based on a program written by J. B. Campbell
780 (2) that computes values of the Bessel functions K of float
781 argument and float order. Modifications include the addition
782 of non-scaled functions, parameterization of machine
783 dependencies, and the use of more accurate approximations
784 for SINH and SIN.
785
786 References: "On Temme's Algorithm for the Modified Bessel
787 Functions of the Third Kind," Campbell, J. B.,
788 TOMS 6(4), Dec. 1980, pp. 581-586.
789
790 "A FORTRAN IV Subroutine for the Modified Bessel
791 Functions of the Third Kind of Real Order and Real
792 Argument," Campbell, J. B., Report NRC/ERB-925,
793 National Research Council, Canada.
794
795 Latest modification: May 30, 1989
796
797 Modified by: W. J. Cody and L. Stoltz
798 Applied Mathematics Division
799 Argonne National Laboratory
800 Argonne, IL 60439
801
802 -------------------------------------------------------------------
803 */
804 /*---------------------------------------------------------------------
805 * Mathematical constants
806 * A = LOG(2) - Euler's constant
807 * D = SQRT(2/PI)
808 ---------------------------------------------------------------------*/
809 const gnm_float a = GNM_const(.11593151565841244881);
810
811 /*---------------------------------------------------------------------
812 P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant
813 Coefficients converted from hex to decimal and modified
814 by W. J. Cody, 2/26/82 */
815 static const gnm_float p[8] = { GNM_const(.805629875690432845),GNM_const(20.4045500205365151),
816 GNM_const(157.705605106676174),GNM_const(536.671116469207504),GNM_const(900.382759291288778),
817 GNM_const(730.923886650660393),GNM_const(229.299301509425145),GNM_const(.822467033424113231) };
818 static const gnm_float q[7] = { GNM_const(29.4601986247850434),GNM_const(277.577868510221208),
819 GNM_const(1206.70325591027438),GNM_const(2762.91444159791519),GNM_const(3443.74050506564618),
820 GNM_const(2210.63190113378647),GNM_const(572.267338359892221) };
821 /* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */
822 static const gnm_float r[5] = { GNM_const(-.48672575865218401848),GNM_const(13.079485869097804016),
823 GNM_const(-101.96490580880537526),GNM_const(347.65409106507813131),
824 GNM_const(3.495898124521934782e-4) };
825 static const gnm_float s[4] = { GNM_const(-25.579105509976461286),GNM_const(212.57260432226544008),
826 GNM_const(-610.69018684944109624),GNM_const(422.69668805777760407) };
827 /* T - Approximation for SINH(Y)/Y */
828 static const gnm_float t[6] = { GNM_const(1.6125990452916363814e-10),
829 GNM_const(2.5051878502858255354e-8),GNM_const(2.7557319615147964774e-6),
830 GNM_const(1.9841269840928373686e-4),GNM_const(.0083333333333334751799),
831 GNM_const(.16666666666666666446) };
832 /*---------------------------------------------------------------------*/
833 static const gnm_float estm[6] = { 52.0583,5.7607,2.7782,14.4303,185.3004, 9.3715 };
834 static const gnm_float estf[7] = { 41.8341,7.1075,6.4306,42.511,GNM_const(1.35633),84.5096,20.};
835
836 /* Local variables */
837 long iend, i, j, k, m, ii, mplus1;
838 gnm_float x2by4, twox, c, blpha, ratio, wminf;
839 gnm_float d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu;
840 gnm_float dm, ex, bk1, bk2, nu;
841
842 ii = 0; /* -Wall */
843
844 ex = *x;
845 nu = *alpha;
846 *ncalc = imin2(*nb,0) - 2;
847 if (*nb > 0 && (0. <= nu && nu < 1.) && (1 <= *ize && *ize <= 2)) {
848 if(ex <= 0 || (*ize == 1 && ex > xmax_BESS_K)) {
849 if(ex <= 0) {
850 ML_ERROR(ME_RANGE);
851 for(i=0; i < *nb; i++)
852 bk[i] = gnm_pinf;
853 } else /* would only have underflow */
854 for(i=0; i < *nb; i++)
855 bk[i] = 0.;
856 *ncalc = *nb;
857 return;
858 }
859 k = 0;
860 if (nu < sqxmin_BESS_K) {
861 nu = 0.;
862 } else if (nu > .5) {
863 k = 1;
864 nu -= 1.;
865 }
866 twonu = nu + nu;
867 iend = *nb + k - 1;
868 c = nu * nu;
869 d3 = -c;
870 if (ex <= 1.) {
871 /* ------------------------------------------------------------
872 Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA
873 Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA
874 ------------------------------------------------------------ */
875 d1 = 0.; d2 = p[0];
876 t1 = 1.; t2 = q[0];
877 for (i = 2; i <= 7; i += 2) {
878 d1 = c * d1 + p[i - 1];
879 d2 = c * d2 + p[i];
880 t1 = c * t1 + q[i - 1];
881 t2 = c * t2 + q[i];
882 }
883 d1 = nu * d1;
884 t1 = nu * t1;
885 f1 = gnm_log(ex);
886 f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1;
887 q0 = gnm_exp(-nu * (a - nu * (p[7] + nu * (d1-d2) / (t1-t2)) - f1));
888 f1 = nu * f0;
889 p0 = gnm_exp(f1);
890 /* -----------------------------------------------------------
891 Calculation of F0 =
892 ----------------------------------------------------------- */
893 d1 = r[4];
894 t1 = 1.;
895 for (i = 0; i < 4; ++i) {
896 d1 = c * d1 + r[i];
897 t1 = c * t1 + s[i];
898 }
899 /* d2 := gnm_sinh(f1)/ nu = gnm_sinh(f1)/(f1/f0)
900 * = f0 * gnm_sinh(f1)/f1 */
901 if (gnm_abs(f1) <= .5) {
902 f1 *= f1;
903 d2 = 0.;
904 for (i = 0; i < 6; ++i) {
905 d2 = f1 * d2 + t[i];
906 }
907 d2 = f0 + f0 * f1 * d2;
908 } else {
909 d2 = gnm_sinh(f1) / nu;
910 }
911 f0 = d2 - nu * d1 / (t1 * p0);
912 if (ex <= 1e-10) {
913 /* ---------------------------------------------------------
914 X <= 1.0E-10
915 Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X)
916 --------------------------------------------------------- */
917 bk[0] = f0 + ex * f0;
918 if (*ize == 1) {
919 bk[0] -= ex * bk[0];
920 }
921 ratio = p0 / f0;
922 c = ex * GNM_MAX;
923 if (k != 0) {
924 /* ---------------------------------------------------
925 Calculation of K(ALPHA,X)
926 and X*K(ALPHA+1,X)/K(ALPHA,X), ALPHA >= 1/2
927 --------------------------------------------------- */
928 *ncalc = -1;
929 if (bk[0] >= c / ratio) {
930 return;
931 }
932 bk[0] = ratio * bk[0] / ex;
933 twonu += 2.;
934 ratio = twonu;
935 }
936 *ncalc = 1;
937 if (*nb == 1)
938 return;
939
940 /* -----------------------------------------------------
941 Calculate K(ALPHA+L,X)/K(ALPHA+L-1,X),
942 L = 1, 2, ... , NB-1
943 ----------------------------------------------------- */
944 *ncalc = -1;
945 for (i = 1; i < *nb; ++i) {
946 if (ratio >= c)
947 return;
948
949 bk[i] = ratio / ex;
950 twonu += 2.;
951 ratio = twonu;
952 }
953 *ncalc = 1;
954 goto L420;
955 } else {
956 /* ------------------------------------------------------
957 10^-10 < X <= 1.0
958 ------------------------------------------------------ */
959 c = 1.;
960 x2by4 = ex * ex / 4.;
961 p0 = .5 * p0;
962 q0 = .5 * q0;
963 d1 = -1.;
964 d2 = 0.;
965 bk1 = 0.;
966 bk2 = 0.;
967 f1 = f0;
968 f2 = p0;
969 do {
970 d1 += 2.;
971 d2 += 1.;
972 d3 = d1 + d3;
973 c = x2by4 * c / d2;
974 f0 = (d2 * f0 + p0 + q0) / d3;
975 p0 /= d2 - nu;
976 q0 /= d2 + nu;
977 t1 = c * f0;
978 t2 = c * (p0 - d2 * f0);
979 bk1 += t1;
980 bk2 += t2;
981 } while (gnm_abs(t1 / (f1 + bk1)) > GNM_EPSILON ||
982 gnm_abs(t2 / (f2 + bk2)) > GNM_EPSILON);
983 bk1 = f1 + bk1;
984 bk2 = 2. * (f2 + bk2) / ex;
985 if (*ize == 2) {
986 d1 = gnm_exp(ex);
987 bk1 *= d1;
988 bk2 *= d1;
989 }
990 wminf = estf[0] * ex + estf[1];
991 }
992 } else if (GNM_EPSILON * ex > 1.) {
993 /* -------------------------------------------------
994 X > 1./EPS
995 ------------------------------------------------- */
996 *ncalc = *nb;
997 bk1 = 1. / (M_SQRT_2dPI * gnm_sqrt(ex));
998 for (i = 0; i < *nb; ++i)
999 bk[i] = bk1;
1000 return;
1001
1002 } else {
1003 /* -------------------------------------------------------
1004 X > 1.0
1005 ------------------------------------------------------- */
1006 twox = ex + ex;
1007 blpha = 0.;
1008 ratio = 0.;
1009 if (ex <= 4.) {
1010 /* ----------------------------------------------------------
1011 Calculation of K(ALPHA+1,X)/K(ALPHA,X), 1.0 <= X <= 4.0
1012 ----------------------------------------------------------*/
1013 d2 = gnm_trunc(estm[0] / ex + estm[1]);
1014 m = (long) d2;
1015 d1 = d2 + d2;
1016 d2 -= .5;
1017 d2 *= d2;
1018 for (i = 2; i <= m; ++i) {
1019 d1 -= 2.;
1020 d2 -= d1;
1021 ratio = (d3 + d2) / (twox + d1 - ratio);
1022 }
1023 /* -----------------------------------------------------------
1024 Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward
1025 recurrence and K(ALPHA,X) from the wronskian
1026 -----------------------------------------------------------*/
1027 d2 = gnm_trunc(estm[2] * ex + estm[3]);
1028 m = (long) d2;
1029 c = gnm_abs(nu);
1030 d3 = c + c;
1031 d1 = d3 - 1.;
1032 f1 = GNM_MIN;
1033 f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * GNM_MIN;
1034 for (i = 3; i <= m; ++i) {
1035 d2 -= 1.;
1036 f2 = (d3 + d2 + d2) * f0;
1037 blpha = (1. + d1 / d2) * (f2 + blpha);
1038 f2 = f2 / ex + f1;
1039 f1 = f0;
1040 f0 = f2;
1041 }
1042 f1 = (d3 + 2.) * f0 / ex + f1;
1043 d1 = 0.;
1044 t1 = 1.;
1045 for (i = 1; i <= 7; ++i) {
1046 d1 = c * d1 + p[i - 1];
1047 t1 = c * t1 + q[i - 1];
1048 }
1049 p0 = gnm_exp(c * (a + c * (p[7] - c * d1 / t1) - gnm_log(ex))) / ex;
1050 f2 = (c + .5 - ratio) * f1 / ex;
1051 bk1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0;
1052 if (*ize == 1) {
1053 bk1 *= gnm_exp(-ex);
1054 }
1055 wminf = estf[2] * ex + estf[3];
1056 } else {
1057 /* ---------------------------------------------------------
1058 Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by
1059 backward recurrence, for X > 4.0
1060 ----------------------------------------------------------*/
1061 dm = gnm_trunc(estm[4] / ex + estm[5]);
1062 m = (long) dm;
1063 d2 = dm - .5;
1064 d2 *= d2;
1065 d1 = dm + dm;
1066 for (i = 2; i <= m; ++i) {
1067 dm -= 1.;
1068 d1 -= 2.;
1069 d2 -= d1;
1070 ratio = (d3 + d2) / (twox + d1 - ratio);
1071 blpha = (ratio + ratio * blpha) / dm;
1072 }
1073 bk1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * gnm_sqrt(ex));
1074 if (*ize == 1)
1075 bk1 *= gnm_exp(-ex);
1076 wminf = estf[4] * (ex - gnm_abs(ex - estf[6])) + estf[5];
1077 }
1078 /* ---------------------------------------------------------
1079 Calculation of K(ALPHA+1,X)
1080 from K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X)
1081 --------------------------------------------------------- */
1082 bk2 = bk1 + bk1 * (nu + .5 - ratio) / ex;
1083 }
1084 /*--------------------------------------------------------------------
1085 Calculation of 'NCALC', K(ALPHA+I,X), I = 0, 1, ... , NCALC-1,
1086 & K(ALPHA+I,X)/K(ALPHA+I-1,X), I = NCALC, NCALC+1, ... , NB-1
1087 -------------------------------------------------------------------*/
1088 *ncalc = *nb;
1089 bk[0] = bk1;
1090 if (iend == 0)
1091 return;
1092
1093 j = 1 - k;
1094 if (j >= 0)
1095 bk[j] = bk2;
1096
1097 if (iend == 1)
1098 return;
1099
1100 m = imin2((long) (wminf - nu),iend);
1101 for (i = 2; i <= m; ++i) {
1102 t1 = bk1;
1103 bk1 = bk2;
1104 twonu += 2.;
1105 if (ex < 1.) {
1106 if (bk1 >= GNM_MAX / twonu * ex)
1107 break;
1108 } else {
1109 if (bk1 / ex >= GNM_MAX / twonu)
1110 break;
1111 }
1112 bk2 = twonu / ex * bk1 + t1;
1113 ii = i;
1114 ++j;
1115 if (j >= 0) {
1116 bk[j] = bk2;
1117 }
1118 }
1119
1120 m = ii;
1121 if (m == iend) {
1122 return;
1123 }
1124 ratio = bk2 / bk1;
1125 mplus1 = m + 1;
1126 *ncalc = -1;
1127 for (i = mplus1; i <= iend; ++i) {
1128 twonu += 2.;
1129 ratio = twonu / ex + 1./ratio;
1130 ++j;
1131 if (j >= 1) {
1132 bk[j] = ratio;
1133 } else {
1134 if (bk2 >= GNM_MAX / ratio)
1135 return;
1136
1137 bk2 *= ratio;
1138 }
1139 }
1140 *ncalc = imax2(1, mplus1 - k);
1141 if (*ncalc == 1)
1142 bk[0] = bk2;
1143 if (*nb == 1)
1144 return;
1145
1146 L420:
1147 for (i = *ncalc; i < *nb; ++i) { /* i == *ncalc */
1148 #ifndef IEEE_754
1149 if (bk[i-1] >= GNM_MAX / bk[i])
1150 return;
1151 #endif
1152 bk[i] *= bk[i-1];
1153 (*ncalc)++;
1154 }
1155 }
1156 }
1157
1158 /* ------------------------------------------------------------------------ */
1159
1160 static gboolean
1161 bessel_ij_series_domain (gnm_float x, gnm_float v)
1162 {
1163 // The taylor series is valid for all possible values of x and v,
1164 // but it isn't efficient and practical for all.
1165 //
1166 // Since bessel_j is used for computing BesselY when v is not an
1167 // integer, we need the domain to be independent of the sign of v
1168 // for such v.
1169 //
1170 // That is a bit of a problem because when v < -0.5 and close
1171 // to an integer, |v+k| will get close to 0 and the terms will
1172 // suddenly jump in size. The jump will not be larger than a
1173 // factor of 2^53, so as long as [-v]! is much larger than that
1174 // we do not have a problem.
1175
1176 if (v < 0 && v == gnm_floor (v))
1177 return FALSE;
1178
1179 // For non-negative v, the factorials will dominate after the
1180 // k'th term if x*x/4 < c*k(v+k). Let's require two bit decreases
1181 // (c=0.25) from k=10. We ignore the sign of v, see above.
1182
1183 return (x * x / 4 < 0.25 * 10 * (gnm_abs (v) + 10));
1184 }
1185
1186
1187 static GnmQuad
1188 bessel_ij_series (gnm_float x, gnm_float v, gboolean qj)
1189 {
1190 GnmQuad qv, qxh, qfv, qs, qt;
1191 int efv;
1192 void *state = gnm_quad_start ();
1193 gnm_float e, s;
1194
1195 gnm_quad_init (&qxh, x / 2);
1196 gnm_quad_init (&qv, v);
1197
1198 gnm_quad_pow (&qt, &e, &qxh, &qv);
1199
1200 switch (qfactf (v, &qfv, &efv)) {
1201 case 0:
1202 gnm_quad_div (&qt, &qt, &qfv);
1203 e -= efv;
1204 break;
1205 case 1:
1206 qt = gnm_quad_zero;
1207 e = 0;
1208 break;
1209 default:
1210 gnm_quad_init (&qt, gnm_nan);
1211 e = 0;
1212 break;
1213 }
1214
1215 qs = qt;
1216 s = gnm_quad_value (&qs);
1217 if (gnm_finite (s) && s != 0) {
1218 int k, mink = 5;
1219 GnmQuad qxh2;
1220
1221 gnm_quad_mul (&qxh2, &qxh, &qxh);
1222
1223 if (v < 0) {
1224 // Terms can get suddenly big for k ~ -v.
1225 gnm_float ltn0 = -v * (1 - M_LN2gnum + gnm_log (x / -v));
1226 if (ltn0 - (gnm_log (s) + e * M_LN2gnum) < gnm_log (GNM_EPSILON) - 10)
1227 mink = (int)(-v) + 5;
1228 }
1229
1230 for (k = 1; k < 200; k++) {
1231 GnmQuad qa, qb;
1232 gnm_float t;
1233
1234 gnm_quad_mul (&qt, &qt, &qxh2);
1235 gnm_quad_init (&qa, k);
1236 gnm_quad_add (&qb, &qv, &qa);
1237 gnm_quad_init (&qa, qj ? -k : k);
1238 gnm_quad_mul (&qa, &qa, &qb);
1239 gnm_quad_div (&qt, &qt, &qa);
1240 t = gnm_quad_value (&qt);
1241 if (t == 0)
1242 break;
1243 gnm_quad_add (&qs, &qs, &qt);
1244 s = gnm_quad_value (&qs);
1245 if (k >= mink &&
1246 gnm_abs (t) <= GNM_EPSILON / (1 << 20) * gnm_abs (s))
1247 break;
1248 }
1249 }
1250
1251 // We scale here at the end to avoid intermediate overflow
1252 // and underflow.
1253
1254 // Clamp won't affect whether we get 0 or inf.
1255 e = CLAMP (e, G_MININT, G_MAXINT);
1256 qs.h = gnm_ldexp (qs.h, (int)e);
1257 qs.l = gnm_ldexp (qs.l, (int)e);
1258
1259 gnm_quad_end (state);
1260
1261 return qs;
1262 }
1263
1264 /* ------------------------------------------------------------------------ */
1265
1266 static const gnm_float legendre20_roots[(20 + 1) / 2] = {
1267 GNM_const(0.0765265211334973),
1268 GNM_const(0.2277858511416451),
1269 GNM_const(0.3737060887154195),
1270 GNM_const(0.5108670019508271),
1271 GNM_const(0.6360536807265150),
1272 GNM_const(0.7463319064601508),
1273 GNM_const(0.8391169718222188),
1274 GNM_const(0.9122344282513259),
1275 GNM_const(0.9639719272779138),
1276 GNM_const(0.9931285991850949)
1277 };
1278
1279 static const gnm_float legendre20_wts[(20 + 1) / 2] = {
1280 GNM_const(0.1527533871307258),
1281 GNM_const(0.1491729864726037),
1282 GNM_const(0.1420961093183820),
1283 GNM_const(0.1316886384491766),
1284 GNM_const(0.1181945319615184),
1285 GNM_const(0.1019301198172404),
1286 GNM_const(0.0832767415767048),
1287 GNM_const(0.0626720483341091),
1288 GNM_const(0.0406014298003869),
1289 GNM_const(0.0176140071391521)
1290 };
1291
1292 static const gnm_float legendre33_roots[(33 + 1) / 2] = {
1293 GNM_const(0.0000000000000000),
1294 GNM_const(0.0936310658547334),
1295 GNM_const(0.1864392988279916),
1296 GNM_const(0.2776090971524970),
1297 GNM_const(0.3663392577480734),
1298 GNM_const(0.4518500172724507),
1299 GNM_const(0.5333899047863476),
1300 GNM_const(0.6102423458363790),
1301 GNM_const(0.6817319599697428),
1302 GNM_const(0.7472304964495622),
1303 GNM_const(0.8061623562741665),
1304 GNM_const(0.8580096526765041),
1305 GNM_const(0.9023167677434336),
1306 GNM_const(0.9386943726111684),
1307 GNM_const(0.9668229096899927),
1308 GNM_const(0.9864557262306425),
1309 GNM_const(0.9974246942464552)
1310 };
1311
1312 static const gnm_float legendre33_wts[(33 + 1) / 2] = {
1313 GNM_const(0.0937684461602100),
1314 GNM_const(0.0933564260655961),
1315 GNM_const(0.0921239866433168),
1316 GNM_const(0.0900819586606386),
1317 GNM_const(0.0872482876188443),
1318 GNM_const(0.0836478760670387),
1319 GNM_const(0.0793123647948867),
1320 GNM_const(0.0742798548439541),
1321 GNM_const(0.0685945728186567),
1322 GNM_const(0.0623064825303175),
1323 GNM_const(0.0554708466316636),
1324 GNM_const(0.0481477428187117),
1325 GNM_const(0.0404015413316696),
1326 GNM_const(0.0323003586323290),
1327 GNM_const(0.0239155481017495),
1328 GNM_const(0.0153217015129347),
1329 GNM_const(0.0066062278475874)
1330 };
1331
1332 static const gnm_float legendre45_roots[(45 + 1) / 2] = {
1333 GNM_const(0.0000000000000000),
1334 GNM_const(0.0689869801631442),
1335 GNM_const(0.1376452059832530),
1336 GNM_const(0.2056474897832637),
1337 GNM_const(0.2726697697523776),
1338 GNM_const(0.3383926542506022),
1339 GNM_const(0.4025029438585419),
1340 GNM_const(0.4646951239196351),
1341 GNM_const(0.5246728204629161),
1342 GNM_const(0.5821502125693532),
1343 GNM_const(0.6368533944532233),
1344 GNM_const(0.6885216807712006),
1345 GNM_const(0.7369088489454904),
1346 GNM_const(0.7817843125939062),
1347 GNM_const(0.8229342205020863),
1348 GNM_const(0.8601624759606642),
1349 GNM_const(0.8932916717532418),
1350 GNM_const(0.9221639367190004),
1351 GNM_const(0.9466416909956291),
1352 GNM_const(0.9666083103968947),
1353 GNM_const(0.9819687150345405),
1354 GNM_const(0.9926499984472037),
1355 GNM_const(0.9986036451819367)
1356 };
1357
1358 static const gnm_float legendre45_wts[(45 + 1) / 2] = {
1359 GNM_const(0.0690418248292320),
1360 GNM_const(0.0688773169776613),
1361 GNM_const(0.0683845773786697),
1362 GNM_const(0.0675659541636075),
1363 GNM_const(0.0664253484498425),
1364 GNM_const(0.0649681957507234),
1365 GNM_const(0.0632014400738199),
1366 GNM_const(0.0611335008310665),
1367 GNM_const(0.0587742327188417),
1368 GNM_const(0.0561348787597865),
1369 GNM_const(0.0532280167312690),
1370 GNM_const(0.0500674992379520),
1371 GNM_const(0.0466683877183734),
1372 GNM_const(0.0430468807091650),
1373 GNM_const(0.0392202367293025),
1374 GNM_const(0.0352066922016090),
1375 GNM_const(0.0310253749345155),
1376 GNM_const(0.0266962139675777),
1377 GNM_const(0.0222398475505787),
1378 GNM_const(0.0176775352579376),
1379 GNM_const(0.0130311049915828),
1380 GNM_const(0.0083231892962182),
1381 GNM_const(0.0035826631552836)
1382 };
1383
1384 typedef gnm_complex (*ComplexIntegrand) (gnm_float x, const gnm_float *args);
1385
1386 static gnm_complex
1387 complex_legendre_integral (size_t N,
1388 gnm_float L, gnm_float H,
1389 ComplexIntegrand f, const gnm_float *args)
1390 {
1391 const gnm_float *roots;
1392 const gnm_float *wts;
1393 gnm_float m = (L + H) / 2;
1394 gnm_float s = (H - L) / 2;
1395 size_t i;
1396 gnm_complex I = GNM_C0;
1397
1398 switch (N) {
1399 case 20:
1400 roots = legendre20_roots;
1401 wts = legendre20_wts;
1402 break;
1403 case 33:
1404 roots = legendre33_roots;
1405 wts = legendre33_wts;
1406 break;
1407 case 45:
1408 roots = legendre45_roots;
1409 wts = legendre45_wts;
1410 break;
1411 default:
1412 g_assert_not_reached ();
1413 }
1414 if (N & 1)
1415 g_assert (roots[0] == 0.0);
1416
1417 for (i = 0; i < (N + 1) / 2; i++) {
1418 gnm_float r = roots[i];
1419 gnm_float w = wts[i];
1420 int neg;
1421 for (neg = 0; neg <= 1; neg++) {
1422 gnm_float u = neg ? m - s * r : m + s * r;
1423 gnm_complex dI = f (u, args);
1424 I = GNM_CADD (I, GNM_CSCALE (dI, w));
1425 if (i == 0 && (N & 1))
1426 break;
1427 }
1428 }
1429 return GNM_CSCALE (I, s);
1430 }
1431
1432 // Trapezoid rule integraion for a complex function defined on a finite
1433 // interval. This breaks the interval into N uniform pieces.
1434 static void
1435 complex_trapezoid_integral (gnm_complex *res, size_t N,
1436 gnm_float L, gnm_float H,
1437 ComplexIntegrand f, const gnm_float *args)
1438 {
1439 gnm_float s = (H - L) / N;
1440 size_t i;
1441
1442 *res = GNM_C0;
1443 for (i = 0; i <= N; i++) {
1444 gnm_float u = L + i * s;
1445 gnm_complex dI = f (u, args);
1446 if (i == 0 || i == N)
1447 dI = GNM_CSCALE (dI, 0.5);
1448 *res = GNM_CADD (*res, dI);
1449 }
1450 *res = GNM_CSCALE (*res, s);
1451 }
1452
1453 // Shrink integration range to exclude vanishing outer parts.
1454 static void
1455 complex_shink_integral_range (gnm_float *L, gnm_float *H, gnm_float refx,
1456 ComplexIntegrand f,
1457 const gnm_float *args)
1458 {
1459 gnm_complex y;
1460 gnm_float refy, limL = refx, limH = refx;
1461 gboolean first;
1462 gboolean debug = FALSE;
1463
1464 g_return_if_fail (*L <= *H);
1465 g_return_if_fail (*L <= refx && refx <= *H);
1466
1467 y = f (refx, args);
1468 refy = GNM_CABS (y) * GNM_EPSILON;
1469
1470 g_return_if_fail (!gnm_isnan (refy));
1471
1472 if (debug)
1473 g_printerr ("Initial range: (%g,%g) refx=%g refy=%g\n",
1474 *L, *H, refx, refy);
1475
1476 first = TRUE;
1477 while (limL - *L > GNM_EPSILON) {
1478 gnm_float testx = first ? *L : (limL + *L) / 2;
1479 gnm_float testy;
1480
1481 y = f (testx, args);
1482 testy = GNM_CABS (y);
1483
1484 first = FALSE;
1485 if (testy <= refy) {
1486 *L = testx;
1487 if (testy >= refy / 16)
1488 break;
1489 continue;
1490 } else
1491 limL = testx;
1492 }
1493
1494 first = TRUE;
1495 while (*H - limH > GNM_EPSILON) {
1496 gnm_float testx = first ? *H : (*H + limH) / 2;
1497 gnm_float testy;
1498
1499 y = f (testx, args);
1500 testy = GNM_CABS (y);
1501
1502 first = FALSE;
1503 if (testy <= refy) {
1504 *H = testx;
1505 if (testy >= refy / 16)
1506 break;
1507 continue;
1508 } else
1509 limH = testx;
1510 }
1511
1512 if (debug)
1513 g_printerr ("Shrunk range: (%g,%g)\n", *L, *H);
1514 }
1515
1516 typedef gnm_float (*Interpolant) (gnm_float x, const gnm_float *args);
1517
1518 static gnm_float
1519 chebyshev_interpolant (size_t N, gnm_float L, gnm_float H, gnm_float x0,
1520 Interpolant f, const gnm_float *args)
1521 {
1522 size_t i, j, k;
1523 gnm_float *coeffs = g_new (gnm_float, N);
1524 gnm_float m = (L + H) / 2;
1525 gnm_float s = (H - L) / 2;
1526 gnm_float x0n = (x0 - m) / s;
1527 gnm_float res, dip1, dip2, di;
1528
1529 for (j = 0; j < N; j++) {
1530 gnm_float cj = 0;
1531 for (k = 0; k < N; k++) {
1532 gnm_float zj = gnm_cospi ((k + 0.5) / N);
1533 gnm_float xj = m + s * zj;
1534 gnm_float fxj = f (xj, args);
1535 cj += fxj * gnm_cospi (j * (k + 0.5) / N);
1536 }
1537 coeffs[j] = cj * 2.0 / N;
1538 }
1539
1540 dip1 = 0.0;
1541 di = 0.0;
1542
1543 for (i = N - 1; i >= 1; i--) {
1544 dip2 = dip1;
1545 dip1 = di;
1546 di = 2.0 * x0n * dip1 - dip2 + coeffs[i];
1547 }
1548 res = x0n * di - dip1 + 0.5 * coeffs[0];
1549
1550 g_free (coeffs);
1551
1552 if (0) g_printerr ("--> f(%g) = %.20g\n", x0, res);
1553
1554 return res;
1555 }
1556
1557
1558 // Coefficient for the debye polynomial u_n
1559 // Lowest coefficent will be for x^n
1560 // Highest coefficent will be for x^(3n)
1561 // Every second coefficent is zero and left out.
1562 // The count of coefficents will thus be n+1.
1563 static const gnm_float *
1564 debye_u_coeffs (size_t n)
1565 {
1566 static gnm_float **coeffs = NULL;
1567 static size_t nalloc = 0;
1568
1569 if (n >= nalloc) {
1570 size_t i;
1571 coeffs = g_renew (gnm_float *, coeffs, n + 1);
1572 for (i = nalloc; i <= n; i++) {
1573 gnm_float *c = coeffs[i] = g_new0 (gnm_float, i + 1);
1574 gnm_float *l;
1575 size_t j;
1576
1577 if (i == 0) {
1578 c[0] = 1;
1579 continue;
1580 } else if (i == 1) {
1581 c[0] = 1 / 8.0;
1582 c[1] = -5 / 24.0;
1583 continue;
1584 }
1585
1586 l = coeffs[i - 1];
1587
1588 for (j = i; j <= 3 * i; j += 2) {
1589 gnm_float k = 0;
1590
1591 // .5tt u_[i-1]'
1592 if (j < 3 * i)
1593 k += 0.5 * (j-1) * l[((j-1)-(i-1)) / 2];
1594
1595 // -.5tttt u[i-1]'
1596 if (j > i)
1597 k -= 0.5 * (j-3) * l[((j-3)-(i-1)) / 2];
1598
1599 // 1/8*Int[u[i-1]]
1600 if (j < 3 * i)
1601 k += 0.125 * l[((j-1)-(i-1)) / 2] / j;
1602
1603
1604 // -5/8*Int[tt*u[i-1]]
1605 if (j > i)
1606 k -= 0.625 * l[((j-3)-(i-1)) / 2] / j;
1607
1608 c[(j - i) / 2] = k;
1609
1610 if (0)
1611 g_printerr ("Debye u_%d : %g * x^%d\n",
1612 (int)i, k, (int)j);
1613 }
1614 }
1615 nalloc = n + 1;
1616 }
1617
1618 return coeffs[n];
1619 }
1620
1621 static gnm_complex
1622 debye_u (size_t n, gnm_float p, gboolean qip)
1623 {
1624 const gnm_float *coeffs = debye_u_coeffs (n);
1625 gnm_float pn = gnm_pow (p, n);
1626 gnm_float pp = qip ? -p * p : p * p;
1627 gnm_float s = 0;
1628 int i;
1629
1630 for (i = 3 * n; i >= (int)n; i -= 2)
1631 s = s * pp + coeffs[(i - n) / 2];
1632
1633 switch (qip ? n % 4 : 0) {
1634 case 0: return GNM_CREAL (s * pn);
1635 case 1: return GNM_CMAKE (0, s * pn);
1636 case 2: return GNM_CREAL (s * -pn);
1637 case 3: return GNM_CMAKE (0, s * -pn);
1638 default:
1639 g_assert_not_reached ();
1640 }
1641 }
1642
1643
1644 static gnm_float
1645 gnm_sinv_m_v_cosv (gnm_float v, gnm_float sinv)
1646 {
1647 // Deviation: Matviyenko uses direct formula for this for all v.
1648
1649 if (v >= 1)
1650 return sinv - v * gnm_cos (v);
1651 else {
1652 gnm_float r = 0, t = -v;
1653 gnm_float vv = v * v;
1654 int i;
1655
1656 for (i = 3; i < 100; i += 2) {
1657 t = -t * vv / (i * (i == 3 ? 1 : i - 3));
1658 r += t;
1659 if (gnm_abs (t) <= gnm_abs (r) * (GNM_EPSILON / 16))
1660 break;
1661 }
1662
1663 if (0) {
1664 gnm_float ref = sinv - v * gnm_cos (v);
1665 g_printerr ("XXX: %g %g %g -- %g\n",
1666 v, ref, r, r - ref);
1667 }
1668
1669 return r;
1670 }
1671 }
1672
1673 static gnm_float
1674 gnm_sinhumu (gnm_float u)
1675 {
1676 if (!gnm_finite (u))
1677 return u;
1678 else if (gnm_abs (u) >= 1)
1679 return gnm_sinh (u) - u;
1680 else {
1681 gnm_float uu = u * u;
1682 size_t i;
1683 gnm_float s = 0;
1684 gnm_float t = u;
1685
1686 for (i = 3; i < 100; i += 2) {
1687 t *= uu / ((i - 1) * i);
1688 s += t;
1689 if (gnm_abs (t) <= gnm_abs (s) * (GNM_EPSILON / 16))
1690 break;
1691 }
1692 return s;
1693 }
1694 }
1695
1696 /* ------------------------------------------------------------------------ */
1697
1698 static gnm_complex
1699 debye_u_sum (gnm_float x, gnm_float nu,
1700 size_t N, gboolean qalt, gboolean qip)
1701 {
1702 size_t n;
1703 gnm_float f;
1704 gnm_float sqdiff = gnm_abs (x * x - nu * nu);
1705 gnm_float diff2 = gnm_sqrt (sqdiff);
1706 gnm_float p = nu / diff2;
1707 gnm_complex sum = GNM_C0;
1708
1709 (void)debye_u_coeffs (N);
1710
1711 f = 1;
1712 for (n = 0; n <= N; n++) {
1713 gnm_complex t;
1714 if (nu == 0) {
1715 // lim(p/nu,nu->0) = 1/x
1716 gnm_float q = debye_u_coeffs (n)[0] / gnm_pow (x, n);
1717 if (qip && (n & 2)) q = -q;
1718 if (qalt && (n & 1)) q = -q;
1719 if (qip && (n & 1))
1720 t = GNM_CMAKE (0, q);
1721 else
1722 t = GNM_CREAL (q);
1723 } else {
1724 t = debye_u (n, p, qip);
1725 t = GNM_CSCALE (t, f);
1726 f /= nu;
1727 if (qalt) f = -f;
1728 }
1729 sum = GNM_CADD (sum, t);
1730 }
1731
1732 return sum;
1733 }
1734
1735 static void
1736 debye_29_eta1 (gnm_float x, gnm_float nu,
1737 gnm_float *r1a, gnm_float *r1b, gnm_float *rpi)
1738 {
1739 static const gnm_float c[] = {
1740 /* 2 */ 1 / (gnm_float)2,
1741 /* 4 */ 1 / (gnm_float)24,
1742 /* 6 */ 1 / (gnm_float)80,
1743 /* 8 */ 5 / (gnm_float)896,
1744 /* 10 */ 7 / (gnm_float)2304,
1745 /* 12 */ 21 / (gnm_float)11264,
1746 /* 14 */ 33 / (gnm_float)26624,
1747 /* 16 */ 143 / (gnm_float)163840,
1748 /* 18 */ 715 / (gnm_float)1114112,
1749 /* 20 */ 2431 / (gnm_float)4980736
1750 };
1751
1752 gnm_float q = nu / x;
1753
1754 // Deviation: we improve this formula for small nu / x by computing
1755 // eta1 in three parts such that eta1 = (*r1a + *r1b + Pi * *rpi)
1756 // with *r1a small and no rounding errors on the others.
1757
1758 if (q < 0.1) {
1759 gnm_float r = 0;
1760 gnm_float qq = q * q;
1761 unsigned ci;
1762 *rpi = -nu / 2 - 0.25;
1763
1764 for (ci = G_N_ELEMENTS(c); ci-- > 0; )
1765 r = r * qq + c[ci];
1766 *r1a = r * q * nu;
1767 *r1b = x;
1768 } else {
1769 *r1a = gnm_sqrt (x * x - nu * nu) - nu * gnm_acos (nu / x);
1770 *r1b = 0;
1771 *rpi = -0.25;
1772 }
1773 }
1774
1775 static gnm_complex
1776 debye_29 (gnm_float x, gnm_float nu, size_t N)
1777 {
1778 gnm_float sqdiff = x * x - nu * nu;
1779 gnm_float eta1a, eta1b, eta1pi;
1780 gnm_float f1 = gnm_sqrt (2 / M_PIgnum) / gnm_pow (sqdiff, 0.25);
1781 gnm_complex sum, f12;
1782
1783 debye_29_eta1 (x, nu, &eta1a, &eta1b, &eta1pi);
1784
1785 f12 = GNM_CPOLAR (f1, eta1a);
1786 if (eta1b)
1787 f12 = GNM_CMUL (f12, GNM_CPOLAR (1, eta1b));
1788 f12 = GNM_CMUL (f12, GNM_CPOLARPI (1, eta1pi));
1789 sum = debye_u_sum (x, nu, N, TRUE, TRUE);
1790 return GNM_CMUL (f12, sum);
1791 }
1792
1793 static gnm_float
1794 debye_32 (gnm_float x, gnm_float nu, gnm_float eta2, size_t N)
1795 {
1796 gnm_float sqdiff = nu * nu - x * x;
1797 gnm_float f = gnm_exp (-eta2) /
1798 (gnm_sqrt (2 * M_PIgnum) * gnm_pow (sqdiff, 0.25));
1799 gnm_float res;
1800 gnm_complex sum;
1801
1802 sum = debye_u_sum (x, nu, N, FALSE, FALSE);
1803 res = f * sum.re;
1804
1805 if (0)
1806 g_printerr ("D32(%g,%g) = %.20g\n", x, nu, res);
1807
1808 return res;
1809 }
1810
1811 static gnm_float
1812 debye_33 (gnm_float x, gnm_float nu, gnm_float eta2, size_t N)
1813 {
1814 gnm_float sqdiff = nu * nu - x * x;
1815 gnm_float f;
1816 gnm_float res;
1817 gnm_complex sum;
1818 gnm_float c = gnm_sqrt (2 / M_PIgnum);
1819
1820 if (eta2 < gnm_log (GNM_MAX) - 0.01)
1821 f = - c * gnm_exp (eta2) / gnm_pow (sqdiff, 0.25);
1822 else {
1823 // Near-overflow panic
1824 f = - gnm_exp (gnm_log (c) + eta2 - 0.25 * gnm_log (sqdiff));
1825 }
1826
1827 sum = debye_u_sum (x, nu, N, TRUE, FALSE);
1828 res = f * sum.re;
1829
1830 if (0)
1831 g_printerr ("D33(%g,%g) = %.20g\n", x, nu, res);
1832
1833 return res;
1834 }
1835
1836 static gnm_float
1837 integral_83_coshum1 (gnm_float d, gnm_float sinv, gnm_float cosv,
1838 gnm_float sinbeta, gnm_float cosbeta)
1839 {
1840 gnm_float r = 0, todd, teven, dd, cotv;
1841 int i;
1842
1843 if (gnm_abs (d) > 0.1)
1844 return (d * cosbeta - (sinv - sinbeta)) / sinv;
1845
1846 cotv = cosv / sinv;
1847 dd = d * d;
1848 teven = 1;
1849 todd = d;
1850 for (i = 2; i < 100; i++) {
1851 gnm_float t;
1852 if (i & 1) {
1853 todd *= -dd / (i == 3 ? 3 : i * (i - 3));
1854 t = todd * cotv;
1855 } else {
1856 teven *= -dd / (i * (i - 3));
1857 t = teven;
1858 }
1859 r += t;
1860 if (gnm_abs (t) <= gnm_abs (r) * (GNM_EPSILON / 16))
1861 break;
1862 }
1863
1864 if (0) {
1865 gnm_float ref = (d * cosbeta - (sinv - sinbeta)) / sinv;
1866 g_printerr ("coshum1(d=%g): %g %g\n", d, ref, r);
1867 }
1868
1869 return r;
1870 }
1871
1872 static gnm_float
1873 integral_83_cosdiff (gnm_float d, gnm_float v,
1874 gnm_float sinbeta, gnm_float cosbeta)
1875 {
1876 gnm_float s = 0;
1877 gnm_float t = 1;
1878 size_t i;
1879 gboolean debug = FALSE;
1880
1881 g_return_val_if_fail (gnm_abs (d) < 1, gnm_nan);
1882
1883 for (i = 1; i < 100; i += 2) {
1884 t *= -d / i;
1885 s += sinbeta * t;
1886 t *= d / (i + 1);
1887 s += cosbeta * t;
1888 if (gnm_abs (t) <= gnm_abs (s) * (GNM_EPSILON / 16))
1889 break;
1890 }
1891
1892 if (debug) {
1893 gnm_float ref = gnm_cos (v) - cosbeta;
1894 g_printerr ("cosdiff(d=%g): %g %g\n", d, ref, s);
1895 }
1896
1897 return s;
1898 }
1899
1900 static gnm_complex
1901 integral_83_integrand (gnm_float v, gnm_float const *args)
1902 {
1903 gnm_float x = args[0];
1904 gnm_float nu = args[1];
1905 gnm_float beta = args[2];
1906 gnm_float du_dv, phi1, xphi1;
1907 gnm_float sinv = gnm_sin (v);
1908 gboolean debug = FALSE;
1909
1910 if (sinv <= 0) {
1911 // Either end
1912 du_dv = gnm_nan;
1913 phi1 = gnm_ninf;
1914 } else {
1915 gnm_float vmbeta = v - beta;
1916 gnm_float cosv = gnm_cos (v);
1917 gnm_float cosbeta = nu / x;
1918 gnm_float sinbeta = gnm_sqrt (1 - cosbeta * cosbeta);
1919 gnm_float coshum1 = integral_83_coshum1 (vmbeta, sinv, cosv,
1920 sinbeta, cosbeta);
1921 gnm_float sinhu = gnm_sqrt (coshum1 * (coshum1 + 2));
1922 // Deviation: use coshum1 here to avoid cancellation
1923 gnm_float u = gnm_log1p (sinhu + coshum1);
1924 gnm_float du = (gnm_sin (v - beta) -
1925 (v - beta) * cosbeta * cosv);
1926 // Erratum: fix sign of u. See Watson 8.32
1927 if (v < beta)
1928 u = -u, sinhu = -sinhu;
1929 if (gnm_abs (vmbeta) < 0.1) {
1930 phi1 = integral_83_cosdiff (vmbeta, v, sinbeta, cosbeta) * sinhu +
1931 gnm_sinhumu (u) * cosbeta;
1932 } else {
1933 phi1 = cosv * sinhu - cosbeta * u;
1934 }
1935 du_dv = du ? du / (sinhu * sinv * sinv) : 0;
1936 if (debug) {
1937 g_printerr ("beta = %g\n", beta);
1938 g_printerr ("v-beta = %g\n", v-beta);
1939 g_printerr ("phi1 = %g\n", phi1);
1940 g_printerr ("du_dv = %g\n", du_dv);
1941 g_printerr ("coshum1 = %g\n", coshum1);
1942 g_printerr ("sinhu = %g\n", sinhu);
1943 }
1944 }
1945
1946 xphi1 = x * phi1;
1947 if (xphi1 == gnm_ninf) {
1948 // "exp" wins.
1949 return GNM_C0;
1950 } else {
1951 gnm_float exphi1 = gnm_exp (xphi1);
1952 return GNM_CMAKE (du_dv * exphi1, exphi1);
1953 }
1954 }
1955
1956 static gnm_complex
1957 integral_83_alt_integrand (gnm_float t, gnm_float const *args)
1958 {
1959 // v = t^vpow; dv/dt = vpow*t^(vpow-1)
1960 gnm_float vpow = args[3];
1961 return GNM_CSCALE (integral_83_integrand (gnm_pow (t, vpow), args),
1962 vpow * gnm_pow (t, vpow - 1));
1963 }
1964
1965 static gnm_complex
1966 integral_83 (gnm_float x, gnm_float nu, size_t N, gnm_float vpow)
1967 {
1968 // -i/Pi * exp(i*(x*sin(beta)-nu*beta)) *
1969 // Integrate[(du/dv+i)*exp(x*phi1),{v,0,Pi}]
1970 //
1971 // beta = acos(nu/x)
1972 // u = acosh((sin(beta)+(v-beta)*cos(beta))/sin(v))
1973 // du/dv = (sin(v-beta)-(v-beta)*cos(beta)*cos(v)) /
1974 // (sinh(u)*sin^2(v))
1975 // phi1 = cos(v)*sinh(u) - cos(beta)*u
1976 //
1977 // When vpow is not 1 we change variable as v=t^vpow numerically.
1978 // (I.e., the trapezoid points will be evenly spaced in t-space
1979 // instead of v-space.)
1980
1981 // x >= 9 && nu < x - 1.5*crbt(x)
1982
1983 gnm_complex I, f1;
1984 gnm_float beta = gnm_acos (nu / x);
1985 gnm_float xsinbeta = gnm_sqrt (x * x - nu * nu);
1986 gnm_float refx = beta;
1987 gnm_float L = 0;
1988 gnm_float H = M_PIgnum;
1989 gnm_float args[4] = { x, nu, beta, vpow };
1990 ComplexIntegrand integrand;
1991
1992 complex_shink_integral_range (&L, &H, refx,
1993 integral_83_integrand, args);
1994
1995 if (vpow != 1) {
1996 L = gnm_pow (L, 1 / vpow);
1997 H = gnm_pow (H, 1 / vpow);
1998 integrand = integral_83_alt_integrand;
1999 } else {
2000 // We could use the indirect path above, but let's go direct
2001 integrand = integral_83_integrand;
2002 }
2003
2004 complex_trapezoid_integral (&I, N, L, H, integrand, args);
2005
2006 f1 = GNM_CPOLAR (1, xsinbeta - nu * beta);
2007 I = GNM_CMUL (I, f1);
2008 return GNM_CMUL (I, GNM_CMAKE (0, -1 / M_PIgnum));
2009 }
2010
2011 static gnm_complex
2012 integral_105_126_integrand (gnm_float u, gnm_float const *args)
2013 {
2014 gnm_float x = args[0];
2015 gnm_float nu = args[1];
2016 return GNM_CREAL (gnm_exp (x * gnm_sinh (u) - nu * u));
2017 }
2018
2019 static gnm_complex
2020 integral_105_126 (gnm_float x, gnm_float nu, gboolean qH0)
2021 {
2022 // -i/Pi * Integrate[Exp[x*Sinh[u]-nu*u],{u,-Infinity,H}]
2023 // where H is either 0 or alpha, see below.
2024
2025 // Deviation: the analysis doesn't seem to consider the case where
2026 // nu < x which occurs for (126). In that case the integrand takes
2027 // its maximum at 0.
2028
2029 gnm_float args[2] = { x, nu };
2030 gnm_complex I;
2031 gnm_float refx = (nu < x) ? 0 : -gnm_acosh (nu / x);
2032 // For the nu~x case, we have sinh(u)-u = u^3/6 + ...
2033 gnm_float L = refx - MAX (gnm_cbrt (6 * 50 / ((nu + x) / 2)), 50.0 / MIN (nu, x));
2034 gnm_float H = qH0 ? 0 : -refx;
2035
2036 complex_shink_integral_range (&L, &H, refx,
2037 integral_105_126_integrand, args);
2038
2039 I = complex_legendre_integral (45, L, H,
2040 integral_105_126_integrand, args);
2041
2042 return GNM_CMAKE (0, I.re / -M_PIgnum);
2043 }
2044
2045 static gnm_complex
2046 integral_106_integrand (gnm_float v, gnm_float const *args)
2047 {
2048 gnm_float x = args[0];
2049 gnm_float nu = args[1];
2050
2051 gnm_float sinv = gnm_sin (v);
2052 gnm_float coshalpha = nu / x;
2053 gnm_float coshu = coshalpha * (v ? (v / sinv) : 1);
2054 // FIXME: u and sinhu are dubious, numerically
2055 gnm_float u = gnm_acosh (coshu);
2056 gnm_float sinhu = gnm_sinh (u);
2057 gnm_float xphi3 = x * sinhu * gnm_cos (v) - nu * u;
2058 gnm_float exphi3 = gnm_exp (xphi3);
2059
2060 gnm_float num = nu * gnm_sinv_m_v_cosv (v, sinv);
2061 gnm_float den = x * sinv * sinv * sinhu;
2062 gnm_float du_dv = v ? num / den : 0;
2063
2064 return GNM_CMAKE (exphi3 * du_dv, exphi3);
2065 }
2066
2067 static gnm_complex
2068 integral_106 (gnm_float x, gnm_float nu)
2069 {
2070 // -i/Pi * Integrate[Exp[x*phi3[v]]*(i+du/dv),{v,0,Pi}]
2071 //
2072 // alpha = acosh(nu/x)
2073 // u(v) = acosh(nu/x * v/sin(v))
2074 // du/dv = cosh(alpha)*(sin(v)-v*cos(v))/(sin^2(v)*sinh(u(v)))
2075 // phi3(v) = sinh(u)*cos(v) - cosh(alpha)*u(v)
2076
2077 // Note: 2 < x < nu.
2078
2079 gnm_complex I;
2080 gnm_float L = 0, H = M_PIgnum;
2081 gnm_float args[2] = { x, nu };
2082
2083 complex_shink_integral_range (&L, &H, 0, integral_106_integrand, args);
2084
2085 I = complex_legendre_integral (45, L, H,
2086 integral_106_integrand, args);
2087
2088 return GNM_CMUL (I, GNM_CMAKE (0, -1 / M_PIgnum));
2089 }
2090
2091 static gnm_float
2092 integral_127_u (gnm_float v)
2093 {
2094 static const gnm_float c[] = {
2095 GNM_const(0.57735026918962576451),
2096 GNM_const(0.025660011963983367312),
2097 GNM_const(0.0014662863979419067035),
2098 GNM_const(0.000097752426529460446901),
2099 GNM_const(7.4525058224720927532e-6),
2100 GNM_const(6.1544207267743329429e-7),
2101 GNM_const(5.2905118464628039046e-8),
2102 GNM_const(4.6529126736818620163e-9),
2103 GNM_const(4.1606321535886269061e-10),
2104 GNM_const(3.7712142304302013266e-11),
2105 GNM_const(3.4567362099184451359e-12),
2106 GNM_const(3.1977726302920313260e-13),
2107 GNM_const(2.9808441172607163378e-14),
2108 GNM_const(2.7965280211260193677e-15)
2109 };
2110 unsigned ci;
2111 gnm_float vv, u = 0;
2112
2113 if (v >= 1)
2114 return gnm_acosh (v / gnm_sin (v));
2115
2116 // Above formula will suffer from cancellation
2117 vv = v * v;
2118 for (ci = G_N_ELEMENTS(c); ci-- > 0; )
2119 u = u * vv + c[ci];
2120 u *= v;
2121
2122 if (0) {
2123 gnm_float ref = gnm_acosh (v / gnm_sin (v));
2124 g_printerr ("XXX: %g %g\n", ref, u);
2125 }
2126
2127 return u;
2128 }
2129
2130 static gnm_float
2131 integral_127_u_m_sinhu_cos_v (gnm_float v, gnm_float u, gnm_float sinhu)
2132 {
2133 static const gnm_float c[] = {
2134 /* 3 */ GNM_const(0.25660011963983367312),
2135 /* 5 */ GNM_const(0.0),
2136 /* 7 */ GNM_const(0.00097752426529460446901),
2137 /* 9 */ GNM_const(0.000072409204836637368075),
2138 /* 11 */ GNM_const(7.4478039260541292877e-6),
2139 /* 13 */ GNM_const(7.4130822294291683120e-7),
2140 /* 15 */ GNM_const(7.4423844019777464899e-8),
2141 /* 17 */ GNM_const(7.4866591579915856176e-9),
2142 /* 19 */ GNM_const(7.5416412192891756316e-10),
2143 /* 21 */ GNM_const(7.6048685642328096017e-11),
2144 /* 23 */ GNM_const(7.6748139912232122716e-12),
2145 /* 25 */ GNM_const(7.7502621827532506438e-13),
2146 /* 27 */ GNM_const(7.8302824791617646275e-14),
2147 /* 29 */ GNM_const(7.9141968028287716142e-15),
2148 /* 31 */ GNM_const(8.0015150114119176413e-16),
2149 /* 33 */ GNM_const(8.0918754232915038797e-17),
2150 /* 35 */ GNM_const(8.1850043476015809121e-18)
2151 };
2152 unsigned ci;
2153 gnm_float vv, r = 0;
2154
2155 if (v >= 1)
2156 return u - sinhu * gnm_cos (v);
2157
2158 // Above formula will suffer from cancellation
2159 vv = v * v;
2160 for (ci = G_N_ELEMENTS(c); ci-- > 0; )
2161 r = r * vv + c[ci];
2162 r *= v * vv;
2163
2164 if (0) {
2165 gnm_float ref = u - sinhu * gnm_cos (v);
2166 g_printerr ("XXX: %g %g %g\n", ref, r, ref - r);
2167 }
2168
2169 return r;
2170 }
2171
2172 static gnm_complex
2173 integral_127_integrand (gnm_float v, gnm_float const *args)
2174 {
2175 gnm_float x = args[0];
2176 gnm_float nu = args[1];
2177
2178 gnm_float u = integral_127_u (v);
2179 // Deviation: Matviyenko uses taylor expansion for sinh for u < 1.
2180 // There is no need, assuming a reasonable sinh implementation.
2181 gnm_float sinhu = gnm_sinh (u);
2182 gnm_float diff = integral_127_u_m_sinhu_cos_v (v, u, sinhu);
2183 gnm_float sinv = gnm_sin (v);
2184 gnm_float num = gnm_sinv_m_v_cosv (v, sinv);
2185 gnm_float den = sinv * sinv * sinhu;
2186 gnm_float du_dv = v ? num / den : 0;
2187
2188 gnm_complex xphi4, exphi4, i_du_dv;
2189
2190 xphi4 = GNM_CMAKE (x * -diff + (x - nu) * u, (x - nu) * v);
2191 exphi4 = GNM_CEXP (xphi4);
2192 i_du_dv = GNM_CMAKE (du_dv, 1);
2193 return GNM_CMUL (exphi4, i_du_dv);
2194 }
2195
2196 static gnm_complex
2197 integral_127 (gnm_float x, gnm_float nu)
2198 {
2199 // -i/Pi * Integrate[Exp[x*phi4[v]]*(i+du/dv),{v,0,Pi}]
2200 //
2201 // tau = 1-nu/x
2202 // u(v) = acosh(v/sin(v))
2203 // du/dv = (sin(v)-v*cos(v))/(sin^2(v)*sinh(u(v)))
2204 // phi4(v) = sinh(u)*cos(v) - u(v) + tau(u(v) + i * v)
2205
2206 gnm_complex I;
2207 gnm_float L = 0, H = M_PIgnum;
2208 gnm_float args[2] = { x, nu };
2209
2210 complex_shink_integral_range (&L, &H, 0,
2211 integral_127_integrand, args);
2212
2213 I = complex_legendre_integral (33, L, H,
2214 integral_127_integrand, args);
2215
2216 return GNM_CMUL (I, GNM_CMAKE (0, -1 / M_PIgnum));
2217 }
2218
2219 static void
2220 jy_via_j_series (gnm_float x, gnm_float nu, gnm_float *pj, gnm_float *py)
2221 {
2222 void *state = gnm_quad_start ();
2223 GnmQuad qnu, qJnu, qJmnu, qYnu, qCos, qSin, qInvSin;
2224
2225 gnm_quad_init (&qnu, nu);
2226 gnm_quad_cospi (&qCos, &qnu);
2227 gnm_quad_sinpi (&qSin, &qnu);
2228 gnm_quad_div (&qInvSin, &gnm_quad_one, &qSin);
2229
2230 qJnu = bessel_ij_series (x, nu, TRUE);
2231 *pj = gnm_quad_value (&qJnu);
2232
2233 qJmnu = bessel_ij_series (x, -nu, TRUE);
2234
2235 gnm_quad_mul (&qYnu, &qJnu, &qCos);
2236 gnm_quad_sub (&qYnu, &qYnu, &qJmnu);
2237 gnm_quad_mul (&qYnu, &qYnu, &qInvSin);
2238
2239 *py = gnm_quad_value (&qYnu);
2240
2241 gnm_quad_end (state);
2242 }
2243
2244
2245 static gnm_float
2246 cb_y_helper (gnm_float nu, const gnm_float *args)
2247 {
2248 gnm_float x = args[0];
2249 gnm_float Ynu;
2250 if (nu == gnm_floor (nu)) {
2251 g_return_val_if_fail (gnm_abs (nu) < G_MAXINT, gnm_nan);
2252 Ynu = gnm_yn ((int)nu, x);
2253 } else {
2254 gnm_float Jnu;
2255 jy_via_j_series (x, nu, &Jnu, &Ynu);
2256 }
2257 return Ynu;
2258 }
2259
2260 static gnm_complex
2261 hankel1_B1 (gnm_float x, gnm_float nu, size_t N)
2262 {
2263 return debye_29 (x, nu, N);
2264 }
2265
2266 static gnm_complex
2267 hankel1_B2 (gnm_float x, gnm_float nu, size_t N)
2268 {
2269 gnm_float q = nu / x;
2270 gnm_float d = gnm_sqrt (q * q - 1);
2271 gnm_float eta2 = nu * gnm_log (q + d) - gnm_sqrt (nu * nu - x * x);
2272
2273 return GNM_CMAKE (debye_32 (x, nu, eta2, N),
2274 debye_33 (x, nu, eta2, N));
2275 }
2276
2277 static gnm_complex
2278 hankel1_A1 (gnm_float x, gnm_float nu)
2279 {
2280 gnm_float rnu = gnm_floor (nu + 0.49); // Close enough
2281 gnm_float Jnu, Ynu;
2282 gboolean use_yn = (gnm_abs (rnu) < 100000 - 1);
2283
2284 if (gnm_abs (nu - rnu) > 5e-4) {
2285 jy_via_j_series (x, nu, &Jnu, &Ynu);
2286 } else if (use_yn && nu == rnu) {
2287 Jnu = gnm_jn ((int)rnu, x);
2288 Ynu = gnm_yn ((int)rnu, x);
2289 } else {
2290 GnmQuad qJnu = bessel_ij_series (x, nu, TRUE);
2291 size_t N = 6;
2292 gnm_float dnu = 1e-3;
2293 gnm_float args[1] = { x };
2294 gnm_float nul = rnu - dnu, nur = rnu + dnu;
2295 if (use_yn)
2296 N |= 1; // Odd, so we use rnu
2297 Ynu = chebyshev_interpolant (N, nul, nur, nu,
2298 cb_y_helper, args);
2299
2300 Jnu = gnm_quad_value (&qJnu);
2301 }
2302
2303 return GNM_CMAKE (Jnu, Ynu);
2304 }
2305
2306 static gnm_complex
2307 hankel1_A2 (gnm_float x, gnm_float nu)
2308 {
2309 return GNM_CADD (integral_105_126 (x, nu, FALSE),
2310 integral_106 (x, nu));
2311 }
2312
2313 static gnm_complex
2314 hankel1_A3 (gnm_float x, gnm_float nu, gnm_float g)
2315 {
2316 // Deviation: Matviyenko says to change variable to v = t^4 for g < 5.
2317 // That works wonders for BesselJ[9,12.5], but is too sudden for
2318 // BesselJ[1,10]. Instead, gradually move from power of 1 to power
2319 // of 4. The was the power is lowered is ad hoc.
2320 //
2321 // Also, we up the number of points from 37 to 47.
2322
2323 if (g > 5)
2324 return integral_83 (x, nu, 25, 1);
2325 else if (g > 4)
2326 return integral_83 (x, nu, 47, 2);
2327 else if (g > 3)
2328 return integral_83 (x, nu, 47, 3);
2329 else
2330 return integral_83 (x, nu, 47, 4);
2331 }
2332
2333 static gnm_complex
2334 hankel1_A4 (gnm_float x, gnm_float nu)
2335 {
2336 // Deviation: when Matviyenko says that (126) is the same as (105)
2337 // with alpha=0, he is glossing over the finer points. When he should
2338 // have said is that the alpha in the limit is replaced by 0 and
2339 // that the cosh(alpha) inside is replaced textually by nu/x. (We may
2340 // have nu<x and alpha isn't even defined in that case.)
2341 return GNM_CADD (integral_105_126 (x, nu, TRUE),
2342 integral_127 (x, nu));
2343 }
2344
2345 /* ------------------------------------------------------------------------ */
2346
2347 // This follows "On the Evaluation of Bessel Functions" by Gregory Matviyenko.
2348 // Research report YALEU/DCS/RR-903, Yale, Dept. of Comp. Sci., 1992.
2349 //
2350 // Note: there are a few deviations are fixes in this code. These are marked
2351 // with "deviation" or "erratum".
2352
2353 static gnm_complex
2354 hankel1 (gnm_float x, gnm_float nu)
2355 {
2356 gnm_float cbrtx, g;
2357
2358 if (gnm_isnan (x) || gnm_isnan (nu))
2359 return GNM_CNAN;
2360
2361 g_return_val_if_fail (x >= 0, GNM_CNAN);
2362
2363 // Deviation: make this work for negative nu also.
2364 if (nu < 0) {
2365 gnm_complex Hmnu = hankel1 (x, -nu);
2366 if (0) g_printerr ("H_{%g,%g} = %.20g + %.20g * i\n", -nu, x, Hmnu.re, Hmnu.im);
2367 return GNM_CMUL (Hmnu, GNM_CPOLARPI (1, -nu));
2368 }
2369
2370 cbrtx = gnm_cbrt (x);
2371 g = gnm_abs (x - nu) / cbrtx;
2372
2373 if (x >= 17 && g >= 6.5) {
2374 // Algorithm B
2375 size_t N;
2376 if (g < 7)
2377 N = 17;
2378 else if (g < 10)
2379 N = 13;
2380 else if (g < 23)
2381 N = 9;
2382 else
2383 N = 5;
2384
2385 if (nu < x)
2386 return hankel1_B1 (x, nu, N);
2387 else
2388 return hankel1_B2 (x, nu, N);
2389 } else {
2390 // Algorithm A
2391 // Deviation: we use the series on a wider domain as our
2392 // series code uses quad precision.
2393 if (bessel_ij_series_domain (x, nu))
2394 return hankel1_A1 (x, nu);
2395 else if (nu > x && g > 1.5)
2396 return hankel1_A2 (x, nu);
2397 else if (x >= 9 && nu < x && g > 1.5)
2398 return hankel1_A3 (x, nu, g);
2399 else
2400 return hankel1_A4 (x, nu);
2401 }
2402 }
2403
2404 /* ------------------------------------------------------------------------ */
2405
2406 static gboolean
2407 bessel_jy_phase_domain (gnm_float x, gnm_float nu)
2408 {
2409 gnm_float anu = gnm_abs (nu);
2410 gnm_float ax = gnm_abs (x);
2411
2412 if (anu < 2)
2413 return ax > 1000000;
2414
2415 if (ax < 20)
2416 return anu < ax / 5;
2417 if (ax < 30)
2418 return anu < ax / 3;
2419 if (ax < 50)
2420 return anu < ax / 2;
2421 if (ax < 100)
2422 return anu < ax / 1.5;
2423 if (ax < 250)
2424 return anu < ax / 1.2;
2425
2426 return anu < ax / 1.1;
2427 }
2428
2429
2430 static gnm_float
2431 gnm_bessel_M (gnm_float z, gnm_float nu)
2432 {
2433 gnm_float s = 1;
2434 gnm_float tn_z2n = 1;
2435 gnm_float z2 = z * z;
2436 gnm_float nu2 = nu * nu;
2437 int n, NMAX = 400;
2438 gboolean debug = FALSE;
2439
2440 if (debug) g_printerr ("M[%g,%g]\n", nu, z);
2441
2442 // log2(1.1^400) = 55.00
2443
2444 for (n = 1; n < NMAX; n++) {
2445 gnm_float nmh = n - 0.5;
2446 gnm_float f = (nu2 - nmh * nmh) * (nmh / n);
2447 gnm_float r = f / z2;
2448 if (gnm_abs (r) > 1) {
2449 if (debug) g_printerr ("Ratio %g\n", r);
2450 break;
2451 }
2452 tn_z2n *= r;
2453 s += tn_z2n;
2454 if (debug) g_printerr ("Step %d: %g (%g)\n", n, s, tn_z2n);
2455 if (gnm_abs (tn_z2n) < GNM_EPSILON * gnm_abs (s)) {
2456 if (debug) g_printerr ("Precision ok\n");
2457 break;
2458 }
2459 }
2460
2461 return gnm_sqrt (s / (z * (M_PIgnum / 2)));
2462 }
2463
2464
2465 static void
2466 gnm_quad_reduce_pi (GnmQuad *res, GnmQuad const *a, int p, int *pk)
2467 {
2468 gnm_float k;
2469 GnmQuad qk, qa, qb;
2470 unsigned ui;
2471 static const GnmQuad qh = { 0.5, 0 };
2472 static const gnm_float pi_parts[] = {
2473 +0x1.921fb54442d18p+1,
2474 +0x1.1a62633145c04p-53,
2475 +0x1.707344a40938p-104,
2476 +0x1.114cf98e80414p-155,
2477 +0x1.bea63b139b224p-206,
2478 +0x1.14a08798e3404p-258,
2479 +0x1.bbdf2a33679a4p-311,
2480 +0x1.a431b302b0a6cp-362,
2481 +0x1.f25f14374fe1p-414,
2482 +0x1.ab6b6a8e122fp-465
2483 };
2484
2485 if (a->h < 0) {
2486 GnmQuad ma;
2487 ma.h = -a->h;
2488 ma.l = -a->l;
2489 gnm_quad_reduce_pi (res, &ma, p, pk);
2490 res->h = -res->h;
2491 res->l = -res->l;
2492 *pk = (p >= 0) ? (-*pk) & ((1 << (p + 1)) - 1) : 0;
2493 return;
2494 }
2495
2496 if (a->h > gnm_ldexp (1.0, GNM_MANT_DIG))
2497 g_warning ("Reduced accuracy for very large trigonometric arguments");
2498
2499 gnm_quad_div (&qk, a, &gnm_quad_pi);
2500 qk.h = gnm_ldexp (qk.h, p);
2501 qk.l = gnm_ldexp (qk.l, p);
2502
2503 gnm_quad_add (&qk, &qk, &qh);
2504 gnm_quad_floor (&qk, &qk);
2505
2506 k = gnm_quad_value (&qk);
2507 *pk = (int)(gnm_fmod (k, 1 << (1 + p)));
2508
2509 k = gnm_ldexp (k, -p);
2510 qa = *a;
2511 for (ui = 0; ui < G_N_ELEMENTS(pi_parts); ui++) {
2512 gnm_quad_mul12 (&qb, pi_parts[ui], k);
2513 gnm_quad_sub (&qa, &qa, &qb);
2514 }
2515
2516 *res = qa;
2517 }
2518
2519
2520 static gnm_float
2521 gnm_bessel_phi (gnm_float z, gnm_float nu, int *n_pi_4)
2522 {
2523 void *state = gnm_quad_start ();
2524 GnmQuad qs = gnm_quad_zero;
2525 GnmQuad tn_z2n[400], sn_z2n[400];
2526 GnmQuad qz, qnu, qzm2, qnu2, nuh, qrz;
2527 int n, N, NMAX = 400, j, dk;
2528 gnm_float rnu;
2529 gnm_float lt = GNM_MAX;
2530 gboolean debug = FALSE;
2531
2532 if (debug) g_printerr ("Phi[%g,%g]\n", nu, z);
2533
2534 go_quad_init (&qz, z);
2535 go_quad_init (&qnu, nu);
2536
2537 // qzm2 = 1 / (z * z)
2538 go_quad_mul12 (&qzm2, z, z);
2539 go_quad_div (&qzm2, &gnm_quad_one, &qzm2);
2540
2541 // qnu2 = nu * nu
2542 go_quad_mul12 (&qnu2, nu, nu);
2543
2544 (void)gnm_frexp (z / nu, &N);
2545 N = GNM_MANT_DIG / N + 1;
2546 N = MIN (N, (int)G_N_ELEMENTS (tn_z2n));
2547
2548 tn_z2n[0] = sn_z2n[0] = gnm_quad_one;
2549
2550 for (n = 1; n < NMAX; n++) {
2551 GnmQuad qnmh, qnmh2, qf, qd, qn;
2552 gnm_float lt2;
2553
2554 gnm_quad_init (&qn, n);
2555
2556 // qnmh = n - 0.5
2557 gnm_quad_init (&qnmh, n - 0.5);
2558
2559 // qnmh2 = (n - 0.5)^2
2560 gnm_quad_mul (&qnmh2, &qnmh, &qnmh);
2561
2562 // qf = (nu^2 - (n - 0.5)^2) * (n - 0.5) / n
2563 gnm_quad_sub (&qf, &qnu2, &qnmh2);
2564 gnm_quad_mul (&qf, &qf, &qnmh);
2565 gnm_quad_div (&qf, &qf, &qn);
2566
2567 // tn_z2n[n] = tn_z2n[n-1] * f / z^2
2568 gnm_quad_mul (&tn_z2n[n], &tn_z2n[n - 1], &qf);
2569 gnm_quad_mul (&tn_z2n[n], &tn_z2n[n], &qzm2);
2570
2571 sn_z2n[n] = gnm_quad_zero;
2572 for (j = 1; j <= n; j++) {
2573 GnmQuad qp;
2574
2575 gnm_quad_mul (&qp, &tn_z2n[j], &sn_z2n[n - j]);
2576 gnm_quad_sub (&sn_z2n[n], &sn_z2n[n], &qp);
2577 }
2578
2579 gnm_quad_init (&qd, 1 - 2 * n);
2580 gnm_quad_div (&qd, &sn_z2n[n], &qd);
2581
2582 // Break out when the tn ratios go the wrong way. The
2583 // sn ratios can have hickups.
2584 lt2 = gnm_abs (gnm_quad_value (&tn_z2n[n]));
2585 if (lt2 > lt) {
2586 if (debug) g_printerr ("t_n ratio %g\n", lt2 / lt);
2587 break;
2588 }
2589 lt = lt2;
2590
2591 gnm_quad_add (&qs, &qs, &qd);
2592 if (debug) g_printerr ("Step %d: %g (%g)\n", n, gnm_quad_value (&qs), gnm_quad_value (&qd));
2593
2594 if (gnm_abs (gnm_quad_value (&qd)) < GNM_EPSILON * GNM_EPSILON * gnm_abs (gnm_quad_value (&qs))) {
2595 if (debug) g_printerr ("Precision ok\n");
2596 break;
2597 }
2598 }
2599 gnm_quad_mul (&qs, &qz, &qs);
2600
2601 // Add z
2602 gnm_quad_reduce_pi (&qrz, &qz, 2, n_pi_4);
2603 gnm_quad_add (&qs, &qs, &qrz);
2604
2605 // Subtract Pi/4
2606 (*n_pi_4)--;
2607
2608 // Subtract nu/2
2609 rnu = rint (-2 * nu);
2610 *n_pi_4 += (int)fmod (rnu, 8);
2611 gnm_quad_init (&nuh, (-2 * nu - rnu) / 4);
2612 gnm_quad_mul (&nuh, &nuh, &gnm_quad_pi);
2613 gnm_quad_add (&qs, &qs, &nuh);
2614
2615 gnm_quad_reduce_pi (&qs, &qs, 2, &dk);
2616 *n_pi_4 += dk;
2617
2618 *n_pi_4 &= 7;
2619
2620 gnm_quad_end (state);
2621
2622 return gnm_quad_value (&qs);
2623 }
2624
2625 /* ------------------------------------------------------------------------ */
2626
2627 static gnm_float
2628 cos_x_plus_n_pi_4 (gnm_float x, int n_pi_4)
2629 {
2630 static const gnm_float SQH = M_SQRT2gnum / 2;
2631
2632 switch (n_pi_4 & 7) {
2633 case 0: return gnm_cos (x);
2634 case 1: return (gnm_cos (x) - gnm_sin (x)) * SQH;
2635 case 2: return -gnm_sin (x);
2636 case 3: return (gnm_cos (x) + gnm_sin (x)) * -SQH;
2637 case 4: return -gnm_cos (x);
2638 case 5: return (gnm_sin (x) - gnm_cos (x)) * SQH;
2639 case 6: return gnm_sin (x);
2640 case 7: return (gnm_cos (x) + gnm_sin (x)) * SQH;
2641 default: g_assert_not_reached ();
2642 }
2643 }
2644
2645 /* ------------------------------------------------------------------------ */
2646
2647 static gnm_float
2648 gnm_bessel_j_phase (gnm_float x, gnm_float nu)
2649 {
2650 int n_pi_4;
2651 gnm_float M = gnm_bessel_M (x, nu);
2652 gnm_float phi = gnm_bessel_phi (x, nu, &n_pi_4);
2653
2654 if (0) g_printerr ("M=%g phi=%.20g + %d * Pi/4\n", M, phi, n_pi_4);
2655
2656 return M * cos_x_plus_n_pi_4 (phi, n_pi_4);
2657 }
2658
2659 static gnm_float
2660 gnm_bessel_y_phase (gnm_float x, gnm_float nu)
2661 {
2662 int n_pi_4;
2663 gnm_float M = gnm_bessel_M (x, nu);
2664 gnm_float phi = gnm_bessel_phi (x, nu, &n_pi_4);
2665 // Adding 6 means we get sin instead.
2666 return M * cos_x_plus_n_pi_4 (phi, n_pi_4 + 6);
2667 }
2668
2669 /* ------------------------------------------------------------------------ */
2670
2671 gnm_float
2672 gnm_bessel_i (gnm_float x, gnm_float alpha)
2673 {
2674 if (gnm_isnan (x) || gnm_isnan (alpha))
2675 return x + alpha;
2676
2677 if (bessel_ij_series_domain (x, alpha)) {
2678 GnmQuad qi = bessel_ij_series (x, alpha, FALSE);
2679 return gnm_quad_value (&qi);
2680 }
2681
2682 if (x < 0) {
2683 if (alpha != gnm_floor (alpha))
2684 return gnm_nan;
2685
2686 return gnm_fmod (alpha, 2) == 0
2687 ? bessel_i (-x, alpha, 1) /* Even for even alpha */
2688 : 0 - bessel_i (-x, alpha, 1); /* Odd for odd alpha */
2689 } else
2690 return bessel_i (x, alpha, 1);
2691 }
2692
2693 gnm_float
2694 gnm_bessel_j (gnm_float x, gnm_float alpha)
2695 {
2696 if (gnm_isnan (x) || gnm_isnan (alpha))
2697 return x + alpha;
2698
2699 if (x < 0) {
2700 if (alpha != gnm_floor (alpha))
2701 return gnm_nan;
2702
2703 return gnm_fmod (alpha, 2) == 0
2704 ? gnm_bessel_j (-x, alpha) /* Even for even alpha */
2705 : 0 - gnm_bessel_j (-x, alpha); /* Odd for odd alpha */
2706 } else if (bessel_jy_phase_domain (x, alpha)) {
2707 return gnm_bessel_j_phase (x, alpha);
2708 } else {
2709 return GNM_CRE (hankel1 (x, alpha));
2710 }
2711 }
2712
2713 gnm_float
2714 gnm_bessel_k (gnm_float x, gnm_float alpha)
2715 {
2716 return bessel_k (x, alpha, 1);
2717 }
2718
2719 gnm_float
2720 gnm_bessel_y (gnm_float x, gnm_float alpha)
2721 {
2722 if (gnm_isnan (x) || gnm_isnan (alpha))
2723 return x + alpha;
2724
2725 if (x < 0) {
2726 if (alpha != gnm_floor (alpha))
2727 return gnm_nan;
2728
2729 return gnm_fmod (alpha, 2) == 0
2730 ? gnm_bessel_y (-x, alpha) /* Even for even alpha */
2731 : 0 - gnm_bessel_y (-x, alpha); /* Odd for odd alpha */
2732 } else if (bessel_jy_phase_domain (x, alpha)) {
2733 return gnm_bessel_y_phase (x, alpha);
2734 } else {
2735 return GNM_CIM (hankel1 (x, alpha));
2736 }
2737 }
2738
2739 /* ------------------------------------------------------------------------- */
gnumeric mirror