← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cephes/j0.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cephes/j0.c
 
1
/*                                                      j0.c
 
2
 *
 
3
 *      Bessel function of order zero
 
4
 *
 
5
 *
 
6
 *
 
7
 * SYNOPSIS:
 
8
 *
 
9
 * double x, y, j0();
 
10
 *
 
11
 * y = j0( x );
 
12
 *
 
13
 *
 
14
 *
 
15
 * DESCRIPTION:
 
16
 *
 
17
 * Returns Bessel function of order zero of the argument.
 
18
 *
 
19
 * The domain is divided into the intervals [0, 5] and
 
20
 * (5, infinity). In the first interval the following rational
 
21
 * approximation is used:
 
22
 *
 
23
 *
 
24
 *        2         2
 
25
 * (w - r  ) (w - r  ) P (w) / Q (w)
 
26
 *       1         2    3       8
 
27
 *
 
28
 *            2
 
29
 * where w = x  and the two r's are zeros of the function.
 
30
 *
 
31
 * In the second interval, the Hankel asymptotic expansion
 
32
 * is employed with two rational functions of degree 6/6
 
33
 * and 7/7.
 
34
 *
 
35
 *
 
36
 *
 
37
 * ACCURACY:
 
38
 *
 
39
 *                      Absolute error:
 
40
 * arithmetic   domain     # trials      peak         rms
 
41
 *    DEC       0, 30       10000       4.4e-17     6.3e-18
 
42
 *    IEEE      0, 30       60000       4.2e-16     1.1e-16
 
43
 *
 
44
 */
 
45
/*                                                     y0.c
 
46
 *
 
47
 *      Bessel function of the second kind, order zero
 
48
 *
 
49
 *
 
50
 *
 
51
 * SYNOPSIS:
 
52
 *
 
53
 * double x, y, y0();
 
54
 *
 
55
 * y = y0( x );
 
56
 *
 
57
 *
 
58
 *
 
59
 * DESCRIPTION:
 
60
 *
 
61
 * Returns Bessel function of the second kind, of order
 
62
 * zero, of the argument.
 
63
 *
 
64
 * The domain is divided into the intervals [0, 5] and
 
65
 * (5, infinity). In the first interval a rational approximation
 
66
 * R(x) is employed to compute
 
67
 *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
 
68
 * Thus a call to j0() is required.
 
69
 *
 
70
 * In the second interval, the Hankel asymptotic expansion
 
71
 * is employed with two rational functions of degree 6/6
 
72
 * and 7/7.
 
73
 *
 
74
 *
 
75
 *
 
76
 * ACCURACY:
 
77
 *
 
78
 *  Absolute error, when y0(x) < 1; else relative error:
 
79
 *
 
80
 * arithmetic   domain     # trials      peak         rms
 
81
 *    DEC       0, 30        9400       7.0e-17     7.9e-18
 
82
 *    IEEE      0, 30       30000       1.3e-15     1.6e-16
 
83
 *
 
84
 */
 
85
 
 
86
/*
 
87
Cephes Math Library Release 2.1:  January, 1989
 
88
Copyright 1984, 1987, 1989 by Stephen L. Moshier
 
89
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
90
*/
 
91
 
 
92
/* Note: all coefficients satisfy the relative error criterion
 
93
 * except YP, YQ which are designed for absolute error. */
 
94
 
 
95
#include "mconf.h"
 
96
 
 
97
#ifdef UNK
 
98
static double PP[7] = {
 
99
  7.96936729297347051624E-4,
 
100
  8.28352392107440799803E-2,
 
101
  1.23953371646414299388E0,
 
102
  5.44725003058768775090E0,
 
103
  8.74716500199817011941E0,
 
104
  5.30324038235394892183E0,
 
105
  9.99999999999999997821E-1,
 
106
};
 
107
static double PQ[7] = {
 
108
  9.24408810558863637013E-4,
 
109
  8.56288474354474431428E-2,
 
110
  1.25352743901058953537E0,
 
111
  5.47097740330417105182E0,
 
112
  8.76190883237069594232E0,
 
113
  5.30605288235394617618E0,
 
114
  1.00000000000000000218E0,
 
115
};
 
116
#endif
 
117
#ifdef DEC
 
118
static unsigned short PP[28] = {
 
119
0035520,0164604,0140733,0054470,
 
120
0037251,0122605,0115356,0107170,
 
121
0040236,0124412,0071500,0056303,
 
122
0040656,0047737,0045720,0045263,
 
123
0041013,0172143,0045004,0142103,
 
124
0040651,0132045,0026241,0026406,
 
125
0040200,0000000,0000000,0000000,
 
126
};
 
127
static unsigned short PQ[28] = {
 
128
0035562,0052006,0070034,0134666,
 
129
0037257,0057055,0055242,0123424,
 
130
0040240,0071626,0046630,0032371,
 
131
0040657,0011077,0032013,0012731,
 
132
0041014,0030307,0050331,0006414,
 
133
0040651,0145457,0065021,0150304,
 
134
0040200,0000000,0000000,0000000,
 
135
};
 
136
#endif
 
137
#ifdef IBMPC
 
138
static unsigned short PP[28] = {
 
139
0x6b27,0x983b,0x1d30,0x3f4a,
 
140
0xd1cf,0xb35d,0x34b0,0x3fb5,
 
141
0x0b98,0x4e68,0xd521,0x3ff3,
 
142
0x0956,0xe97a,0xc9fb,0x4015,
 
143
0x9888,0x6940,0x7e8c,0x4021,
 
144
0x25a1,0xa594,0x3684,0x4015,
 
145
0x0000,0x0000,0x0000,0x3ff0,
 
146
};
 
147
static unsigned short PQ[28] = {
 
148
0x9737,0xce03,0x4a80,0x3f4e,
 
149
0x54e3,0xab54,0xebc5,0x3fb5,
 
150
0x069f,0xc9b3,0x0e72,0x3ff4,
 
151
0x62bb,0xe681,0xe247,0x4015,
 
152
0x21a1,0xea1b,0x8618,0x4021,
 
153
0x3a19,0xed42,0x3965,0x4015,
 
154
0x0000,0x0000,0x0000,0x3ff0,
 
155
};
 
156
#endif
 
157
#ifdef MIEEE
 
158
static unsigned short PP[28] = {
 
159
0x3f4a,0x1d30,0x983b,0x6b27,
 
160
0x3fb5,0x34b0,0xb35d,0xd1cf,
 
161
0x3ff3,0xd521,0x4e68,0x0b98,
 
162
0x4015,0xc9fb,0xe97a,0x0956,
 
163
0x4021,0x7e8c,0x6940,0x9888,
 
164
0x4015,0x3684,0xa594,0x25a1,
 
165
0x3ff0,0x0000,0x0000,0x0000,
 
166
};
 
167
static unsigned short PQ[28] = {
 
168
0x3f4e,0x4a80,0xce03,0x9737,
 
169
0x3fb5,0xebc5,0xab54,0x54e3,
 
170
0x3ff4,0x0e72,0xc9b3,0x069f,
 
171
0x4015,0xe247,0xe681,0x62bb,
 
172
0x4021,0x8618,0xea1b,0x21a1,
 
173
0x4015,0x3965,0xed42,0x3a19,
 
174
0x3ff0,0x0000,0x0000,0x0000,
 
175
};
 
176
#endif
 
177
 
 
178
#ifdef UNK
 
179
static double QP[8] = {
 
180
-1.13663838898469149931E-2,
 
181
-1.28252718670509318512E0,
 
182
-1.95539544257735972385E1,
 
183
-9.32060152123768231369E1,
 
184
-1.77681167980488050595E2,
 
185
-1.47077505154951170175E2,
 
186
-5.14105326766599330220E1,
 
187
-6.05014350600728481186E0,
 
188
};
 
189
static double QQ[7] = {
 
190
/*  1.00000000000000000000E0,*/
 
191
  6.43178256118178023184E1,
 
192
  8.56430025976980587198E2,
 
193
  3.88240183605401609683E3,
 
194
  7.24046774195652478189E3,
 
195
  5.93072701187316984827E3,
 
196
  2.06209331660327847417E3,
 
197
  2.42005740240291393179E2,
 
198
};
 
199
#endif
 
200
#ifdef DEC
 
201
static unsigned short QP[32] = {
 
202
0136472,0035021,0142451,0141115,
 
203
0140244,0024731,0150620,0105642,
 
204
0141234,0067177,0124161,0060141,
 
205
0141672,0064572,0151557,0043036,
 
206
0142061,0127141,0003127,0043517,
 
207
0142023,0011727,0060271,0144544,
 
208
0141515,0122142,0126620,0143150,
 
209
0140701,0115306,0106715,0007344,
 
210
};
 
211
static unsigned short QQ[28] = {
 
212
/*0040200,0000000,0000000,0000000,*/
 
213
0041600,0121272,0004741,0026544,
 
214
0042526,0015605,0105654,0161771,
 
215
0043162,0123155,0165644,0062645,
 
216
0043342,0041675,0167576,0130756,
 
217
0043271,0052720,0165631,0154214,
 
218
0043000,0160576,0034614,0172024,
 
219
0042162,0000570,0030500,0051235,
 
220
};
 
221
#endif
 
222
#ifdef IBMPC
 
223
static unsigned short QP[32] = {
 
224
0x384a,0x38a5,0x4742,0xbf87,
 
225
0x1174,0x3a32,0x853b,0xbff4,
 
226
0x2c0c,0xf50e,0x8dcf,0xc033,
 
227
0xe8c4,0x5a6d,0x4d2f,0xc057,
 
228
0xe8ea,0x20ca,0x35cc,0xc066,
 
229
0x392d,0xec17,0x627a,0xc062,
 
230
0x18cd,0x55b2,0xb48c,0xc049,
 
231
0xa1dd,0xd1b9,0x3358,0xc018,
 
232
};
 
233
static unsigned short QQ[28] = {
 
234
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
235
0x25ac,0x413c,0x1457,0x4050,
 
236
0x9c7f,0xb175,0xc370,0x408a,
 
237
0x8cb5,0xbd74,0x54cd,0x40ae,
 
238
0xd63e,0xbdef,0x4877,0x40bc,
 
239
0x3b11,0x1d73,0x2aba,0x40b7,
 
240
0x9e82,0xc731,0x1c2f,0x40a0,
 
241
0x0a54,0x0628,0x402f,0x406e,
 
242
};
 
243
#endif
 
244
#ifdef MIEEE
 
245
static unsigned short QP[32] = {
 
246
0xbf87,0x4742,0x38a5,0x384a,
 
247
0xbff4,0x853b,0x3a32,0x1174,
 
248
0xc033,0x8dcf,0xf50e,0x2c0c,
 
249
0xc057,0x4d2f,0x5a6d,0xe8c4,
 
250
0xc066,0x35cc,0x20ca,0xe8ea,
 
251
0xc062,0x627a,0xec17,0x392d,
 
252
0xc049,0xb48c,0x55b2,0x18cd,
 
253
0xc018,0x3358,0xd1b9,0xa1dd,
 
254
};
 
255
static unsigned short QQ[28] = {
 
256
/*0x3ff0,0x0000,0x0000,0x0000,*/
 
257
0x4050,0x1457,0x413c,0x25ac,
 
258
0x408a,0xc370,0xb175,0x9c7f,
 
259
0x40ae,0x54cd,0xbd74,0x8cb5,
 
260
0x40bc,0x4877,0xbdef,0xd63e,
 
261
0x40b7,0x2aba,0x1d73,0x3b11,
 
262
0x40a0,0x1c2f,0xc731,0x9e82,
 
263
0x406e,0x402f,0x0628,0x0a54,
 
264
};
 
265
#endif
 
266
 
 
267
 
 
268
#ifdef UNK
 
269
static double YP[8] = {
 
270
 1.55924367855235737965E4,
 
271
-1.46639295903971606143E7,
 
272
 5.43526477051876500413E9,
 
273
-9.82136065717911466409E11,
 
274
 8.75906394395366999549E13,
 
275
-3.46628303384729719441E15,
 
276
 4.42733268572569800351E16,
 
277
-1.84950800436986690637E16,
 
278
};
 
279
static double YQ[7] = {
 
280
/* 1.00000000000000000000E0,*/
 
281
 1.04128353664259848412E3,
 
282
 6.26107330137134956842E5,
 
283
 2.68919633393814121987E8,
 
284
 8.64002487103935000337E10,
 
285
 2.02979612750105546709E13,
 
286
 3.17157752842975028269E15,
 
287
 2.50596256172653059228E17,
 
288
};
 
289
#endif
 
290
#ifdef DEC
 
291
static unsigned short YP[32] = {
 
292
0043563,0120677,0042264,0046166,
 
293
0146137,0140371,0113444,0042260,
 
294
0050241,0175707,0100502,0063344,
 
295
0152144,0125737,0007265,0164526,
 
296
0053637,0051621,0163035,0060546,
 
297
0155105,0004416,0107306,0060023,
 
298
0056035,0045133,0030132,0000024,
 
299
0155603,0065132,0144061,0131732,
 
300
};
 
301
static unsigned short YQ[28] = {
 
302
/*0040200,0000000,0000000,0000000,*/
 
303
0042602,0024422,0135557,0162663,
 
304
0045030,0155665,0044075,0160135,
 
305
0047200,0035432,0105446,0104005,
 
306
0051240,0167331,0056063,0022743,
 
307
0053223,0127746,0025764,0012160,
 
308
0055064,0044206,0177532,0145545,
 
309
0056536,0111375,0163715,0127201,
 
310
};
 
311
#endif
 
312
#ifdef IBMPC
 
313
static unsigned short YP[32] = {
 
314
0x898f,0xe896,0x7437,0x40ce,
 
315
0x8896,0x32e4,0xf81f,0xc16b,
 
316
0x4cdd,0xf028,0x3f78,0x41f4,
 
317
0xbd2b,0xe1d6,0x957b,0xc26c,
 
318
0xac2d,0x3cc3,0xea72,0x42d3,
 
319
0xcc02,0xd1d8,0xa121,0xc328,
 
320
0x4003,0x660b,0xa94b,0x4363,
 
321
0x367b,0x5906,0x6d4b,0xc350,
 
322
};
 
323
static unsigned short YQ[28] = {
 
324
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
325
0xfcb6,0x576d,0x4522,0x4090,
 
326
0xbc0c,0xa907,0x1b76,0x4123,
 
327
0xd101,0x5164,0x0763,0x41b0,
 
328
0x64bc,0x2b86,0x1ddb,0x4234,
 
329
0x828e,0xc57e,0x75fc,0x42b2,
 
330
0x596d,0xdfeb,0x8910,0x4326,
 
331
0xb5d0,0xbcf9,0xd25f,0x438b,
 
332
};
 
333
#endif
 
334
#ifdef MIEEE
 
335
static unsigned short YP[32] = {
 
336
0x40ce,0x7437,0xe896,0x898f,
 
337
0xc16b,0xf81f,0x32e4,0x8896,
 
338
0x41f4,0x3f78,0xf028,0x4cdd,
 
339
0xc26c,0x957b,0xe1d6,0xbd2b,
 
340
0x42d3,0xea72,0x3cc3,0xac2d,
 
341
0xc328,0xa121,0xd1d8,0xcc02,
 
342
0x4363,0xa94b,0x660b,0x4003,
 
343
0xc350,0x6d4b,0x5906,0x367b,
 
344
};
 
345
static unsigned short YQ[28] = {
 
346
/*0x3ff0,0x0000,0x0000,0x0000,*/
 
347
0x4090,0x4522,0x576d,0xfcb6,
 
348
0x4123,0x1b76,0xa907,0xbc0c,
 
349
0x41b0,0x0763,0x5164,0xd101,
 
350
0x4234,0x1ddb,0x2b86,0x64bc,
 
351
0x42b2,0x75fc,0xc57e,0x828e,
 
352
0x4326,0x8910,0xdfeb,0x596d,
 
353
0x438b,0xd25f,0xbcf9,0xb5d0,
 
354
};
 
355
#endif
 
356
 
 
357
#ifdef UNK
 
358
/*  5.783185962946784521175995758455807035071 */
 
359
static double DR1 = 5.78318596294678452118E0;
 
360
/* 30.47126234366208639907816317502275584842 */
 
361
static double DR2 = 3.04712623436620863991E1;
 
362
#endif
 
363
 
 
364
#ifdef DEC
 
365
static unsigned short R1[] = {0040671,0007734,0001061,0056734};
 
366
#define DR1 *(double *)R1
 
367
static unsigned short R2[] = {0041363,0142445,0030416,0165567};
 
368
#define DR2 *(double *)R2
 
369
#endif
 
370
 
 
371
#ifdef IBMPC
 
372
static unsigned short R1[] = {0x2bbb,0x8046,0x21fb,0x4017};
 
373
#define DR1 *(double *)R1
 
374
static unsigned short R2[] = {0xdd6f,0xa621,0x78a4,0x403e};
 
375
#define DR2 *(double *)R2
 
376
#endif
 
377
 
 
378
#ifdef MIEEE
 
379
static unsigned short R1[] = {0x4017,0x21fb,0x8046,0x2bbb};
 
380
#define DR1 *(double *)R1
 
381
static unsigned short R2[] = {0x403e,0x78a4,0xa621,0xdd6f};
 
382
#define DR2 *(double *)R2
 
383
#endif
 
384
 
 
385
#ifdef UNK
 
386
static double RP[4] = {
 
387
-4.79443220978201773821E9,
 
388
 1.95617491946556577543E12,
 
389
-2.49248344360967716204E14,
 
390
 9.70862251047306323952E15,
 
391
};
 
392
static double RQ[8] = {
 
393
/* 1.00000000000000000000E0,*/
 
394
 4.99563147152651017219E2,
 
395
 1.73785401676374683123E5,
 
396
 4.84409658339962045305E7,
 
397
 1.11855537045356834862E10,
 
398
 2.11277520115489217587E12,
 
399
 3.10518229857422583814E14,
 
400
 3.18121955943204943306E16,
 
401
 1.71086294081043136091E18,
 
402
};
 
403
#endif
 
404
#ifdef DEC
 
405
static unsigned short RP[16] = {
 
406
0150216,0161235,0064344,0014450,
 
407
0052343,0135216,0035624,0144153,
 
408
0154142,0130247,0003310,0003667,
 
409
0055411,0173703,0047772,0176635,
 
410
};
 
411
static unsigned short RQ[32] = {
 
412
/*0040200,0000000,0000000,0000000,*/
 
413
0042371,0144025,0032265,0136137,
 
414
0044451,0133131,0132420,0151466,
 
415
0046470,0144641,0072540,0030636,
 
416
0050446,0126600,0045042,0044243,
 
417
0052365,0172633,0110301,0071063,
 
418
0054215,0032424,0062272,0043513,
 
419
0055742,0005013,0171731,0072335,
 
420
0057275,0170646,0036663,0013134,
 
421
};
 
422
#endif
 
423
#ifdef IBMPC
 
424
static unsigned short RP[16] = {
 
425
0x8325,0xad1c,0xdc53,0xc1f1,
 
426
0x990d,0xc772,0x7751,0x427c,
 
427
0x00f7,0xe0d9,0x5614,0xc2ec,
 
428
0x5fb4,0x69ff,0x3ef8,0x4341,
 
429
};
 
430
static unsigned short RQ[32] = {
 
431
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
432
0xb78c,0xa696,0x3902,0x407f,
 
433
0x1a67,0x36a2,0x36cb,0x4105,
 
434
0x0634,0x2eac,0x1934,0x4187,
 
435
0x4914,0x0944,0xd5b0,0x4204,
 
436
0x2e46,0x7218,0xbeb3,0x427e,
 
437
0x48e9,0x8c97,0xa6a2,0x42f1,
 
438
0x2e9c,0x7e7b,0x4141,0x435c,
 
439
0x62cc,0xc7b6,0xbe34,0x43b7,
 
440
};
 
441
#endif
 
442
#ifdef MIEEE
 
443
static unsigned short RP[16] = {
 
444
0xc1f1,0xdc53,0xad1c,0x8325,
 
445
0x427c,0x7751,0xc772,0x990d,
 
446
0xc2ec,0x5614,0xe0d9,0x00f7,
 
447
0x4341,0x3ef8,0x69ff,0x5fb4,
 
448
};
 
449
static unsigned short RQ[32] = {
 
450
/*0x3ff0,0x0000,0x0000,0x0000,*/
 
451
0x407f,0x3902,0xa696,0xb78c,
 
452
0x4105,0x36cb,0x36a2,0x1a67,
 
453
0x4187,0x1934,0x2eac,0x0634,
 
454
0x4204,0xd5b0,0x0944,0x4914,
 
455
0x427e,0xbeb3,0x7218,0x2e46,
 
456
0x42f1,0xa6a2,0x8c97,0x48e9,
 
457
0x435c,0x4141,0x7e7b,0x2e9c,
 
458
0x43b7,0xbe34,0xc7b6,0x62cc,
 
459
};
 
460
#endif
 
461
 
 
462
#ifndef ANSIPROT
 
463
double j0(), polevl(), p1evl(), log(), sin(), cos(), sqrt();
 
464
#endif
 
465
extern double TWOOPI, SQ2OPI, PIO4;
 
466
 
 
467
double j0(x)
 
468
double x;
 
469
{
 
470
double w, z, p, q, xn;
 
471
 
 
472
if( x < 0 )
 
473
        x = -x;
 
474
 
 
475
if( x <= 5.0 )
 
476
        {
 
477
        z = x * x;
 
478
        if( x < 1.0e-5 )
 
479
                return( 1.0 - z/4.0 );
 
480
 
 
481
        p = (z - DR1) * (z - DR2);
 
482
        p = p * polevl( z, RP, 3)/p1evl( z, RQ, 8 );
 
483
        return( p );
 
484
        }
 
485
 
 
486
w = 5.0/x;
 
487
q = 25.0/(x*x);
 
488
p = polevl( q, PP, 6)/polevl( q, PQ, 6 );
 
489
q = polevl( q, QP, 7)/p1evl( q, QQ, 7 );
 
490
xn = x - PIO4;
 
491
p = p * cos(xn) - w * q * sin(xn);
 
492
return( p * SQ2OPI / sqrt(x) );
 
493
}
 
494
 
 
495
/*                                                      y0() 2  */
 
496
/* Bessel function of second kind, order zero   */
 
497
 
 
498
/* Rational approximation coefficients YP[], YQ[] are used here.
 
499
 * The function computed is  y0(x)  -  2 * log(x) * j0(x) / PI,
 
500
 * whose value at x = 0 is  2 * ( log(0.5) + EUL ) / PI
 
501
 * = 0.073804295108687225.
 
502
 */
 
503
 
 
504
/*
 
505
#define PIO4 .78539816339744830962
 
506
#define SQ2OPI .79788456080286535588
 
507
*/
 
508
extern double MAXNUM;
 
509
 
 
510
double y0(x)
 
511
double x;
 
512
{
 
513
double w, z, p, q, xn;
 
514
 
 
515
if( x <= 5.0 )
 
516
        {
 
517
        if( x <= 0.0 )
 
518
                {
 
519
                mtherr( "y0", DOMAIN );
 
520
                return( -MAXNUM );
 
521
                }
 
522
        z = x * x;
 
523
        w = polevl( z, YP, 7) / p1evl( z, YQ, 7 );
 
524
        w += TWOOPI * log(x) * j0(x);
 
525
        return( w );
 
526
        }
 
527
 
 
528
w = 5.0/x;
 
529
z = 25.0 / (x * x);
 
530
p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
 
531
q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
 
532
xn = x - PIO4;
 
533
p = p * sin(xn) + w * q * cos(xn);
 
534
return( p * SQ2OPI / sqrt(x) );
 
535
}