← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/cephes/k1.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/k1.c
 
1
/*                                                      k1.c
 
2
 *
 
3
 *      Modified Bessel function, third kind, order one
 
4
 *
 
5
 *
 
6
 *
 
7
 * SYNOPSIS:
 
8
 *
 
9
 * double x, y, k1();
 
10
 *
 
11
 * y = k1( x );
 
12
 *
 
13
 *
 
14
 *
 
15
 * DESCRIPTION:
 
16
 *
 
17
 * Computes the modified Bessel function of the third kind
 
18
 * of order one of the argument.
 
19
 *
 
20
 * The range is partitioned into the two intervals [0,2] and
 
21
 * (2, infinity).  Chebyshev polynomial expansions are employed
 
22
 * in each interval.
 
23
 *
 
24
 *
 
25
 *
 
26
 * ACCURACY:
 
27
 *
 
28
 *                      Relative error:
 
29
 * arithmetic   domain     # trials      peak         rms
 
30
 *    DEC       0, 30        3300       8.9e-17     2.2e-17
 
31
 *    IEEE      0, 30       30000       1.2e-15     1.6e-16
 
32
 *
 
33
 * ERROR MESSAGES:
 
34
 *
 
35
 *   message         condition      value returned
 
36
 * k1 domain          x <= 0          MAXNUM
 
37
 *
 
38
 */
 
39
/*                                                     k1e.c
 
40
 *
 
41
 *      Modified Bessel function, third kind, order one,
 
42
 *      exponentially scaled
 
43
 *
 
44
 *
 
45
 *
 
46
 * SYNOPSIS:
 
47
 *
 
48
 * double x, y, k1e();
 
49
 *
 
50
 * y = k1e( x );
 
51
 *
 
52
 *
 
53
 *
 
54
 * DESCRIPTION:
 
55
 *
 
56
 * Returns exponentially scaled modified Bessel function
 
57
 * of the third kind of order one of the argument:
 
58
 *
 
59
 *      k1e(x) = exp(x) * k1(x).
 
60
 *
 
61
 *
 
62
 *
 
63
 * ACCURACY:
 
64
 *
 
65
 *                      Relative error:
 
66
 * arithmetic   domain     # trials      peak         rms
 
67
 *    IEEE      0, 30       30000       7.8e-16     1.2e-16
 
68
 * See k1().
 
69
 *
 
70
 */
 
71
 
 
72
/*
 
73
Cephes Math Library Release 2.0:  April, 1987
 
74
Copyright 1984, 1987 by Stephen L. Moshier
 
75
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
76
*/
 
77
 
 
78
#include "mconf.h"
 
79
 
 
80
/* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
 
81
 * in the interval [0,2].
 
82
 * 
 
83
 * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
 
84
 */
 
85
 
 
86
#ifdef UNK
 
87
static double A[] =
 
88
{
 
89
-7.02386347938628759343E-18,
 
90
-2.42744985051936593393E-15,
 
91
-6.66690169419932900609E-13,
 
92
-1.41148839263352776110E-10,
 
93
-2.21338763073472585583E-8,
 
94
-2.43340614156596823496E-6,
 
95
-1.73028895751305206302E-4,
 
96
-6.97572385963986435018E-3,
 
97
-1.22611180822657148235E-1,
 
98
-3.53155960776544875667E-1,
 
99
 1.52530022733894777053E0
 
100
};
 
101
#endif
 
102
 
 
103
#ifdef DEC
 
104
static unsigned short A[] = {
 
105
0122001,0110501,0164746,0151255,
 
106
0124056,0165213,0150034,0147377,
 
107
0126073,0124026,0167207,0001044,
 
108
0130033,0030735,0141061,0033116,
 
109
0131676,0020350,0121341,0107175,
 
110
0133443,0046631,0062031,0070716,
 
111
0135065,0067427,0026435,0164022,
 
112
0136344,0112234,0165752,0006222,
 
113
0137373,0015622,0017016,0155636,
 
114
0137664,0150333,0125730,0067240,
 
115
0040303,0036411,0130200,0043120
 
116
};
 
117
#endif
 
118
 
 
119
#ifdef IBMPC
 
120
static unsigned short A[] = {
 
121
0xda56,0x3d3c,0x3228,0xbc60,
 
122
0x99e0,0x7a03,0xdd51,0xbce5,
 
123
0xe045,0xddd0,0x7502,0xbd67,
 
124
0x26ca,0xb846,0x663b,0xbde3,
 
125
0x31d0,0x145c,0xc41d,0xbe57,
 
126
0x2e3a,0x2c83,0x69b3,0xbec4,
 
127
0xbd02,0xe5a3,0xade2,0xbf26,
 
128
0x4192,0x9d7d,0x9293,0xbf7c,
 
129
0xdb74,0x43c1,0x6372,0xbfbf,
 
130
0x0dd4,0x757b,0x9a1b,0xbfd6,
 
131
0x08ca,0x3610,0x67a1,0x3ff8
 
132
};
 
133
#endif
 
134
 
 
135
#ifdef MIEEE
 
136
static unsigned short A[] = {
 
137
0xbc60,0x3228,0x3d3c,0xda56,
 
138
0xbce5,0xdd51,0x7a03,0x99e0,
 
139
0xbd67,0x7502,0xddd0,0xe045,
 
140
0xbde3,0x663b,0xb846,0x26ca,
 
141
0xbe57,0xc41d,0x145c,0x31d0,
 
142
0xbec4,0x69b3,0x2c83,0x2e3a,
 
143
0xbf26,0xade2,0xe5a3,0xbd02,
 
144
0xbf7c,0x9293,0x9d7d,0x4192,
 
145
0xbfbf,0x6372,0x43c1,0xdb74,
 
146
0xbfd6,0x9a1b,0x757b,0x0dd4,
 
147
0x3ff8,0x67a1,0x3610,0x08ca
 
148
};
 
149
#endif
 
150
 
 
151
 
 
152
 
 
153
/* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
 
154
 * in the interval [2,infinity].
 
155
 *
 
156
 * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
 
157
 */
 
158
 
 
159
#ifdef UNK
 
160
static double B[] =
 
161
{
 
162
-5.75674448366501715755E-18,
 
163
 1.79405087314755922667E-17,
 
164
-5.68946255844285935196E-17,
 
165
 1.83809354436663880070E-16,
 
166
-6.05704724837331885336E-16,
 
167
 2.03870316562433424052E-15,
 
168
-7.01983709041831346144E-15,
 
169
 2.47715442448130437068E-14,
 
170
-8.97670518232499435011E-14,
 
171
 3.34841966607842919884E-13,
 
172
-1.28917396095102890680E-12,
 
173
 5.13963967348173025100E-12,
 
174
-2.12996783842756842877E-11,
 
175
 9.21831518760500529508E-11,
 
176
-4.19035475934189648750E-10,
 
177
 2.01504975519703286596E-9,
 
178
-1.03457624656780970260E-8,
 
179
 5.74108412545004946722E-8,
 
180
-3.50196060308781257119E-7,
 
181
 2.40648494783721712015E-6,
 
182
-1.93619797416608296024E-5,
 
183
 1.95215518471351631108E-4,
 
184
-2.85781685962277938680E-3,
 
185
 1.03923736576817238437E-1,
 
186
 2.72062619048444266945E0
 
187
};
 
188
#endif
 
189
 
 
190
#ifdef DEC
 
191
static unsigned short B[] = {
 
192
0121724,0061352,0013041,0150076,
 
193
0022245,0074324,0016172,0173232,
 
194
0122603,0030250,0135670,0165221,
 
195
0023123,0165362,0023561,0060124,
 
196
0123456,0112436,0141654,0073623,
 
197
0024022,0163557,0077564,0006753,
 
198
0124374,0165221,0131014,0026524,
 
199
0024737,0017512,0144250,0175451,
 
200
0125312,0021456,0123136,0076633,
 
201
0025674,0077720,0020125,0102607,
 
202
0126265,0067543,0007744,0043701,
 
203
0026664,0152702,0033002,0074202,
 
204
0127273,0055234,0120016,0071733,
 
205
0027712,0133200,0042441,0075515,
 
206
0130346,0057000,0015456,0074470,
 
207
0031012,0074441,0051636,0111155,
 
208
0131461,0136444,0177417,0002101,
 
209
0032166,0111743,0032176,0021410,
 
210
0132674,0001224,0076555,0027060,
 
211
0033441,0077430,0135226,0106663,
 
212
0134242,0065610,0167155,0113447,
 
213
0035114,0131304,0043664,0102163,
 
214
0136073,0045065,0171465,0122123,
 
215
0037324,0152767,0147401,0017732,
 
216
0040456,0017275,0050061,0062120,
 
217
};
 
218
#endif
 
219
 
 
220
#ifdef IBMPC
 
221
static unsigned short B[] = {
 
222
0x3a08,0x42c4,0x8c5d,0xbc5a,
 
223
0x5ed3,0x838f,0xaf1a,0x3c74,
 
224
0x1d52,0x1777,0x6615,0xbc90,
 
225
0x2c0b,0x44ee,0x7d5e,0x3caa,
 
226
0x8ef2,0xd875,0xd2a3,0xbcc5,
 
227
0x81bd,0xefee,0x5ced,0x3ce2,
 
228
0x85ab,0x3641,0x9d52,0xbcff,
 
229
0x1f65,0x5915,0xe3e9,0x3d1b,
 
230
0xcfb3,0xd4cb,0x4465,0xbd39,
 
231
0xb0b1,0x040a,0x8ffa,0x3d57,
 
232
0x88f8,0x61fc,0xadec,0xbd76,
 
233
0x4f10,0x46c0,0x9ab8,0x3d96,
 
234
0xce7b,0x9401,0x6b53,0xbdb7,
 
235
0x2f6a,0x08a4,0x56d0,0x3dd9,
 
236
0xcf27,0x0365,0xcbc0,0xbdfc,
 
237
0xd24e,0x2a73,0x4f24,0x3e21,
 
238
0xe088,0x9fe1,0x37a4,0xbe46,
 
239
0xc461,0x668f,0xd27c,0x3e6e,
 
240
0xa5c6,0x8fad,0x8052,0xbe97,
 
241
0xd1b6,0x1752,0x2fe3,0x3ec4,
 
242
0xb2e5,0x1dcd,0x4d71,0xbef4,
 
243
0x908e,0x88f6,0x9658,0x3f29,
 
244
0xb48a,0xbe66,0x6946,0xbf67,
 
245
0x23fb,0xf9e0,0x9abe,0x3fba,
 
246
0x2c8a,0xaa06,0xc3d7,0x4005
 
247
};
 
248
#endif
 
249
 
 
250
#ifdef MIEEE
 
251
static unsigned short B[] = {
 
252
0xbc5a,0x8c5d,0x42c4,0x3a08,
 
253
0x3c74,0xaf1a,0x838f,0x5ed3,
 
254
0xbc90,0x6615,0x1777,0x1d52,
 
255
0x3caa,0x7d5e,0x44ee,0x2c0b,
 
256
0xbcc5,0xd2a3,0xd875,0x8ef2,
 
257
0x3ce2,0x5ced,0xefee,0x81bd,
 
258
0xbcff,0x9d52,0x3641,0x85ab,
 
259
0x3d1b,0xe3e9,0x5915,0x1f65,
 
260
0xbd39,0x4465,0xd4cb,0xcfb3,
 
261
0x3d57,0x8ffa,0x040a,0xb0b1,
 
262
0xbd76,0xadec,0x61fc,0x88f8,
 
263
0x3d96,0x9ab8,0x46c0,0x4f10,
 
264
0xbdb7,0x6b53,0x9401,0xce7b,
 
265
0x3dd9,0x56d0,0x08a4,0x2f6a,
 
266
0xbdfc,0xcbc0,0x0365,0xcf27,
 
267
0x3e21,0x4f24,0x2a73,0xd24e,
 
268
0xbe46,0x37a4,0x9fe1,0xe088,
 
269
0x3e6e,0xd27c,0x668f,0xc461,
 
270
0xbe97,0x8052,0x8fad,0xa5c6,
 
271
0x3ec4,0x2fe3,0x1752,0xd1b6,
 
272
0xbef4,0x4d71,0x1dcd,0xb2e5,
 
273
0x3f29,0x9658,0x88f6,0x908e,
 
274
0xbf67,0x6946,0xbe66,0xb48a,
 
275
0x3fba,0x9abe,0xf9e0,0x23fb,
 
276
0x4005,0xc3d7,0xaa06,0x2c8a
 
277
};
 
278
#endif
 
279
 
 
280
#ifndef ANSIPROT
 
281
double chbevl(), exp(), i1(), log(), sqrt();
 
282
#endif
 
283
extern double PI;
 
284
extern double MINLOG, MAXNUM;
 
285
 
 
286
double k1(x)
 
287
double x;
 
288
{
 
289
double y, z;
 
290
 
 
291
z = 0.5 * x;
 
292
if( z <= 0.0 )
 
293
        {
 
294
        mtherr( "k1", DOMAIN );
 
295
        return( MAXNUM );
 
296
        }
 
297
 
 
298
if( x <= 2.0 )
 
299
        {
 
300
        y = x * x - 2.0;
 
301
        y =  log(z) * i1(x)  +  chbevl( y, A, 11 ) / x;
 
302
        return( y );
 
303
        }
 
304
 
 
305
return(  exp(-x) * chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x) );
 
306
}
 
307
 
 
308
 
 
309
 
 
310
 
 
311
double k1e( x )
 
312
double x;
 
313
{
 
314
double y;
 
315
 
 
316
if( x <= 0.0 )
 
317
        {
 
318
        mtherr( "k1e", DOMAIN );
 
319
        return( MAXNUM );
 
320
        }
 
321
 
 
322
if( x <= 2.0 )
 
323
        {
 
324
        y = x * x - 2.0;
 
325
        y =  log( 0.5 * x ) * i1(x)  +  chbevl( y, A, 11 ) / x;
 
326
        return( y * exp(x) );
 
327
        }
 
328
 
 
329
return(  chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x) );
 
330
}