forked from freemint/fdlibm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathe_j0.c
145 lines (135 loc) · 4.25 KB
/
e_j0.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
/* @(#)e_j0.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
for performance improvement on pipelined processors.
*/
/* __ieee754_j0(x)
* Bessel function of the first kind of order zero.
* Method -- j0(x):
* 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
* 2. Reduce x to |x| since j0(x)=j0(-x), and
* for x in (0,2)
* j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
* (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
* for x in (2,inf)
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* as follow:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (cos(x) + sin(x))
* sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
*
* 3 Special cases
* j0(nan)= nan
* j0(0) = 1
* j0(inf) = 0
*/
#ifndef __FDLIBM_H__
#include "fdlibm.h"
#endif
double __ieee754_j0(double x)
{
double z, s, c, ss, cc, r, u, v, r1, r2, s1, s2, z2, z4;
int32_t hx, ix;
static const double hugeval = 1e300;
static const double one = 1.0;
static const double zero = 0.0;
static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
/* R0/S0 on [0, 2.00] */
static const double R02 = 1.56249999999999947958e-02; /* 0x3F8FFFFF, 0xFFFFFFFD */
static const double R03 = -1.89979294238854721751e-04; /* 0xBF28E6A5, 0xB61AC6E9 */
static const double R04 = 1.82954049532700665670e-06; /* 0x3EBEB1D1, 0x0C503919 */
static const double R05 = -4.61832688532103189199e-09; /* 0xBE33D5E7, 0x73D63FCE */
static const double S01 = 1.56191029464890010492e-02; /* 0x3F8FFCE8, 0x82C8C2A4 */
static const double S02 = 1.16926784663337450260e-04; /* 0x3F1EA6D2, 0xDD57DBF4 */
static const double S03 = 5.13546550207318111446e-07; /* 0x3EA13B54, 0xCE84D5A9 */
static const double S04 = 1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
GET_HIGH_WORD(hx, x);
ix = hx & IC(0x7fffffff);
if (ix >= IC(0x7ff00000))
return one / (x * x);
x = __ieee754_fabs(x);
if (ix >= IC(0x40000000))
{ /* |x| >= 2.0 */
__ieee754_sincos(x, &s, &c);
ss = s - c;
cc = s + c;
if (ix < IC(0x7fe00000))
{ /* make sure x+x not overflow */
z = -__ieee754_cos(x + x);
if ((s * c) < zero)
cc = z / ss;
else
ss = z / cc;
}
/*
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if (ix > IC(0x48000000))
z = (invsqrtpi * cc) / __ieee754_sqrt(x);
else
{
u = __j0_y0_pzero(x);
v = __j0_y0_qzero(x);
z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt(x);
}
return z;
}
if (ix < IC(0x3f200000))
{ /* |x| < 2**-13 */
math_force_eval(hugeval + x); /* raise inexact if x != 0 */
if (ix < IC(0x3e400000))
return one; /* |x|<2**-27 */
else
return one - 0.25 * x * x;
}
z = x * x;
#ifdef DO_NOT_USE_THIS
r = z * (R02 + z * (R03 + z * (R04 + z * R05)));
s = one + z * (S01 + z * (S02 + z * (S03 + z * S04)));
#else
r1 = z * R02;
z2 = z * z;
r2 = R03 + z * R04;
z4 = z2 * z2;
r = r1 + z2 * r2 + z4 * R05;
s1 = one + z * S01;
s2 = S02 + z * S03;
s = s1 + z2 * s2 + z4 * S04;
#endif
if (ix < IC(0x3FF00000))
{ /* |x| < 1.00 */
return one + z * (-0.25 + (r / s));
} else
{
u = 0.5 * x;
return ((one + u) * (one - u) + z * (r / s));
}
}
/* wrapper j0 */
double __j0(double x)
{
if (isgreater(__ieee754_fabs(x), X_TLOSS) && _LIB_VERSION != _IEEE_ && _LIB_VERSION != _POSIX_)
/* j0(|x|>X_TLOSS) */
return __kernel_standard(x, x, 0.0, KMATHERR_J0_TLOSS);
return __ieee754_j0(x);
}
__typeof(__j0) j0 __attribute__((weak, alias("__j0")));
#ifdef __NO_LONG_DOUBLE_MATH
__typeof(__j0l) __j0l __attribute__((alias("__j0")));
__typeof(__j0l) j0l __attribute__((weak, alias("__j0l")));
#endif