← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/special/cephes/j1.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/special/cephes/j1.c
 
1
/*                                                      j1.c
 
2
 *
 
3
 *      Bessel function of order one
 
4
 *
 
5
 *
 
6
 *
 
7
 * SYNOPSIS:
 
8
 *
 
9
 * double x, y, j1();
 
10
 *
 
11
 * y = j1( x );
 
12
 *
 
13
 *
 
14
 *
 
15
 * DESCRIPTION:
 
16
 *
 
17
 * Returns Bessel function of order one of the argument.
 
18
 *
 
19
 * The domain is divided into the intervals [0, 8] and
 
20
 * (8, infinity). In the first interval a 24 term Chebyshev
 
21
 * expansion is used. In the second, the asymptotic
 
22
 * trigonometric representation is employed using two
 
23
 * rational functions of degree 5/5.
 
24
 *
 
25
 *
 
26
 *
 
27
 * ACCURACY:
 
28
 *
 
29
 *                      Absolute error:
 
30
 * arithmetic   domain      # trials      peak         rms
 
31
 *    DEC       0, 30       10000       4.0e-17     1.1e-17
 
32
 *    IEEE      0, 30       30000       2.6e-16     1.1e-16
 
33
 *
 
34
 *
 
35
 */
 
36
/*                                                     y1.c
 
37
 *
 
38
 *      Bessel function of second kind of order one
 
39
 *
 
40
 *
 
41
 *
 
42
 * SYNOPSIS:
 
43
 *
 
44
 * double x, y, y1();
 
45
 *
 
46
 * y = y1( x );
 
47
 *
 
48
 *
 
49
 *
 
50
 * DESCRIPTION:
 
51
 *
 
52
 * Returns Bessel function of the second kind of order one
 
53
 * of the argument.
 
54
 *
 
55
 * The domain is divided into the intervals [0, 8] and
 
56
 * (8, infinity). In the first interval a 25 term Chebyshev
 
57
 * expansion is used, and a call to j1() is required.
 
58
 * In the second, the asymptotic trigonometric representation
 
59
 * is employed using two rational functions of degree 5/5.
 
60
 *
 
61
 *
 
62
 *
 
63
 * ACCURACY:
 
64
 *
 
65
 *                      Absolute error:
 
66
 * arithmetic   domain      # trials      peak         rms
 
67
 *    DEC       0, 30       10000       8.6e-17     1.3e-17
 
68
 *    IEEE      0, 30       30000       1.0e-15     1.3e-16
 
69
 *
 
70
 * (error criterion relative when |y1| > 1).
 
71
 *
 
72
 */
 
73
 
 
74
 
 
75
/*
 
76
Cephes Math Library Release 2.8:  June, 2000
 
77
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
 
78
*/
 
79
 
 
80
/*
 
81
#define PIO4 .78539816339744830962
 
82
#define THPIO4 2.35619449019234492885
 
83
#define SQ2OPI .79788456080286535588
 
84
*/
 
85
 
 
86
#include "mconf.h"
 
87
 
 
88
#ifdef UNK
 
89
static double RP[4] = {
 
90
-8.99971225705559398224E8,
 
91
 4.52228297998194034323E11,
 
92
-7.27494245221818276015E13,
 
93
 3.68295732863852883286E15,
 
94
};
 
95
static double RQ[8] = {
 
96
/* 1.00000000000000000000E0,*/
 
97
 6.20836478118054335476E2,
 
98
 2.56987256757748830383E5,
 
99
 8.35146791431949253037E7,
 
100
 2.21511595479792499675E10,
 
101
 4.74914122079991414898E12,
 
102
 7.84369607876235854894E14,
 
103
 8.95222336184627338078E16,
 
104
 5.32278620332680085395E18,
 
105
};
 
106
#endif
 
107
#ifdef DEC
 
108
static unsigned short RP[16] = {
 
109
0147526,0110742,0063322,0077052,
 
110
0051722,0112720,0065034,0061530,
 
111
0153604,0052227,0033147,0105650,
 
112
0055121,0055025,0032276,0022015,
 
113
};
 
114
static unsigned short RQ[32] = {
 
115
/*0040200,0000000,0000000,0000000,*/
 
116
0042433,0032610,0155604,0033473,
 
117
0044572,0173320,0067270,0006616,
 
118
0046637,0045246,0162225,0006606,
 
119
0050645,0004773,0157577,0053004,
 
120
0052612,0033734,0001667,0176501,
 
121
0054462,0054121,0173147,0121367,
 
122
0056237,0002777,0121451,0176007,
 
123
0057623,0136253,0131601,0044710,
 
124
};
 
125
#endif
 
126
#ifdef IBMPC
 
127
static unsigned short RP[16] = {
 
128
0x4fc5,0x4cda,0xd23c,0xc1ca,
 
129
0x8c6b,0x0d43,0x52ba,0x425a,
 
130
0xf175,0xe6cc,0x8a92,0xc2d0,
 
131
0xc482,0xa697,0x2b42,0x432a,
 
132
};
 
133
static unsigned short RQ[32] = {
 
134
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
135
0x86e7,0x1b70,0x66b1,0x4083,
 
136
0x01b2,0x0dd7,0x5eda,0x410f,
 
137
0xa1b1,0xdc92,0xe954,0x4193,
 
138
0xeac1,0x7bef,0xa13f,0x4214,
 
139
0xffa8,0x8076,0x46fb,0x4291,
 
140
0xf45f,0x3ecc,0x4b0a,0x4306,
 
141
0x3f81,0xf465,0xe0bf,0x4373,
 
142
0x2939,0x7670,0x7795,0x43d2,
 
143
};
 
144
#endif
 
145
#ifdef MIEEE
 
146
static unsigned short RP[16] = {
 
147
0xc1ca,0xd23c,0x4cda,0x4fc5,
 
148
0x425a,0x52ba,0x0d43,0x8c6b,
 
149
0xc2d0,0x8a92,0xe6cc,0xf175,
 
150
0x432a,0x2b42,0xa697,0xc482,
 
151
};
 
152
static unsigned short RQ[32] = {
 
153
/*0x3ff0,0x0000,0x0000,0x0000,*/
 
154
0x4083,0x66b1,0x1b70,0x86e7,
 
155
0x410f,0x5eda,0x0dd7,0x01b2,
 
156
0x4193,0xe954,0xdc92,0xa1b1,
 
157
0x4214,0xa13f,0x7bef,0xeac1,
 
158
0x4291,0x46fb,0x8076,0xffa8,
 
159
0x4306,0x4b0a,0x3ecc,0xf45f,
 
160
0x4373,0xe0bf,0xf465,0x3f81,
 
161
0x43d2,0x7795,0x7670,0x2939,
 
162
};
 
163
#endif
 
164
 
 
165
#ifdef UNK
 
166
static double PP[7] = {
 
167
 7.62125616208173112003E-4,
 
168
 7.31397056940917570436E-2,
 
169
 1.12719608129684925192E0,
 
170
 5.11207951146807644818E0,
 
171
 8.42404590141772420927E0,
 
172
 5.21451598682361504063E0,
 
173
 1.00000000000000000254E0,
 
174
};
 
175
static double PQ[7] = {
 
176
 5.71323128072548699714E-4,
 
177
 6.88455908754495404082E-2,
 
178
 1.10514232634061696926E0,
 
179
 5.07386386128601488557E0,
 
180
 8.39985554327604159757E0,
 
181
 5.20982848682361821619E0,
 
182
 9.99999999999999997461E-1,
 
183
};
 
184
#endif
 
185
#ifdef DEC
 
186
static unsigned short PP[28] = {
 
187
0035507,0144542,0061543,0024326,
 
188
0037225,0145105,0017766,0022661,
 
189
0040220,0043766,0010254,0133255,
 
190
0040643,0113047,0142611,0151521,
 
191
0041006,0144344,0055351,0074261,
 
192
0040646,0156520,0120574,0006416,
 
193
0040200,0000000,0000000,0000000,
 
194
};
 
195
static unsigned short PQ[28] = {
 
196
0035425,0142330,0115041,0165514,
 
197
0037214,0177352,0145105,0052026,
 
198
0040215,0072515,0141207,0073255,
 
199
0040642,0056427,0137222,0106405,
 
200
0041006,0062716,0166427,0165450,
 
201
0040646,0133352,0035425,0123304,
 
202
0040200,0000000,0000000,0000000,
 
203
};
 
204
#endif
 
205
#ifdef IBMPC
 
206
static unsigned short PP[28] = {
 
207
0x651b,0x4c6c,0xf92c,0x3f48,
 
208
0xc4b6,0xa3fe,0xb948,0x3fb2,
 
209
0x96d6,0xc215,0x08fe,0x3ff2,
 
210
0x3a6a,0xf8b1,0x72c4,0x4014,
 
211
0x2f16,0x8b5d,0xd91c,0x4020,
 
212
0x81a2,0x142f,0xdbaa,0x4014,
 
213
0x0000,0x0000,0x0000,0x3ff0,
 
214
};
 
215
static unsigned short PQ[28] = {
 
216
0x3d69,0x1344,0xb89b,0x3f42,
 
217
0xaa83,0x5948,0x9fdd,0x3fb1,
 
218
0xeed6,0xb850,0xaea9,0x3ff1,
 
219
0x51a1,0xf7d2,0x4ba2,0x4014,
 
220
0xfd65,0xdda2,0xccb9,0x4020,
 
221
0xb4d9,0x4762,0xd6dd,0x4014,
 
222
0x0000,0x0000,0x0000,0x3ff0,
 
223
};
 
224
#endif
 
225
#ifdef MIEEE
 
226
static unsigned short PP[28] = {
 
227
0x3f48,0xf92c,0x4c6c,0x651b,
 
228
0x3fb2,0xb948,0xa3fe,0xc4b6,
 
229
0x3ff2,0x08fe,0xc215,0x96d6,
 
230
0x4014,0x72c4,0xf8b1,0x3a6a,
 
231
0x4020,0xd91c,0x8b5d,0x2f16,
 
232
0x4014,0xdbaa,0x142f,0x81a2,
 
233
0x3ff0,0x0000,0x0000,0x0000,
 
234
};
 
235
static unsigned short PQ[28] = {
 
236
0x3f42,0xb89b,0x1344,0x3d69,
 
237
0x3fb1,0x9fdd,0x5948,0xaa83,
 
238
0x3ff1,0xaea9,0xb850,0xeed6,
 
239
0x4014,0x4ba2,0xf7d2,0x51a1,
 
240
0x4020,0xccb9,0xdda2,0xfd65,
 
241
0x4014,0xd6dd,0x4762,0xb4d9,
 
242
0x3ff0,0x0000,0x0000,0x0000,
 
243
};
 
244
#endif
 
245
 
 
246
#ifdef UNK
 
247
static double QP[8] = {
 
248
 5.10862594750176621635E-2,
 
249
 4.98213872951233449420E0,
 
250
 7.58238284132545283818E1,
 
251
 3.66779609360150777800E2,
 
252
 7.10856304998926107277E2,
 
253
 5.97489612400613639965E2,
 
254
 2.11688757100572135698E2,
 
255
 2.52070205858023719784E1,
 
256
};
 
257
static double QQ[7] = {
 
258
/* 1.00000000000000000000E0,*/
 
259
 7.42373277035675149943E1,
 
260
 1.05644886038262816351E3,
 
261
 4.98641058337653607651E3,
 
262
 9.56231892404756170795E3,
 
263
 7.99704160447350683650E3,
 
264
 2.82619278517639096600E3,
 
265
 3.36093607810698293419E2,
 
266
};
 
267
#endif
 
268
#ifdef DEC
 
269
static unsigned short QP[32] = {
 
270
0037121,0037723,0055605,0151004,
 
271
0040637,0066656,0031554,0077264,
 
272
0041627,0122714,0153170,0161466,
 
273
0042267,0061712,0036520,0140145,
 
274
0042461,0133315,0131573,0071176,
 
275
0042425,0057525,0147500,0013201,
 
276
0042123,0130122,0061245,0154131,
 
277
0041311,0123772,0064254,0172650,
 
278
};
 
279
static unsigned short QQ[28] = {
 
280
/*0040200,0000000,0000000,0000000,*/
 
281
0041624,0074603,0002112,0101670,
 
282
0042604,0007135,0010162,0175565,
 
283
0043233,0151510,0157757,0172010,
 
284
0043425,0064506,0112006,0104276,
 
285
0043371,0164125,0032271,0164242,
 
286
0043060,0121425,0122750,0136013,
 
287
0042250,0005773,0053472,0146267,
 
288
};
 
289
#endif
 
290
#ifdef IBMPC
 
291
static unsigned short QP[32] = {
 
292
0xba40,0x6b70,0x27fa,0x3faa,
 
293
0x8fd6,0xc66d,0xedb5,0x4013,
 
294
0x1c67,0x9acf,0xf4b9,0x4052,
 
295
0x180d,0x47aa,0xec79,0x4076,
 
296
0x6e50,0xb66f,0x36d9,0x4086,
 
297
0x02d0,0xb9e8,0xabea,0x4082,
 
298
0xbb0b,0x4c54,0x760a,0x406a,
 
299
0x9eb5,0x4d15,0x34ff,0x4039,
 
300
};
 
301
static unsigned short QQ[28] = {
 
302
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
303
0x5077,0x6089,0x8f30,0x4052,
 
304
0x5f6f,0xa20e,0x81cb,0x4090,
 
305
0xfe81,0x1bfd,0x7a69,0x40b3,
 
306
0xd118,0xd280,0xad28,0x40c2,
 
307
0x3d14,0xa697,0x3d0a,0x40bf,
 
308
0x1781,0xb4bd,0x1462,0x40a6,
 
309
0x5997,0x6ae7,0x017f,0x4075,
 
310
};
 
311
#endif
 
312
#ifdef MIEEE
 
313
static unsigned short QP[32] = {
 
314
0x3faa,0x27fa,0x6b70,0xba40,
 
315
0x4013,0xedb5,0xc66d,0x8fd6,
 
316
0x4052,0xf4b9,0x9acf,0x1c67,
 
317
0x4076,0xec79,0x47aa,0x180d,
 
318
0x4086,0x36d9,0xb66f,0x6e50,
 
319
0x4082,0xabea,0xb9e8,0x02d0,
 
320
0x406a,0x760a,0x4c54,0xbb0b,
 
321
0x4039,0x34ff,0x4d15,0x9eb5,
 
322
};
 
323
static unsigned short QQ[28] = {
 
324
/*0x3ff0,0x0000,0x0000,0x0000,*/
 
325
0x4052,0x8f30,0x6089,0x5077,
 
326
0x4090,0x81cb,0xa20e,0x5f6f,
 
327
0x40b3,0x7a69,0x1bfd,0xfe81,
 
328
0x40c2,0xad28,0xd280,0xd118,
 
329
0x40bf,0x3d0a,0xa697,0x3d14,
 
330
0x40a6,0x1462,0xb4bd,0x1781,
 
331
0x4075,0x017f,0x6ae7,0x5997,
 
332
};
 
333
#endif
 
334
 
 
335
#ifdef UNK
 
336
static double YP[6] = {
 
337
 1.26320474790178026440E9,
 
338
-6.47355876379160291031E11,
 
339
 1.14509511541823727583E14,
 
340
-8.12770255501325109621E15,
 
341
 2.02439475713594898196E17,
 
342
-7.78877196265950026825E17,
 
343
};
 
344
static double YQ[8] = {
 
345
/* 1.00000000000000000000E0,*/
 
346
 5.94301592346128195359E2,
 
347
 2.35564092943068577943E5,
 
348
 7.34811944459721705660E7,
 
349
 1.87601316108706159478E10,
 
350
 3.88231277496238566008E12,
 
351
 6.20557727146953693363E14,
 
352
 6.87141087355300489866E16,
 
353
 3.97270608116560655612E18,
 
354
};
 
355
#endif
 
356
#ifdef DEC
 
357
static unsigned short YP[24] = {
 
358
0047626,0112763,0013715,0133045,
 
359
0152026,0134552,0142033,0024411,
 
360
0053720,0045245,0102210,0077565,
 
361
0155347,0000321,0136415,0102031,
 
362
0056463,0146550,0055633,0032605,
 
363
0157054,0171012,0167361,0054265,
 
364
};
 
365
static unsigned short YQ[32] = {
 
366
/*0040200,0000000,0000000,0000000,*/
 
367
0042424,0111515,0044773,0153014,
 
368
0044546,0005405,0171307,0075774,
 
369
0046614,0023575,0047105,0063556,
 
370
0050613,0143034,0101533,0156026,
 
371
0052541,0175367,0166514,0114257,
 
372
0054415,0014466,0134350,0171154,
 
373
0056164,0017436,0025075,0022101,
 
374
0057534,0103614,0103663,0121772,
 
375
};
 
376
#endif
 
377
#ifdef IBMPC
 
378
static unsigned short YP[24] = {
 
379
0xb6c5,0x62f9,0xd2be,0x41d2,
 
380
0x6521,0x5883,0xd72d,0xc262,
 
381
0x0fef,0xb091,0x0954,0x42da,
 
382
0xb083,0x37a1,0xe01a,0xc33c,
 
383
0x66b1,0x0b73,0x79ad,0x4386,
 
384
0x2b17,0x5dde,0x9e41,0xc3a5,
 
385
};
 
386
static unsigned short YQ[32] = {
 
387
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
388
0x7ac2,0xa93f,0x9269,0x4082,
 
389
0xef7f,0xbe58,0xc160,0x410c,
 
390
0xacee,0xa9c8,0x84ef,0x4191,
 
391
0x7b83,0x906b,0x78c3,0x4211,
 
392
0x9316,0xfda9,0x3f5e,0x428c,
 
393
0x1e4e,0xd71d,0xa326,0x4301,
 
394
0xa488,0xc547,0x83e3,0x436e,
 
395
0x747f,0x90f6,0x90f1,0x43cb,
 
396
};
 
397
#endif
 
398
#ifdef MIEEE
 
399
static unsigned short YP[24] = {
 
400
0x41d2,0xd2be,0x62f9,0xb6c5,
 
401
0xc262,0xd72d,0x5883,0x6521,
 
402
0x42da,0x0954,0xb091,0x0fef,
 
403
0xc33c,0xe01a,0x37a1,0xb083,
 
404
0x4386,0x79ad,0x0b73,0x66b1,
 
405
0xc3a5,0x9e41,0x5dde,0x2b17,
 
406
};
 
407
static unsigned short YQ[32] = {
 
408
/*0x3ff0,0x0000,0x0000,0x0000,*/
 
409
0x4082,0x9269,0xa93f,0x7ac2,
 
410
0x410c,0xc160,0xbe58,0xef7f,
 
411
0x4191,0x84ef,0xa9c8,0xacee,
 
412
0x4211,0x78c3,0x906b,0x7b83,
 
413
0x428c,0x3f5e,0xfda9,0x9316,
 
414
0x4301,0xa326,0xd71d,0x1e4e,
 
415
0x436e,0x83e3,0xc547,0xa488,
 
416
0x43cb,0x90f1,0x90f6,0x747f,
 
417
};
 
418
#endif
 
419
 
 
420
 
 
421
#ifdef UNK
 
422
static double Z1 = 1.46819706421238932572E1;
 
423
static double Z2 = 4.92184563216946036703E1;
 
424
#endif
 
425
 
 
426
#ifdef DEC
 
427
static unsigned short DZ1[] = {0041152,0164532,0006114,0010540};
 
428
static unsigned short DZ2[] = {0041504,0157663,0001625,0020621};
 
429
#define Z1 (*(double *)DZ1)
 
430
#define Z2 (*(double *)DZ2)
 
431
#endif
 
432
 
 
433
#ifdef IBMPC
 
434
static unsigned short DZ1[] = {0x822c,0x4189,0x5d2b,0x402d};
 
435
static unsigned short DZ2[] = {0xa432,0x6072,0x9bf6,0x4048};
 
436
#define Z1 (*(double *)DZ1)
 
437
#define Z2 (*(double *)DZ2)
 
438
#endif
 
439
 
 
440
#ifdef MIEEE
 
441
static unsigned short DZ1[] = {0x402d,0x5d2b,0x4189,0x822c};
 
442
static unsigned short DZ2[] = {0x4048,0x9bf6,0x6072,0xa432};
 
443
#define Z1 (*(double *)DZ1)
 
444
#define Z2 (*(double *)DZ2)
 
445
#endif
 
446
 
 
447
#ifdef ANSIPROT
 
448
extern double polevl ( double, void *, int );
 
449
extern double p1evl ( double, void *, int );
 
450
extern double log ( double );
 
451
extern double sin ( double );
 
452
extern double cos ( double );
 
453
extern double sqrt ( double );
 
454
double j1 ( double );
 
455
#else
 
456
double polevl(), p1evl(), log(), sin(), cos(), sqrt();
 
457
double j1();
 
458
#endif
 
459
extern double TWOOPI, THPIO4, SQ2OPI;
 
460
 
 
461
double j1(x)
 
462
double x;
 
463
{
 
464
double w, z, p, q, xn;
 
465
 
 
466
w = x;
 
467
if( x < 0 )
 
468
        return -j1(-x);
 
469
 
 
470
if( w <= 5.0 )
 
471
        {
 
472
        z = x * x;      
 
473
        w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 );
 
474
        w = w * x * (z - Z1) * (z - Z2);
 
475
        return( w );
 
476
        }
 
477
 
 
478
w = 5.0/x;
 
479
z = w * w;
 
480
p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
 
481
q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
 
482
xn = x - THPIO4;
 
483
p = p * cos(xn) - w * q * sin(xn);
 
484
return( p * SQ2OPI / sqrt(x) );
 
485
}
 
486
 
 
487
 
 
488
extern double INFINITY, NAN;
 
489
 
 
490
double y1(x)
 
491
double x;
 
492
{
 
493
double w, z, p, q, xn;
 
494
 
 
495
if( x <= 5.0 )
 
496
        {
 
497
        if (x == 0.0) {
 
498
                mtherr("y1", SING);
 
499
                return -INFINITY;
 
500
        } else if (x <= 0.0) {
 
501
                mtherr("y1", DOMAIN);
 
502
                return NAN;
 
503
        }
 
504
        z = x * x;
 
505
        w = x * (polevl( z, YP, 5 ) / p1evl( z, YQ, 8 ));
 
506
        w += TWOOPI * ( j1(x) * log(x)  -  1.0/x );
 
507
        return( w );
 
508
        }
 
509
 
 
510
w = 5.0/x;
 
511
z = w * w;
 
512
p = polevl( z, PP, 6)/polevl( z, PQ, 6 );
 
513
q = polevl( z, QP, 7)/p1evl( z, QQ, 7 );
 
514
xn = x - THPIO4;
 
515
p = p * sin(xn) + w * q * cos(xn);
 
516
return( p * SQ2OPI / sqrt(x) );
 
517
}