« back to all changes in this revision
Viewing changes to src/cl/pyopencl-bessel-j.cl
- browse files at revision 6
- compare with another revision
- download diff
- download tarball
- view history from revision 6
- Committer: Package Import Robot
- Author(s): Tomasz Rybak
- Date: 2012-06-21 21:18:03 UTC
- mfrom: (2.1.6 sid)
- Revision ID: package-import@ubuntu.com-20120621211803-424gma7uil7lje4u
* New upstream release.
* Depend on free ocl-icd-libopencl1 instead of non-free libraries.
* Add NEWS entry describing OpenCL library dependencies.
* Change my email.
* Move python-pyopengl from Recommends to Suggests because package
can be used for computations without any graphical environment.
* Depend on free ocl-icd-libopencl1 instead of non-free libraries.
* Add NEWS entry describing OpenCL library dependencies.
* Change my email.
* Move python-pyopengl from Recommends to Suggests because package
can be used for computations without any graphical environment.
- files added:
- .pc/replace-setuptools.patch/Makefile.in
- contrib
- contrib/fortran-to-opencl
- pyopencl/characterize
- pyopencl/compyte
- pyopencl/compyte/ndarray
- src/cl
- files removed:
- debian/docs
- files modified:
- .pc/python-versions.patch/setup.py
added
removed
src/cl/pyopencl-bessel-j.cl
1
// Pieced together from Boost C++ and Cephes by
2
// Andreas Kloeckner (C) 2012
3
//
4
// Pieces from:
5
//
6
// Copyright (c) 2006 Xiaogang Zhang, John Maddock
7
// Use, modification and distribution are subject to the
8
// Boost Software License, Version 1.0. (See
9
// http://www.boost.org/LICENSE_1_0.txt)
10
//
11
// Cephes Math Library Release 2.8: June, 2000
12
// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
13
// What you see here may be used freely, but it comes with no support or
14
// guarantee.
15
16
17
18
19
20
#pragma once
21
22
#include <pyopencl-cephes.cl>
23
#include <pyopencl-airy.cl>
24
25
typedef double T;
26
27
// {{{ evaluate_rational
28
29
T evaluate_rational_backend(__constant const T* num, __constant const T* denom, T z, int count)
30
{
31
T s1, s2;
32
if(z <= 1)
33
{
34
s1 = num[count-1];
35
s2 = denom[count-1];
36
for(int i = (int)count - 2; i >= 0; --i)
37
{
38
s1 *= z;
39
s2 *= z;
40
s1 += num[i];
41
s2 += denom[i];
42
}
43
}
44
else
45
{
46
z = 1 / z;
47
s1 = num[0];
48
s2 = denom[0];
49
for(unsigned i = 1; i < count; ++i)
50
{
51
s1 *= z;
52
s2 *= z;
53
s1 += num[i];
54
s2 += denom[i];
55
}
56
}
57
return s1 / s2;
58
}
59
60
// }}}
61
62
#define evaluate_rational(num, denom, z) \
63
evaluate_rational_backend(num, denom, z, sizeof(num)/sizeof(T))
64
65
// {{{ bessel_j0
66
67
__constant const T bessel_j0_P1[] = {
68
-4.1298668500990866786e+11,
69
2.7282507878605942706e+10,
70
-6.2140700423540120665e+08,
71
6.6302997904833794242e+06,
72
-3.6629814655107086448e+04,
73
1.0344222815443188943e+02,
74
-1.2117036164593528341e-01
75
};
76
__constant const T bessel_j0_Q1[] = {
77
2.3883787996332290397e+12,
78
2.6328198300859648632e+10,
79
1.3985097372263433271e+08,
80
4.5612696224219938200e+05,
81
9.3614022392337710626e+02,
82
1.0,
83
0.0
84
};
85
__constant const T bessel_j0_P2[] = {
86
-1.8319397969392084011e+03,
87
-1.2254078161378989535e+04,
88
-7.2879702464464618998e+03,
89
1.0341910641583726701e+04,
90
1.1725046279757103576e+04,
91
4.4176707025325087628e+03,
92
7.4321196680624245801e+02,
93
4.8591703355916499363e+01
94
};
95
__constant const T bessel_j0_Q2[] = {
96
-3.5783478026152301072e+05,
97
2.4599102262586308984e+05,
98
-8.4055062591169562211e+04,
99
1.8680990008359188352e+04,
100
-2.9458766545509337327e+03,
101
3.3307310774649071172e+02,
102
-2.5258076240801555057e+01,
103
1.0
104
};
105
__constant const T bessel_j0_PC[] = {
106
2.2779090197304684302e+04,
107
4.1345386639580765797e+04,
108
2.1170523380864944322e+04,
109
3.4806486443249270347e+03,
110
1.5376201909008354296e+02,
111
8.8961548424210455236e-01
112
};
113
__constant const T bessel_j0_QC[] = {
114
2.2779090197304684318e+04,
115
4.1370412495510416640e+04,
116
2.1215350561880115730e+04,
117
3.5028735138235608207e+03,
118
1.5711159858080893649e+02,
119
1.0
120
};
121
__constant const T bessel_j0_PS[] = {
122
-8.9226600200800094098e+01,
123
-1.8591953644342993800e+02,
124
-1.1183429920482737611e+02,
125
-2.2300261666214198472e+01,
126
-1.2441026745835638459e+00,
127
-8.8033303048680751817e-03
128
};
129
__constant const T bessel_j0_QS[] = {
130
5.7105024128512061905e+03,
131
1.1951131543434613647e+04,
132
7.2642780169211018836e+03,
133
1.4887231232283756582e+03,
134
9.0593769594993125859e+01,
135
1.0
136
};
137
138
T bessel_j0(T x)
139
{
140
const T x1 = 2.4048255576957727686e+00,
141
x2 = 5.5200781102863106496e+00,
142
x11 = 6.160e+02,
143
x12 = -1.42444230422723137837e-03,
144
x21 = 1.4130e+03,
145
x22 = 5.46860286310649596604e-04;
146
147
T value, factor, r, rc, rs;
148
149
if (x < 0)
150
{
151
x = -x; // even function
152
}
153
if (x == 0)
154
{
155
return 1;
156
}
157
if (x <= 4) // x in (0, 4]
158
{
159
T y = x * x;
160
r = evaluate_rational(bessel_j0_P1, bessel_j0_Q1, y);
161
factor = (x + x1) * ((x - x11/256) - x12);
162
value = factor * r;
163
}
164
else if (x <= 8.0) // x in (4, 8]
165
{
166
T y = 1 - (x * x)/64;
167
r = evaluate_rational(bessel_j0_P2, bessel_j0_Q2, y);
168
factor = (x + x2) * ((x - x21/256) - x22);
169
value = factor * r;
170
}
171
else // x in (8, \infty)
172
{
173
T y = 8 / x;
174
T y2 = y * y;
175
T z = x - 0.25f * M_PI;
176
rc = evaluate_rational(bessel_j0_PC, bessel_j0_QC, y2);
177
rs = evaluate_rational(bessel_j0_PS, bessel_j0_QS, y2);
178
factor = sqrt(2 / (x * M_PI));
179
value = factor * (rc * cos(z) - y * rs * sin(z));
180
}
181
182
return value;
183
}
184
185
// }}}
186
187
// {{{ bessel_j1
188
189
__constant const T bessel_j1_P1[] = {
190
-1.4258509801366645672e+11,
191
6.6781041261492395835e+09,
192
-1.1548696764841276794e+08,
193
9.8062904098958257677e+05,
194
-4.4615792982775076130e+03,
195
1.0650724020080236441e+01,
196
-1.0767857011487300348e-02
197
};
198
__constant const T bessel_j1_Q1[] = {
199
4.1868604460820175290e+12,
200
4.2091902282580133541e+10,
201
2.0228375140097033958e+08,
202
5.9117614494174794095e+05,
203
1.0742272239517380498e+03,
204
1.0,
205
0.0
206
};
207
__constant const T bessel_j1_P2[] = {
208
-1.7527881995806511112e+16,
209
1.6608531731299018674e+15,
210
-3.6658018905416665164e+13,
211
3.5580665670910619166e+11,
212
-1.8113931269860667829e+09,
213
5.0793266148011179143e+06,
214
-7.5023342220781607561e+03,
215
4.6179191852758252278e+00
216
};
217
__constant const T bessel_j1_Q2[] = {
218
1.7253905888447681194e+18,
219
1.7128800897135812012e+16,
220
8.4899346165481429307e+13,
221
2.7622777286244082666e+11,
222
6.4872502899596389593e+08,
223
1.1267125065029138050e+06,
224
1.3886978985861357615e+03,
225
1.0
226
};
227
__constant const T bessel_j1_PC[] = {
228
-4.4357578167941278571e+06,
229
-9.9422465050776411957e+06,
230
-6.6033732483649391093e+06,
231
-1.5235293511811373833e+06,
232
-1.0982405543459346727e+05,
233
-1.6116166443246101165e+03,
234
0.0
235
};
236
__constant const T bessel_j1_QC[] = {
237
-4.4357578167941278568e+06,
238
-9.9341243899345856590e+06,
239
-6.5853394797230870728e+06,
240
-1.5118095066341608816e+06,
241
-1.0726385991103820119e+05,
242
-1.4550094401904961825e+03,
243
1.0
244
};
245
__constant const T bessel_j1_PS[] = {
246
3.3220913409857223519e+04,
247
8.5145160675335701966e+04,
248
6.6178836581270835179e+04,
249
1.8494262873223866797e+04,
250
1.7063754290207680021e+03,
251
3.5265133846636032186e+01,
252
0.0
253
};
254
__constant const T bessel_j1_QS[] = {
255
7.0871281941028743574e+05,
256
1.8194580422439972989e+06,
257
1.4194606696037208929e+06,
258
4.0029443582266975117e+05,
259
3.7890229745772202641e+04,
260
8.6383677696049909675e+02,
261
1.0
262
};
263
264
265
T bessel_j1(T x)
266
{
267
const T x1 = 3.8317059702075123156e+00,
268
x2 = 7.0155866698156187535e+00,
269
x11 = 9.810e+02,
270
x12 = -3.2527979248768438556e-04,
271
x21 = 1.7960e+03,
272
x22 = -3.8330184381246462950e-05;
273
274
T value, factor, r, rc, rs, w;
275
276
w = fabs(x);
277
if (x == 0)
278
{
279
return 0;
280
}
281
if (w <= 4) // w in (0, 4]
282
{
283
T y = x * x;
284
r = evaluate_rational(bessel_j1_P1, bessel_j1_Q1, y);
285
factor = w * (w + x1) * ((w - x11/256) - x12);
286
value = factor * r;
287
}
288
else if (w <= 8) // w in (4, 8]
289
{
290
T y = x * x;
291
r = evaluate_rational(bessel_j1_P2, bessel_j1_Q2, y);
292
factor = w * (w + x2) * ((w - x21/256) - x22);
293
value = factor * r;
294
}
295
else // w in (8, \infty)
296
{
297
T y = 8 / w;
298
T y2 = y * y;
299
T z = w - 0.75f * M_PI;
300
rc = evaluate_rational(bessel_j1_PC, bessel_j1_QC, y2);
301
rs = evaluate_rational(bessel_j1_PS, bessel_j1_QS, y2);
302
factor = sqrt(2 / (w * M_PI));
303
value = factor * (rc * cos(z) - y * rs * sin(z));
304
}
305
306
if (x < 0)
307
{
308
value *= -1; // odd function
309
}
310
return value;
311
}
312
313
// }}}
314
315
// {{{ bessel_recur
316
317
/* Reduce the order by backward recurrence.
318
* AMS55 #9.1.27 and 9.1.73.
319
*/
320
321
#define BESSEL_BIG 1.44115188075855872E+17
322
323
double bessel_recur(double *n, double x, double *newn, int cancel )
324
{
325
double pkm2, pkm1, pk, qkm2, qkm1;
326
/* double pkp1; */
327
double k, ans, qk, xk, yk, r, t, kf;
328
const double big = BESSEL_BIG;
329
int nflag, ctr;
330
331
/* continued fraction for Jn(x)/Jn-1(x) */
332
if( *n < 0.0 )
333
nflag = 1;
334
else
335
nflag = 0;
336
337
fstart:
338
339
#if DEBUG
340
printf( "recur: n = %.6e, newn = %.6e, cfrac = ", *n, *newn );
341
#endif
342
343
pkm2 = 0.0;
344
qkm2 = 1.0;
345
pkm1 = x;
346
qkm1 = *n + *n;
347
xk = -x * x;
348
yk = qkm1;
349
ans = 1.0;
350
ctr = 0;
351
do
352
{
353
yk += 2.0;
354
pk = pkm1 * yk + pkm2 * xk;
355
qk = qkm1 * yk + qkm2 * xk;
356
pkm2 = pkm1;
357
pkm1 = pk;
358
qkm2 = qkm1;
359
qkm1 = qk;
360
if( qk != 0 )
361
r = pk/qk;
362
else
363
r = 0.0;
364
if( r != 0 )
365
{
366
t = fabs( (ans - r)/r );
367
ans = r;
368
}
369
else
370
t = 1.0;
371
372
if( ++ctr > 1000 )
373
{
374
//mtherr( "jv", UNDERFLOW );
375
pk = nan((uint)24);
376
377
goto done;
378
}
379
if( t < DBL_EPSILON )
380
goto done;
381
382
if( fabs(pk) > big )
383
{
384
pkm2 /= big;
385
pkm1 /= big;
386
qkm2 /= big;
387
qkm1 /= big;
388
}
389
}
390
while( t > DBL_EPSILON );
391
392
done:
393
394
#if DEBUG
395
printf( "%.6e\n", ans );
396
#endif
397
398
/* Change n to n-1 if n < 0 and the continued fraction is small
399
*/
400
if( nflag > 0 )
401
{
402
if( fabs(ans) < 0.125 )
403
{
404
nflag = -1;
405
*n = *n - 1.0;
406
goto fstart;
407
}
408
}
409
410
411
kf = *newn;
412
413
/* backward recurrence
414
* 2k
415
* J (x) = --- J (x) - J (x)
416
* k-1 x k k+1
417
*/
418
419
pk = 1.0;
420
pkm1 = 1.0/ans;
421
k = *n - 1.0;
422
r = 2 * k;
423
do
424
{
425
pkm2 = (pkm1 * r - pk * x) / x;
426
/* pkp1 = pk; */
427
pk = pkm1;
428
pkm1 = pkm2;
429
r -= 2.0;
430
/*
431
t = fabs(pkp1) + fabs(pk);
432
if( (k > (kf + 2.5)) && (fabs(pkm1) < 0.25*t) )
433
{
434
k -= 1.0;
435
t = x*x;
436
pkm2 = ( (r*(r+2.0)-t)*pk - r*x*pkp1 )/t;
437
pkp1 = pk;
438
pk = pkm1;
439
pkm1 = pkm2;
440
r -= 2.0;
441
}
442
*/
443
k -= 1.0;
444
}
445
while( k > (kf + 0.5) );
446
447
/* Take the larger of the last two iterates
448
* on the theory that it may have less cancellation error.
449
*/
450
451
if( cancel )
452
{
453
if( (kf >= 0.0) && (fabs(pk) > fabs(pkm1)) )
454
{
455
k += 1.0;
456
pkm2 = pk;
457
}
458
}
459
*newn = k;
460
#if DEBUG
461
printf( "newn %.6e rans %.6e\n", k, pkm2 );
462
#endif
463
return( pkm2 );
464
}
465
466
// }}}
467
468
// {{{ bessel_jvs
469
470
#define BESSEL_MAXGAM 171.624376956302725
471
#define BESSEL_MAXLOG 7.09782712893383996843E2
472
473
/* Ascending power series for Jv(x).
474
* AMS55 #9.1.10.
475
*/
476
477
double bessel_jvs(double n, double x)
478
{
479
double t, u, y, z, k;
480
int ex;
481
int sgngam = 1;
482
483
z = -x * x / 4.0;
484
u = 1.0;
485
y = u;
486
k = 1.0;
487
t = 1.0;
488
489
while( t > DBL_EPSILON )
490
{
491
u *= z / (k * (n+k));
492
y += u;
493
k += 1.0;
494
if( y != 0 )
495
t = fabs( u/y );
496
}
497
#if DEBUG
498
printf( "power series=%.5e ", y );
499
#endif
500
t = frexp( 0.5*x, &ex );
501
ex = ex * n;
502
if( (ex > -1023)
503
&& (ex < 1023)
504
&& (n > 0.0)
505
&& (n < (BESSEL_MAXGAM-1.0)) )
506
{
507
t = pow( 0.5*x, n ) / tgamma( n + 1.0 );
508
#if DEBUG
509
printf( "pow(.5*x, %.4e)/gamma(n+1)=%.5e\n", n, t );
510
#endif
511
y *= t;
512
}
513
else
514
{
515
#if DEBUG
516
z = n * log(0.5*x);
517
k = lgamma( n+1.0 );
518
t = z - k;
519
printf( "log pow=%.5e, lgam(%.4e)=%.5e\n", z, n+1.0, k );
520
#else
521
t = n * log(0.5*x) - lgamma(n + 1.0);
522
#endif
523
if( y < 0 )
524
{
525
sgngam = -sgngam;
526
y = -y;
527
}
528
t += log(y);
529
#if DEBUG
530
printf( "log y=%.5e\n", log(y) );
531
#endif
532
if( t < -BESSEL_MAXLOG )
533
{
534
return( 0.0 );
535
}
536
if( t > BESSEL_MAXLOG )
537
{
538
// mtherr( "Jv", OVERFLOW );
539
return( DBL_MAX);
540
}
541
y = sgngam * exp( t );
542
}
543
return(y);
544
}
545
546
// }}}
547
548
// {{{ bessel_jnt
549
550
__constant const double bessel_jnt_PF2[] = {
551
-9.0000000000000000000e-2,
552
8.5714285714285714286e-2
553
};
554
__constant const double bessel_jnt_PF3[] = {
555
1.3671428571428571429e-1,
556
-5.4920634920634920635e-2,
557
-4.4444444444444444444e-3
558
};
559
__constant const double bessel_jnt_PF4[] = {
560
1.3500000000000000000e-3,
561
-1.6036054421768707483e-1,
562
4.2590187590187590188e-2,
563
2.7330447330447330447e-3
564
};
565
__constant const double bessel_jnt_PG1[] = {
566
-2.4285714285714285714e-1,
567
1.4285714285714285714e-2
568
};
569
__constant const double bessel_jnt_PG2[] = {
570
-9.0000000000000000000e-3,
571
1.9396825396825396825e-1,
572
-1.1746031746031746032e-2
573
};
574
__constant const double bessel_jnt_PG3[] = {
575
1.9607142857142857143e-2,
576
-1.5983694083694083694e-1,
577
6.3838383838383838384e-3
578
};
579
580
double bessel_jnt(double n, double x)
581
{
582
double z, zz, z3;
583
double cbn, n23, cbtwo;
584
double ai, aip, bi, bip; /* Airy functions */
585
double nk, fk, gk, pp, qq;
586
double F[5], G[4];
587
int k;
588
589
cbn = cbrt(n);
590
z = (x - n)/cbn;
591
cbtwo = cbrt( 2.0 );
592
593
/* Airy function */
594
zz = -cbtwo * z;
595
airy( zz, &ai, &aip, &bi, &bip );
596
597
/* polynomials in expansion */
598
zz = z * z;
599
z3 = zz * z;
600
F[0] = 1.0;
601
F[1] = -z/5.0;
602
F[2] = cephes_polevl( z3, bessel_jnt_PF2, 1 ) * zz;
603
F[3] = cephes_polevl( z3, bessel_jnt_PF3, 2 );
604
F[4] = cephes_polevl( z3, bessel_jnt_PF4, 3 ) * z;
605
G[0] = 0.3 * zz;
606
G[1] = cephes_polevl( z3, bessel_jnt_PG1, 1 );
607
G[2] = cephes_polevl( z3, bessel_jnt_PG2, 2 ) * z;
608
G[3] = cephes_polevl( z3, bessel_jnt_PG3, 2 ) * zz;
609
#if DEBUG
610
for( k=0; k<=4; k++ )
611
printf( "F[%d] = %.5E\n", k, F[k] );
612
for( k=0; k<=3; k++ )
613
printf( "G[%d] = %.5E\n", k, G[k] );
614
#endif
615
pp = 0.0;
616
qq = 0.0;
617
nk = 1.0;
618
n23 = cbrt( n * n );
619
620
for( k=0; k<=4; k++ )
621
{
622
fk = F[k]*nk;
623
pp += fk;
624
if( k != 4 )
625
{
626
gk = G[k]*nk;
627
qq += gk;
628
}
629
#if DEBUG
630
printf("fk[%d] %.5E, gk[%d] %.5E\n", k, fk, k, gk );
631
#endif
632
nk /= n23;
633
}
634
635
fk = cbtwo * ai * pp/cbn + cbrt(4.0) * aip * qq/n;
636
return(fk);
637
}
638
639
// }}}
640
641
// {{{ bessel_jnx
642
643
__constant const double bessel_jnx_lambda[] = {
644
1.0,
645
1.041666666666666666666667E-1,
646
8.355034722222222222222222E-2,
647
1.282265745563271604938272E-1,
648
2.918490264641404642489712E-1,
649
8.816272674437576524187671E-1,
650
3.321408281862767544702647E+0,
651
1.499576298686255465867237E+1,
652
7.892301301158651813848139E+1,
653
4.744515388682643231611949E+2,
654
3.207490090890661934704328E+3
655
};
656
__constant const double bessel_jnx_mu[] = {
657
1.0,
658
-1.458333333333333333333333E-1,
659
-9.874131944444444444444444E-2,
660
-1.433120539158950617283951E-1,
661
-3.172272026784135480967078E-1,
662
-9.424291479571202491373028E-1,
663
-3.511203040826354261542798E+0,
664
-1.572726362036804512982712E+1,
665
-8.228143909718594444224656E+1,
666
-4.923553705236705240352022E+2,
667
-3.316218568547972508762102E+3
668
};
669
__constant const double bessel_jnx_P1[] = {
670
-2.083333333333333333333333E-1,
671
1.250000000000000000000000E-1
672
};
673
__constant const double bessel_jnx_P2[] = {
674
3.342013888888888888888889E-1,
675
-4.010416666666666666666667E-1,
676
7.031250000000000000000000E-2
677
};
678
__constant const double bessel_jnx_P3[] = {
679
-1.025812596450617283950617E+0,
680
1.846462673611111111111111E+0,
681
-8.912109375000000000000000E-1,
682
7.324218750000000000000000E-2
683
};
684
__constant const double bessel_jnx_P4[] = {
685
4.669584423426247427983539E+0,
686
-1.120700261622299382716049E+1,
687
8.789123535156250000000000E+0,
688
-2.364086914062500000000000E+0,
689
1.121520996093750000000000E-1
690
};
691
__constant const double bessel_jnx_P5[] = {
692
-2.8212072558200244877E1,
693
8.4636217674600734632E1,
694
-9.1818241543240017361E1,
695
4.2534998745388454861E1,
696
-7.3687943594796316964E0,
697
2.27108001708984375E-1
698
};
699
__constant const double bessel_jnx_P6[] = {
700
2.1257013003921712286E2,
701
-7.6525246814118164230E2,
702
1.0599904525279998779E3,
703
-6.9957962737613254123E2,
704
2.1819051174421159048E2,
705
-2.6491430486951555525E1,
706
5.7250142097473144531E-1
707
};
708
__constant const double bessel_jnx_P7[] = {
709
-1.9194576623184069963E3,
710
8.0617221817373093845E3,
711
-1.3586550006434137439E4,
712
1.1655393336864533248E4,
713
-5.3056469786134031084E3,
714
1.2009029132163524628E3,
715
-1.0809091978839465550E2,
716
1.7277275025844573975E0
717
};
718
719
double bessel_jnx(double n, double x)
720
{
721
double zeta, sqz, zz, zp, np;
722
double cbn, n23, t, z, sz;
723
double pp, qq, z32i, zzi;
724
double ak, bk, akl, bkl;
725
int sign, doa, dob, nflg, k, s, tk, tkp1, m;
726
double u[8];
727
double ai, aip, bi, bip;
728
729
/* Test for x very close to n.
730
* Use expansion for transition region if so.
731
*/
732
cbn = cbrt(n);
733
z = (x - n)/cbn;
734
if( fabs(z) <= 0.7 )
735
return( bessel_jnt(n,x) );
736
737
z = x/n;
738
zz = 1.0 - z*z;
739
if( zz == 0.0 )
740
return(0.0);
741
742
if( zz > 0.0 )
743
{
744
sz = sqrt( zz );
745
t = 1.5 * (log( (1.0+sz)/z ) - sz ); /* zeta ** 3/2 */
746
zeta = cbrt( t * t );
747
nflg = 1;
748
}
749
else
750
{
751
sz = sqrt(-zz);
752
t = 1.5 * (sz - acos(1.0/z));
753
zeta = -cbrt( t * t );
754
nflg = -1;
755
}
756
z32i = fabs(1.0/t);
757
sqz = cbrt(t);
758
759
/* Airy function */
760
n23 = cbrt( n * n );
761
t = n23 * zeta;
762
763
#if DEBUG
764
printf("zeta %.5E, Airy(%.5E)\n", zeta, t );
765
#endif
766
airy( t, &ai, &aip, &bi, &bip );
767
768
/* polynomials in expansion */
769
u[0] = 1.0;
770
zzi = 1.0/zz;
771
u[1] = cephes_polevl( zzi, bessel_jnx_P1, 1 )/sz;
772
u[2] = cephes_polevl( zzi, bessel_jnx_P2, 2 )/zz;
773
u[3] = cephes_polevl( zzi, bessel_jnx_P3, 3 )/(sz*zz);
774
pp = zz*zz;
775
u[4] = cephes_polevl( zzi, bessel_jnx_P4, 4 )/pp;
776
u[5] = cephes_polevl( zzi, bessel_jnx_P5, 5 )/(pp*sz);
777
pp *= zz;
778
u[6] = cephes_polevl( zzi, bessel_jnx_P6, 6 )/pp;
779
u[7] = cephes_polevl( zzi, bessel_jnx_P7, 7 )/(pp*sz);
780
781
#if DEBUG
782
for( k=0; k<=7; k++ )
783
printf( "u[%d] = %.5E\n", k, u[k] );
784
#endif
785
786
pp = 0.0;
787
qq = 0.0;
788
np = 1.0;
789
/* flags to stop when terms get larger */
790
doa = 1;
791
dob = 1;
792
akl = DBL_MAX;
793
bkl = DBL_MAX;
794
795
for( k=0; k<=3; k++ )
796
{
797
tk = 2 * k;
798
tkp1 = tk + 1;
799
zp = 1.0;
800
ak = 0.0;
801
bk = 0.0;
802
for( s=0; s<=tk; s++ )
803
{
804
if( doa )
805
{
806
if( (s & 3) > 1 )
807
sign = nflg;
808
else
809
sign = 1;
810
ak += sign * bessel_jnx_mu[s] * zp * u[tk-s];
811
}
812
813
if( dob )
814
{
815
m = tkp1 - s;
816
if( ((m+1) & 3) > 1 )
817
sign = nflg;
818
else
819
sign = 1;
820
bk += sign * bessel_jnx_lambda[s] * zp * u[m];
821
}
822
zp *= z32i;
823
}
824
825
if( doa )
826
{
827
ak *= np;
828
t = fabs(ak);
829
if( t < akl )
830
{
831
akl = t;
832
pp += ak;
833
}
834
else
835
doa = 0;
836
}
837
838
if( dob )
839
{
840
bk += bessel_jnx_lambda[tkp1] * zp * u[0];
841
bk *= -np/sqz;
842
t = fabs(bk);
843
if( t < bkl )
844
{
845
bkl = t;
846
qq += bk;
847
}
848
else
849
dob = 0;
850
}
851
#if DEBUG
852
printf("a[%d] %.5E, b[%d] %.5E\n", k, ak, k, bk );
853
#endif
854
if( np < DBL_EPSILON )
855
break;
856
np /= n*n;
857
}
858
859
/* normalizing factor ( 4*zeta/(1 - z**2) )**1/4 */
860
t = 4.0 * zeta/zz;
861
t = sqrt( sqrt(t) );
862
863
t *= ai*pp/cbrt(n) + aip*qq/(n23*n);
864
return(t);
865
}
866
867
// }}}
868
869
// {{{ bessel_hankel
870
871
/* Hankel's asymptotic expansion
872
* for large x.
873
* AMS55 #9.2.5.
874
*/
875
876
double bessel_hankel( double n, double x )
877
{
878
double t, u, z, k, sign, conv;
879
double p, q, j, m, pp, qq;
880
int flag;
881
882
m = 4.0*n*n;
883
j = 1.0;
884
z = 8.0 * x;
885
k = 1.0;
886
p = 1.0;
887
u = (m - 1.0)/z;
888
q = u;
889
sign = 1.0;
890
conv = 1.0;
891
flag = 0;
892
t = 1.0;
893
pp = 1.0e38;
894
qq = 1.0e38;
895
896
while( t > DBL_EPSILON )
897
{
898
k += 2.0;
899
j += 1.0;
900
sign = -sign;
901
u *= (m - k * k)/(j * z);
902
p += sign * u;
903
k += 2.0;
904
j += 1.0;
905
u *= (m - k * k)/(j * z);
906
q += sign * u;
907
t = fabs(u/p);
908
if( t < conv )
909
{
910
conv = t;
911
qq = q;
912
pp = p;
913
flag = 1;
914
}
915
/* stop if the terms start getting larger */
916
if( (flag != 0) && (t > conv) )
917
{
918
#if DEBUG
919
printf( "Hankel: convergence to %.4E\n", conv );
920
#endif
921
goto hank1;
922
}
923
}
924
925
hank1:
926
u = x - (0.5*n + 0.25) * M_PI;
927
t = sqrt( 2.0/(M_PI*x) ) * ( pp * cos(u) - qq * sin(u) );
928
#if DEBUG
929
printf( "hank: %.6e\n", t );
930
#endif
931
return( t );
932
}
933
// }}}
934
935
// {{{ bessel_jv
936
937
// SciPy says jn has no advantage over jv, so alias the two.
938
939
#define bessel_jn bessel_jv
940
941
double bessel_jv(double n, double x)
942
{
943
double k, q, t, y, an;
944
int i, sign, nint;
945
946
nint = 0; /* Flag for integer n */
947
sign = 1; /* Flag for sign inversion */
948
an = fabs( n );
949
y = floor( an );
950
if( y == an )
951
{
952
nint = 1;
953
i = an - 16384.0 * floor( an/16384.0 );
954
if( n < 0.0 )
955
{
956
if( i & 1 )
957
sign = -sign;
958
n = an;
959
}
960
if( x < 0.0 )
961
{
962
if( i & 1 )
963
sign = -sign;
964
x = -x;
965
}
966
if( n == 0.0 )
967
return( bessel_j0(x) );
968
if( n == 1.0 )
969
return( sign * bessel_j1(x) );
970
}
971
972
if( (x < 0.0) && (y != an) )
973
{
974
// mtherr( "Jv", DOMAIN );
975
// y = 0.0;
976
y = nan((uint)22);
977
goto done;
978
}
979
980
y = fabs(x);
981
982
if( y < DBL_EPSILON )
983
goto underf;
984
985
k = 3.6 * sqrt(y);
986
t = 3.6 * sqrt(an);
987
if( (y < t) && (an > 21.0) )
988
return( sign * bessel_jvs(n,x) );
989
if( (an < k) && (y > 21.0) )
990
return( sign * bessel_hankel(n,x) );
991
992
if( an < 500.0 )
993
{
994
/* Note: if x is too large, the continued
995
* fraction will fail; but then the
996
* Hankel expansion can be used.
997
*/
998
if( nint != 0 )
999
{
1000
k = 0.0;
1001
q = bessel_recur( &n, x, &k, 1 );
1002
if( k == 0.0 )
1003
{
1004
y = bessel_j0(x)/q;
1005
goto done;
1006
}
1007
if( k == 1.0 )
1008
{
1009
y = bessel_j1(x)/q;
1010
goto done;
1011
}
1012
}
1013
1014
if( an > 2.0 * y )
1015
goto rlarger;
1016
1017
if( (n >= 0.0) && (n < 20.0)
1018
&& (y > 6.0) && (y < 20.0) )
1019
{
1020
/* Recur backwards from a larger value of n
1021
*/
1022
rlarger:
1023
k = n;
1024
1025
y = y + an + 1.0;
1026
if( y < 30.0 )
1027
y = 30.0;
1028
y = n + floor(y-n);
1029
q = bessel_recur( &y, x, &k, 0 );
1030
y = bessel_jvs(y,x) * q;
1031
goto done;
1032
}
1033
1034
if( k <= 30.0 )
1035
{
1036
k = 2.0;
1037
}
1038
else if( k < 90.0 )
1039
{
1040
k = (3*k)/4;
1041
}
1042
if( an > (k + 3.0) )
1043
{
1044
if( n < 0.0 )
1045
k = -k;
1046
q = n - floor(n);
1047
k = floor(k) + q;
1048
if( n > 0.0 )
1049
q = bessel_recur( &n, x, &k, 1 );
1050
else
1051
{
1052
t = k;
1053
k = n;
1054
q = bessel_recur( &t, x, &k, 1 );
1055
k = t;
1056
}
1057
if( q == 0.0 )
1058
{
1059
underf:
1060
y = 0.0;
1061
goto done;
1062
}
1063
}
1064
else
1065
{
1066
k = n;
1067
q = 1.0;
1068
}
1069
1070
/* boundary between convergence of
1071
* power series and Hankel expansion
1072
*/
1073
y = fabs(k);
1074
if( y < 26.0 )
1075
t = (0.0083*y + 0.09)*y + 12.9;
1076
else
1077
t = 0.9 * y;
1078
1079
if( x > t )
1080
y = bessel_hankel(k,x);
1081
else
1082
y = bessel_jvs(k,x);
1083
#if DEBUG
1084
printf( "y = %.16e, recur q = %.16e\n", y, q );
1085
#endif
1086
if( n > 0.0 )
1087
y /= q;
1088
else
1089
y *= q;
1090
}
1091
1092
else
1093
{
1094
/* For large n, use the uniform expansion
1095
* or the transitional expansion.
1096
* But if x is of the order of n**2,
1097
* these may blow up, whereas the
1098
* Hankel expansion will then work.
1099
*/
1100
if( n < 0.0 )
1101
{
1102
//mtherr( "Jv", TLOSS );
1103
//y = 0.0;
1104
y = nan((uint)23);
1105
goto done;
1106
}
1107
t = x/n;
1108
t /= n;
1109
if( t > 0.3 )
1110
y = bessel_hankel(n,x);
1111
else
1112
y = bessel_jnx(n,x);
1113
}
1114
1115
done: return( sign * y);
1116
}
1117
1118
// }}}
1119
1120
// vim: fdm=marker
Loggerhead is a web-based interface for Breezy
Version: 2.0.1