← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to Lib/special/cephes/fresnl.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:
Lib/special/cephes/fresnl.c
1
 
/*                                                      fresnl.c
2
 
 *
3
 
 *      Fresnel integral
4
 
 *
5
 
 *
6
 
 *
7
 
 * SYNOPSIS:
8
 
 *
9
 
 * double x, S, C;
10
 
 * void fresnl();
11
 
 *
12
 
 * fresnl( x, _&S, _&C );
13
 
 *
14
 
 *
15
 
 * DESCRIPTION:
16
 
 *
17
 
 * Evaluates the Fresnel integrals
18
 
 *
19
 
 *           x
20
 
 *           -
21
 
 *          | |
22
 
 * C(x) =   |   cos(pi/2 t**2) dt,
23
 
 *        | |
24
 
 *         -
25
 
 *          0
26
 
 *
27
 
 *           x
28
 
 *           -
29
 
 *          | |
30
 
 * S(x) =   |   sin(pi/2 t**2) dt.
31
 
 *        | |
32
 
 *         -
33
 
 *          0
34
 
 *
35
 
 *
36
 
 * The integrals are evaluated by a power series for x < 1.
37
 
 * For x >= 1 auxiliary functions f(x) and g(x) are employed
38
 
 * such that
39
 
 *
40
 
 * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
41
 
 * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
42
 
 *
43
 
 *
44
 
 *
45
 
 * ACCURACY:
46
 
 *
47
 
 *  Relative error.
48
 
 *
49
 
 * Arithmetic  function   domain     # trials      peak         rms
50
 
 *   IEEE       S(x)      0, 10       10000       2.0e-15     3.2e-16
51
 
 *   IEEE       C(x)      0, 10       10000       1.8e-15     3.3e-16
52
 
 *   DEC        S(x)      0, 10        6000       2.2e-16     3.9e-17
53
 
 *   DEC        C(x)      0, 10        5000       2.3e-16     3.9e-17
54
 
 */
55
 
 
56
 
/*
57
 
Cephes Math Library Release 2.1:  January, 1989
58
 
Copyright 1984, 1987, 1989 by Stephen L. Moshier
59
 
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
60
 
*/
61
 
 
62
 
#include "mconf.h"
63
 
 
64
 
/* S(x) for small x */
65
 
#ifdef UNK
66
 
static double sn[6] = {
67
 
-2.99181919401019853726E3,
68
 
 7.08840045257738576863E5,
69
 
-6.29741486205862506537E7,
70
 
 2.54890880573376359104E9,
71
 
-4.42979518059697779103E10,
72
 
 3.18016297876567817986E11,
73
 
};
74
 
static double sd[6] = {
75
 
/* 1.00000000000000000000E0,*/
76
 
 2.81376268889994315696E2,
77
 
 4.55847810806532581675E4,
78
 
 5.17343888770096400730E6,
79
 
 4.19320245898111231129E8,
80
 
 2.24411795645340920940E10,
81
 
 6.07366389490084639049E11,
82
 
};
83
 
#endif
84
 
#ifdef DEC
85
 
static unsigned short sn[24] = {
86
 
0143072,0176433,0065455,0127034,
87
 
0045055,0007200,0134540,0026661,
88
 
0146560,0035061,0023667,0127545,
89
 
0050027,0166503,0002673,0153756,
90
 
0151045,0002721,0121737,0102066,
91
 
0051624,0013177,0033451,0021271,
92
 
};
93
 
static unsigned short sd[24] = {
94
 
/*0040200,0000000,0000000,0000000,*/
95
 
0042214,0130051,0112070,0101617,
96
 
0044062,0010307,0172346,0152510,
97
 
0045635,0160575,0143200,0136642,
98
 
0047307,0171215,0127457,0052361,
99
 
0050647,0031447,0032621,0013510,
100
 
0052015,0064733,0117362,0012653,
101
 
};
102
 
#endif
103
 
#ifdef IBMPC
104
 
static unsigned short sn[24] = {
105
 
0xb5c3,0x6d65,0x5fa3,0xc0a7,
106
 
0x05b6,0x172c,0xa1d0,0x4125,
107
 
0xf5ed,0x24f6,0x0746,0xc18e,
108
 
0x7afe,0x60b7,0xfda8,0x41e2,
109
 
0xf087,0x347b,0xa0ba,0xc224,
110
 
0x2457,0xe6e5,0x82cf,0x4252,
111
 
};
112
 
static unsigned short sd[24] = {
113
 
/*0x0000,0x0000,0x0000,0x3ff0,*/
114
 
0x1072,0x3287,0x9605,0x4071,
115
 
0xdaa9,0xfe9c,0x4218,0x40e6,
116
 
0x17b4,0xb8d0,0xbc2f,0x4153,
117
 
0xea9e,0xb5e5,0xfe51,0x41b8,
118
 
0x22e9,0xe6b2,0xe664,0x4214,
119
 
0x42b5,0x73de,0xad3b,0x4261,
120
 
};
121
 
#endif
122
 
#ifdef MIEEE
123
 
static unsigned short sn[24] = {
124
 
0xc0a7,0x5fa3,0x6d65,0xb5c3,
125
 
0x4125,0xa1d0,0x172c,0x05b6,
126
 
0xc18e,0x0746,0x24f6,0xf5ed,
127
 
0x41e2,0xfda8,0x60b7,0x7afe,
128
 
0xc224,0xa0ba,0x347b,0xf087,
129
 
0x4252,0x82cf,0xe6e5,0x2457,
130
 
};
131
 
static unsigned short sd[24] = {
132
 
/*0x3ff0,0x0000,0x0000,0x0000,*/
133
 
0x4071,0x9605,0x3287,0x1072,
134
 
0x40e6,0x4218,0xfe9c,0xdaa9,
135
 
0x4153,0xbc2f,0xb8d0,0x17b4,
136
 
0x41b8,0xfe51,0xb5e5,0xea9e,
137
 
0x4214,0xe664,0xe6b2,0x22e9,
138
 
0x4261,0xad3b,0x73de,0x42b5,
139
 
};
140
 
#endif
141
 
 
142
 
/* C(x) for small x */
143
 
#ifdef UNK
144
 
static double cn[6] = {
145
 
-4.98843114573573548651E-8,
146
 
 9.50428062829859605134E-6,
147
 
-6.45191435683965050962E-4,
148
 
 1.88843319396703850064E-2,
149
 
-2.05525900955013891793E-1,
150
 
 9.99999999999999998822E-1,
151
 
};
152
 
static double cd[7] = {
153
 
 3.99982968972495980367E-12,
154
 
 9.15439215774657478799E-10,
155
 
 1.25001862479598821474E-7,
156
 
 1.22262789024179030997E-5,
157
 
 8.68029542941784300606E-4,
158
 
 4.12142090722199792936E-2,
159
 
 1.00000000000000000118E0,
160
 
};
161
 
#endif
162
 
#ifdef DEC
163
 
static unsigned short cn[24] = {
164
 
0132126,0040141,0063733,0013231,
165
 
0034037,0072223,0010200,0075637,
166
 
0135451,0021020,0073264,0036057,
167
 
0036632,0131520,0101316,0060233,
168
 
0137522,0072541,0136124,0132202,
169
 
0040200,0000000,0000000,0000000,
170
 
};
171
 
static unsigned short cd[28] = {
172
 
0026614,0135503,0051776,0032631,
173
 
0030573,0121116,0154033,0126712,
174
 
0032406,0034100,0012442,0106212,
175
 
0034115,0017567,0150520,0164623,
176
 
0035543,0106171,0177336,0146351,
177
 
0037050,0150073,0000607,0171635,
178
 
0040200,0000000,0000000,0000000,
179
 
};
180
 
#endif
181
 
#ifdef IBMPC
182
 
static unsigned short cn[24] = {
183
 
0x62d3,0x2cfb,0xc80c,0xbe6a,
184
 
0x0f74,0x6210,0xee92,0x3ee3,
185
 
0x8786,0x0ed6,0x2442,0xbf45,
186
 
0xcc13,0x1059,0x566a,0x3f93,
187
 
0x9690,0x378a,0x4eac,0xbfca,
188
 
0x0000,0x0000,0x0000,0x3ff0,
189
 
};
190
 
static unsigned short cd[28] = {
191
 
0xc6b3,0x6a7f,0x9768,0x3d91,
192
 
0x75b9,0xdb03,0x7449,0x3e0f,
193
 
0x5191,0x02a4,0xc708,0x3e80,
194
 
0x1d32,0xfa2a,0xa3ee,0x3ee9,
195
 
0xd99d,0x3fdb,0x718f,0x3f4c,
196
 
0xfe74,0x6030,0x1a07,0x3fa5,
197
 
0x0000,0x0000,0x0000,0x3ff0,
198
 
};
199
 
#endif
200
 
#ifdef MIEEE
201
 
static unsigned short cn[24] = {
202
 
0xbe6a,0xc80c,0x2cfb,0x62d3,
203
 
0x3ee3,0xee92,0x6210,0x0f74,
204
 
0xbf45,0x2442,0x0ed6,0x8786,
205
 
0x3f93,0x566a,0x1059,0xcc13,
206
 
0xbfca,0x4eac,0x378a,0x9690,
207
 
0x3ff0,0x0000,0x0000,0x0000,
208
 
};
209
 
static unsigned short cd[28] = {
210
 
0x3d91,0x9768,0x6a7f,0xc6b3,
211
 
0x3e0f,0x7449,0xdb03,0x75b9,
212
 
0x3e80,0xc708,0x02a4,0x5191,
213
 
0x3ee9,0xa3ee,0xfa2a,0x1d32,
214
 
0x3f4c,0x718f,0x3fdb,0xd99d,
215
 
0x3fa5,0x1a07,0x6030,0xfe74,
216
 
0x3ff0,0x0000,0x0000,0x0000,
217
 
};
218
 
#endif
219
 
 
220
 
/* Auxiliary function f(x) */
221
 
#ifdef UNK
222
 
static double fn[10] = {
223
 
  4.21543555043677546506E-1,
224
 
  1.43407919780758885261E-1,
225
 
  1.15220955073585758835E-2,
226
 
  3.45017939782574027900E-4,
227
 
  4.63613749287867322088E-6,
228
 
  3.05568983790257605827E-8,
229
 
  1.02304514164907233465E-10,
230
 
  1.72010743268161828879E-13,
231
 
  1.34283276233062758925E-16,
232
 
  3.76329711269987889006E-20,
233
 
};
234
 
static double fd[10] = {
235
 
/*  1.00000000000000000000E0,*/
236
 
  7.51586398353378947175E-1,
237
 
  1.16888925859191382142E-1,
238
 
  6.44051526508858611005E-3,
239
 
  1.55934409164153020873E-4,
240
 
  1.84627567348930545870E-6,
241
 
  1.12699224763999035261E-8,
242
 
  3.60140029589371370404E-11,
243
 
  5.88754533621578410010E-14,
244
 
  4.52001434074129701496E-17,
245
 
  1.25443237090011264384E-20,
246
 
};
247
 
#endif
248
 
#ifdef DEC
249
 
static unsigned short fn[40] = {
250
 
0037727,0152216,0106601,0016214,
251
 
0037422,0154606,0112710,0071355,
252
 
0036474,0143453,0154253,0166545,
253
 
0035264,0161606,0022250,0073743,
254
 
0033633,0110036,0024653,0136246,
255
 
0032003,0036652,0041164,0036413,
256
 
0027740,0174122,0046305,0036726,
257
 
0025501,0125270,0121317,0167667,
258
 
0023032,0150555,0076175,0047443,
259
 
0020061,0133570,0070130,0027657,
260
 
};
261
 
static unsigned short fd[40] = {
262
 
/*0040200,0000000,0000000,0000000,*/
263
 
0040100,0063767,0054413,0151452,
264
 
0037357,0061566,0007243,0065754,
265
 
0036323,0005365,0033552,0133625,
266
 
0035043,0101123,0000275,0165402,
267
 
0033367,0146614,0110623,0023647,
268
 
0031501,0116644,0125222,0144263,
269
 
0027436,0062051,0117235,0001411,
270
 
0025204,0111543,0056370,0036201,
271
 
0022520,0071351,0015227,0122144,
272
 
0017554,0172240,0112713,0005006,
273
 
};
274
 
#endif
275
 
#ifdef IBMPC
276
 
static unsigned short fn[40] = {
277
 
0x2391,0xd1b0,0xfa91,0x3fda,
278
 
0x0e5e,0xd2b9,0x5b30,0x3fc2,
279
 
0x7dad,0x7b15,0x98e5,0x3f87,
280
 
0x0efc,0xc495,0x9c70,0x3f36,
281
 
0x7795,0xc535,0x7203,0x3ed3,
282
 
0x87a1,0x484e,0x67b5,0x3e60,
283
 
0xa7bb,0x4998,0x1f0a,0x3ddc,
284
 
0xfdf7,0x1459,0x3557,0x3d48,
285
 
0xa9e4,0xaf8f,0x5a2d,0x3ca3,
286
 
0x05f6,0x0e0b,0x36ef,0x3be6,
287
 
};
288
 
static unsigned short fd[40] = {
289
 
/*0x0000,0x0000,0x0000,0x3ff0,*/
290
 
0x7a65,0xeb21,0x0cfe,0x3fe8,
291
 
0x6d7d,0xc1d4,0xec6e,0x3fbd,
292
 
0x56f3,0xa6ed,0x615e,0x3f7a,
293
 
0xbd60,0x6017,0x704a,0x3f24,
294
 
0x64f5,0x9232,0xf9b1,0x3ebe,
295
 
0x5916,0x9552,0x33b4,0x3e48,
296
 
0xa061,0x33d3,0xcc85,0x3dc3,
297
 
0x0790,0x6b9f,0x926c,0x3d30,
298
 
0xf48d,0x2352,0x0e5d,0x3c8a,
299
 
0x6141,0x12b9,0x9e94,0x3bcd,
300
 
};
301
 
#endif
302
 
#ifdef MIEEE
303
 
static unsigned short fn[40] = {
304
 
0x3fda,0xfa91,0xd1b0,0x2391,
305
 
0x3fc2,0x5b30,0xd2b9,0x0e5e,
306
 
0x3f87,0x98e5,0x7b15,0x7dad,
307
 
0x3f36,0x9c70,0xc495,0x0efc,
308
 
0x3ed3,0x7203,0xc535,0x7795,
309
 
0x3e60,0x67b5,0x484e,0x87a1,
310
 
0x3ddc,0x1f0a,0x4998,0xa7bb,
311
 
0x3d48,0x3557,0x1459,0xfdf7,
312
 
0x3ca3,0x5a2d,0xaf8f,0xa9e4,
313
 
0x3be6,0x36ef,0x0e0b,0x05f6,
314
 
};
315
 
static unsigned short fd[40] = {
316
 
/*0x3ff0,0x0000,0x0000,0x0000,*/
317
 
0x3fe8,0x0cfe,0xeb21,0x7a65,
318
 
0x3fbd,0xec6e,0xc1d4,0x6d7d,
319
 
0x3f7a,0x615e,0xa6ed,0x56f3,
320
 
0x3f24,0x704a,0x6017,0xbd60,
321
 
0x3ebe,0xf9b1,0x9232,0x64f5,
322
 
0x3e48,0x33b4,0x9552,0x5916,
323
 
0x3dc3,0xcc85,0x33d3,0xa061,
324
 
0x3d30,0x926c,0x6b9f,0x0790,
325
 
0x3c8a,0x0e5d,0x2352,0xf48d,
326
 
0x3bcd,0x9e94,0x12b9,0x6141,
327
 
};
328
 
#endif
329
 
 
330
 
 
331
 
/* Auxiliary function g(x) */
332
 
#ifdef UNK
333
 
static double gn[11] = {
334
 
  5.04442073643383265887E-1,
335
 
  1.97102833525523411709E-1,
336
 
  1.87648584092575249293E-2,
337
 
  6.84079380915393090172E-4,
338
 
  1.15138826111884280931E-5,
339
 
  9.82852443688422223854E-8,
340
 
  4.45344415861750144738E-10,
341
 
  1.08268041139020870318E-12,
342
 
  1.37555460633261799868E-15,
343
 
  8.36354435630677421531E-19,
344
 
  1.86958710162783235106E-22,
345
 
};
346
 
static double gd[11] = {
347
 
/*  1.00000000000000000000E0,*/
348
 
  1.47495759925128324529E0,
349
 
  3.37748989120019970451E-1,
350
 
  2.53603741420338795122E-2,
351
 
  8.14679107184306179049E-4,
352
 
  1.27545075667729118702E-5,
353
 
  1.04314589657571990585E-7,
354
 
  4.60680728146520428211E-10,
355
 
  1.10273215066240270757E-12,
356
 
  1.38796531259578871258E-15,
357
 
  8.39158816283118707363E-19,
358
 
  1.86958710162783236342E-22,
359
 
};
360
 
#endif
361
 
#ifdef DEC
362
 
static unsigned short gn[44] = {
363
 
0040001,0021435,0120406,0053123,
364
 
0037511,0152523,0037703,0122011,
365
 
0036631,0134302,0122721,0110235,
366
 
0035463,0051712,0043215,0114732,
367
 
0034101,0025677,0147725,0057630,
368
 
0032323,0010342,0067523,0002206,
369
 
0030364,0152247,0110007,0054107,
370
 
0026230,0057654,0035464,0047124,
371
 
0023706,0036401,0167705,0045440,
372
 
0021166,0154447,0105632,0142461,
373
 
0016142,0002353,0011175,0170530,
374
 
};
375
 
static unsigned short gd[44] = {
376
 
/*0040200,0000000,0000000,0000000,*/
377
 
0040274,0145551,0016742,0127005,
378
 
0037654,0166557,0076416,0015165,
379
 
0036717,0140217,0030675,0050111,
380
 
0035525,0110060,0076405,0070502,
381
 
0034125,0176061,0060120,0031730,
382
 
0032340,0001615,0054343,0120501,
383
 
0030375,0041414,0070747,0107060,
384
 
0026233,0031034,0160757,0074526,
385
 
0023710,0003341,0137100,0144664,
386
 
0021167,0126414,0023774,0015435,
387
 
0016142,0002353,0011175,0170530,
388
 
};
389
 
#endif
390
 
#ifdef IBMPC
391
 
static unsigned short gn[44] = {
392
 
0xcaca,0xb420,0x2463,0x3fe0,
393
 
0x7481,0x67f8,0x3aaa,0x3fc9,
394
 
0x3214,0x54ba,0x3718,0x3f93,
395
 
0xb33b,0x48d1,0x6a79,0x3f46,
396
 
0xabf3,0xf9fa,0x2577,0x3ee8,
397
 
0x6091,0x4dea,0x621c,0x3e7a,
398
 
0xeb09,0xf200,0x9a94,0x3dfe,
399
 
0x89cb,0x8766,0x0bf5,0x3d73,
400
 
0xa964,0x3df8,0xc7a0,0x3cd8,
401
 
0x58a6,0xf173,0xdb24,0x3c2e,
402
 
0xbe2b,0x624f,0x409d,0x3b6c,
403
 
};
404
 
static unsigned short gd[44] = {
405
 
/*0x0000,0x0000,0x0000,0x3ff0,*/
406
 
0x55c1,0x23bc,0x996d,0x3ff7,
407
 
0xc34f,0xefa1,0x9dad,0x3fd5,
408
 
0xaa09,0xe637,0xf811,0x3f99,
409
 
0xae28,0x0fa0,0xb206,0x3f4a,
410
 
0x067b,0x2c0a,0xbf86,0x3eea,
411
 
0x7428,0xab1c,0x0071,0x3e7c,
412
 
0xf1c6,0x8e3c,0xa861,0x3dff,
413
 
0xef2b,0x9c3d,0x6643,0x3d73,
414
 
0x1936,0x37c8,0x00dc,0x3cd9,
415
 
0x8364,0x84ff,0xf5a1,0x3c2e,
416
 
0xbe2b,0x624f,0x409d,0x3b6c,
417
 
};
418
 
#endif
419
 
#ifdef MIEEE
420
 
static unsigned short gn[44] = {
421
 
0x3fe0,0x2463,0xb420,0xcaca,
422
 
0x3fc9,0x3aaa,0x67f8,0x7481,
423
 
0x3f93,0x3718,0x54ba,0x3214,
424
 
0x3f46,0x6a79,0x48d1,0xb33b,
425
 
0x3ee8,0x2577,0xf9fa,0xabf3,
426
 
0x3e7a,0x621c,0x4dea,0x6091,
427
 
0x3dfe,0x9a94,0xf200,0xeb09,
428
 
0x3d73,0x0bf5,0x8766,0x89cb,
429
 
0x3cd8,0xc7a0,0x3df8,0xa964,
430
 
0x3c2e,0xdb24,0xf173,0x58a6,
431
 
0x3b6c,0x409d,0x624f,0xbe2b,
432
 
};
433
 
static unsigned short gd[44] = {
434
 
/*0x3ff0,0x0000,0x0000,0x0000,*/
435
 
0x3ff7,0x996d,0x23bc,0x55c1,
436
 
0x3fd5,0x9dad,0xefa1,0xc34f,
437
 
0x3f99,0xf811,0xe637,0xaa09,
438
 
0x3f4a,0xb206,0x0fa0,0xae28,
439
 
0x3eea,0xbf86,0x2c0a,0x067b,
440
 
0x3e7c,0x0071,0xab1c,0x7428,
441
 
0x3dff,0xa861,0x8e3c,0xf1c6,
442
 
0x3d73,0x6643,0x9c3d,0xef2b,
443
 
0x3cd9,0x00dc,0x37c8,0x1936,
444
 
0x3c2e,0xf5a1,0x84ff,0x8364,
445
 
0x3b6c,0x409d,0x624f,0xbe2b,
446
 
};
447
 
#endif
448
 
 
449
 
#ifndef ANSIPROT
450
 
double fabs(), cos(), sin(), polevl(), p1evl();
451
 
#endif
452
 
extern double PI, PIO2, MACHEP;
453
 
 
454
 
int fresnl( xxa, ssa, cca )
455
 
double xxa, *ssa, *cca;
456
 
{
457
 
double f, g, cc, ss, c, s, t, u;
458
 
double x, x2;
459
 
 
460
 
x = fabs(xxa);
461
 
x2 = x * x;
462
 
if( x2 < 2.5625 )
463
 
        {
464
 
        t = x2 * x2;
465
 
        ss = x * x2 * polevl( t, sn, 5)/p1evl( t, sd, 6 );
466
 
        cc = x * polevl( t, cn, 5)/polevl(t, cd, 6 );
467
 
        goto done;
468
 
        }
469
 
 
470
 
 
471
 
 
472
 
 
473
 
 
474
 
 
475
 
if( x > 36974.0 )
476
 
        {
477
 
        cc = 0.5;
478
 
        ss = 0.5;
479
 
        goto done;
480
 
        }
481
 
 
482
 
 
483
 
/*              Asymptotic power series auxiliary functions
484
 
 *              for large argument
485
 
 */
486
 
        x2 = x * x;
487
 
        t = PI * x2;
488
 
        u = 1.0/(t * t);
489
 
        t = 1.0/t;
490
 
        f = 1.0 - u * polevl( u, fn, 9)/p1evl(u, fd, 10);
491
 
        g = t * polevl( u, gn, 10)/p1evl(u, gd, 11);
492
 
 
493
 
        t = PIO2 * x2;
494
 
        c = cos(t);
495
 
        s = sin(t);
496
 
        t = PI * x;
497
 
        cc = 0.5  +  (f * s  -  g * c)/t;
498
 
        ss = 0.5  -  (f * c  +  g * s)/t;
499
 
 
500
 
done:
501
 
if( xxa < 0.0 )
502
 
        {
503
 
        cc = -cc;
504
 
        ss = -ss;
505
 
        }
506
 
 
507
 
*cca = cc;
508
 
*ssa = ss;
509
 
return(0);
510
 
}