Print this page
11210 libm should be cstyle(1ONBLD) clean
@@ -16,13 +16,15 @@
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
+
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2005 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
@@ -38,13 +40,11 @@
#if defined(__i386) && !defined(__amd64)
extern int __swapRP(int);
#endif
-static const float
- zerof = 0.0f,
- onef = 1.0f;
+static const float zerof = 0.0f, onef = 1.0f;
static const double C[] = {
0.0,
-0.125,
0.25,
@@ -100,13 +100,13 @@
0.106695157020407986137501682e3,
1.0,
};
#define pr Cj0y0
-#define ps (Cj0y0+7)
-#define qr (Cj0y0+14)
-#define qs (Cj0y0+21)
+#define ps (Cj0y0 + 7)
+#define qr (Cj0y0 + 14)
+#define qs (Cj0y0 + 21)
static const double Cj0[] = {
-2.500000000000003622131880894830476755537e-0001, /* r0 */
1.095597547334830263234433855932375353303e-0002,
-1.819734750463320921799187258987098087697e-0004,
@@ -134,13 +134,13 @@
2.801559820648939665270492520004836611187e-0018,
2.073101088320349159764410261466350732968e-0021,
};
#define r0 Cj0
-#define s0 (Cj0+4)
-#define r1 (Cj0+8)
-#define s1 (Cj0+17)
+#define s0 (Cj0 + 4)
+#define r1 (Cj0 + 8)
+#define s1 (Cj0 + 17)
static const double Cy0[] = {
-7.380429510868722526754723020704317641941e-0002, /* u0 */
1.772607102684869924301459663049874294814e-0001,
-1.524370666542713828604078090970799356306e-0002,
@@ -160,11 +160,11 @@
1.001495929158861646659010844136682454906e-0008,
4.725338116256021660204443235685358593611e-0012,
};
#define u0 Cy0
-#define v0 (Cy0+13)
+#define v0 (Cy0 + 13)
static const double Cj1y1[] = {
-0.4435757816794127857114720794e7, /* pr0 */
-0.9942246505077641195658377899e7,
-0.6603373248364939109255245434e7,
@@ -190,13 +190,13 @@
0.3789022974577220264142952256e5,
0.8638367769604990967475517183e3,
};
#define pr0 Cj1y1
-#define ps0 (Cj1y1+6)
-#define qr0 (Cj1y1+12)
-#define qs0 (Cj1y1+18)
+#define ps0 (Cj1y1 + 6)
+#define qr0 (Cj1y1 + 12)
+#define qs0 (Cj1y1 + 18)
static const double Cj1[] = {
-6.250000000000002203053200981413218949548e-0002, /* a0 */
1.600998455640072901321605101981501263762e-0003,
-1.963888815948313758552511884390162864930e-0005,
@@ -223,13 +223,13 @@
1.180768373106166527048240364872043816050e-0008,
5.942665743476099355323245707680648588540e-0012,
};
#define a0 Cj1
-#define b0 (Cj1+4)
-#define a1 (Cj1+8)
-#define b1 (Cj1+20)
+#define b0 (Cj1 + 4)
+#define a1 (Cj1 + 8)
+#define b1 (Cj1 + 20)
static const double Cy1[] = {
-1.960570906462389461018983259589655961560e-0001, /* c0 */
4.931824118350661953459180060007970291139e-0002,
-1.626975871565393656845930125424683008677e-0003,
@@ -257,88 +257,101 @@
1.268035774543174837829534603830227216291e-0008,
6.579416271766610825192542295821308730206e-0012,
};
#define c0 Cy1
-#define d0 (Cy1+4)
-#define c1 (Cy1+9)
-#define d1 (Cy1+21)
-
+#define d0 (Cy1 + 4)
+#define c1 (Cy1 + 9)
+#define d1 (Cy1 + 21)
/* core of j0f computation; assumes fx is finite */
static double
__k_j0f(float fx)
{
double x, z, s, c, ss, cc, r, t, p0, q0;
int ix, i;
ix = *(int *)&fx & ~0x80000000;
x = fabs((double)fx);
+
if (ix > 0x41000000) {
/* x > 8; see comments in j0.c */
s = sin(x);
c = cos(x);
+
if (signbit(s) != signbit(c)) {
ss = s - c;
cc = -cos(x + x) / ss;
} else {
cc = s + c;
ss = -cos(x + x) / cc;
}
+
if (ix > 0x501502f9) {
/* x > 1.0e10 */
p0 = one;
q0 = neighth / x;
} else {
t = eight / x;
z = t * t;
- p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
- z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
- (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
- z * (ps[4] + z * (ps[5] + z))))));
+ p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + z *
+ (pr[4] + z * (pr[5] + z * pr[6])))))) / (ps[0] + z *
+ (ps[1] + z * (ps[2] + z * (ps[3] + z * (ps[4] + z *
+ (ps[5] + z))))));
q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
- z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
- (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
- z * (qs[4] + z * (qs[5] + z))))))) * t;
+ z * (qr[4] + z * (qr[5] + z * qr[6])))))) / (qs[0] +
+ z * (qs[1] + z * (qs[2] + z * (qs[3] + z * (qs[4] +
+ z * (qs[5] + z))))))) * t;
}
+
return (isqrtpi * (p0 * cc - q0 * ss) / sqrt(x));
}
+
if (ix <= 0x3727c5ac) {
/* x <= 1.0e-5 */
if (ix <= 0x219392ef) /* x <= 1.0e-18 */
return (one - x);
+
return (one - x * x * quarter);
}
+
z = x * x;
+
if (ix <= 0x3fa3d70a) {
/* x <= 1.28 */
r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
return (one + z * (r / s));
}
+
r = r1[8];
s = s1[8];
+
for (i = 7; i >= 0; i--) {
r = r * z + r1[i];
s = s * z + s1[i];
}
+
return (r / s);
}
float
j0f(float fx)
{
float f;
int ix;
+
#if defined(__i386) && !defined(__amd64)
int rp;
#endif
ix = *(int *)&fx & ~0x80000000;
+
if (ix >= 0x7f800000) { /* nan or inf */
if (ix > 0x7f800000)
return (fx * fx);
+
return (zerof);
}
#if defined(__i386) && !defined(__amd64)
rp = __swapRP(fp_extended);
@@ -358,66 +371,79 @@
double x, z, s, c, ss, cc, t, p0, q0, u, v;
int ix, i;
ix = *(int *)&fx;
x = (double)fx;
+
if (ix > 0x41000000) {
/* x > 8; see comments in j0.c */
s = sin(x);
c = cos(x);
+
if (signbit(s) != signbit(c)) {
ss = s - c;
cc = -cos(x + x) / ss;
} else {
cc = s + c;
ss = -cos(x + x) / cc;
}
+
if (ix > 0x501502f9) {
/* x > 1.0e10 */
p0 = one;
q0 = neighth / x;
} else {
t = eight / x;
z = t * t;
- p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] +
- z * (pr[4] + z * (pr[5] + z * pr[6])))))) /
- (ps[0] + z * (ps[1] + z * (ps[2] + z * (ps[3] +
- z * (ps[4] + z * (ps[5] + z))))));
+ p0 = (pr[0] + z * (pr[1] + z * (pr[2] + z * (pr[3] + z *
+ (pr[4] + z * (pr[5] + z * pr[6])))))) / (ps[0] + z *
+ (ps[1] + z * (ps[2] + z * (ps[3] + z * (ps[4] + z *
+ (ps[5] + z))))));
q0 = ((qr[0] + z * (qr[1] + z * (qr[2] + z * (qr[3] +
- z * (qr[4] + z * (qr[5] + z * qr[6])))))) /
- (qs[0] + z * (qs[1] + z * (qs[2] + z * (qs[3] +
- z * (qs[4] + z * (qs[5] + z))))))) * t;
+ z * (qr[4] + z * (qr[5] + z * qr[6])))))) / (qs[0] +
+ z * (qs[1] + z * (qs[2] + z * (qs[3] + z * (qs[4] +
+ z * (qs[5] + z))))))) * t;
}
+
return (isqrtpi * (p0 * ss + q0 * cc) / sqrt(x));
}
+
if (ix <= 0x219392ef) /* x <= 1.0e-18 */
return (u0[0] + tpi * log(x));
+
z = x * x;
u = u0[12];
+
for (i = 11; i >= 0; i--)
u = u * z + u0[i];
+
v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
return (u / v + tpi * (__k_j0f(fx) * log(x)));
}
float
y0f(float fx)
{
float f;
int ix;
+
#if defined(__i386) && !defined(__amd64)
int rp;
#endif
ix = *(int *)&fx;
+
if ((ix & ~0x80000000) > 0x7f800000) /* nan */
return (fx * fx);
+
if (ix <= 0) { /* zero or negative */
if ((ix << 1) == 0)
return (-onef / zerof);
+
return (zerof / zerof);
}
+
if (ix == 0x7f800000) /* +inf */
return (zerof);
#if defined(__i386) && !defined(__amd64)
rp = __swapRP(fp_extended);
@@ -439,74 +465,87 @@
ix = *(int *)&fx;
sgn = (unsigned)ix >> 31;
ix &= ~0x80000000;
x = fabs((double)fx);
+
if (ix > 0x41000000) {
/* x > 8; see comments in j1.c */
s = sin(x);
c = cos(x);
+
if (signbit(s) != signbit(c)) {
cc = s - c;
ss = cos(x + x) / cc;
} else {
ss = -s - c;
cc = cos(x + x) / ss;
}
+
if (ix > 0x501502f9) {
/* x > 1.0e10 */
p1 = one;
q1 = three8 / x;
} else {
t = eight / x;
z = t * t;
- p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
- (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
- (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
- (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
+ p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * (pr0[3] +
+ z * (pr0[4] + z * pr0[5]))))) / (ps0[0] + z *
+ (ps0[1] + z * (ps0[2] + z * (ps0[3] + z * (ps0[4] +
+ z * (ps0[5] + z))))));
q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
- (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
- (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
- (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
+ (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / (qs0[0] +
+ z * (qs0[1] + z * (qs0[2] + z * (qs0[3] + z *
+ (qs0[4] + z * (qs0[5] + z))))))) * t;
}
+
t = isqrtpi * (p1 * cc - q1 * ss) / sqrt(x);
- return ((sgn)? -t : t);
+ return ((sgn) ? -t : t);
}
+
if (ix <= 0x3727c5ac) {
/* x <= 1.0e-5 */
if (ix <= 0x219392ef) /* x <= 1.0e-18 */
t = half * x;
else
t = x * (half + neighth * x * x);
- return ((sgn)? -t : t);
+
+ return ((sgn) ? -t : t);
}
+
z = x * x;
+
if (ix < 0x3fa3d70a) {
/* x < 1.28 */
r = a0[0] + z * (a0[1] + z * (a0[2] + z * a0[3]));
s = b0[0] + z * (b0[1] + z * (b0[2] + z * b0[3]));
t = x * half + x * (z * (r / s));
} else {
r = a1[11];
+
for (i = 10; i >= 0; i--)
r = r * z + a1[i];
+
s = b1[0] + z * (b1[1] + z * (b1[2] + z * (b1[3] + z * b1[4])));
t = x * (r / s);
}
- return ((sgn)? -t : t);
+
+ return ((sgn) ? -t : t);
}
float
j1f(float fx)
{
float f;
int ix;
+
#if defined(__i386) && !defined(__amd64)
int rp;
#endif
ix = *(int *)&fx & ~0x80000000;
+
if (ix >= 0x7f800000) /* nan or inf */
return (onef / fx);
#if defined(__i386) && !defined(__amd64)
rp = __swapRP(fp_extended);
@@ -526,72 +565,87 @@
double x, z, s, c, ss, cc, u, v, p1, q1, t;
int i, ix;
ix = *(int *)&fx;
x = (double)fx;
+
if (ix > 0x41000000) {
/* x > 8; see comments in j1.c */
s = sin(x);
c = cos(x);
+
if (signbit(s) != signbit(c)) {
cc = s - c;
ss = cos(x + x) / cc;
} else {
ss = -s - c;
cc = cos(x + x) / ss;
}
+
if (ix > 0x501502f9) {
/* x > 1.0e10 */
p1 = one;
q1 = three8 / x;
} else {
t = eight / x;
z = t * t;
- p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z *
- (pr0[3] + z * (pr0[4] + z * pr0[5]))))) /
- (ps0[0] + z * (ps0[1] + z * (ps0[2] + z *
- (ps0[3] + z * (ps0[4] + z * (ps0[5] + z))))));
+ p1 = (pr0[0] + z * (pr0[1] + z * (pr0[2] + z * (pr0[3] +
+ z * (pr0[4] + z * pr0[5]))))) / (ps0[0] + z *
+ (ps0[1] + z * (ps0[2] + z * (ps0[3] + z * (ps0[4] +
+ z * (ps0[5] + z))))));
q1 = ((qr0[0] + z * (qr0[1] + z * (qr0[2] + z *
- (qr0[3] + z * (qr0[4] + z * qr0[5]))))) /
- (qs0[0] + z * (qs0[1] + z * (qs0[2] + z *
- (qs0[3] + z * (qs0[4] + z * (qs0[5] + z))))))) * t;
+ (qr0[3] + z * (qr0[4] + z * qr0[5]))))) / (qs0[0] +
+ z * (qs0[1] + z * (qs0[2] + z * (qs0[3] + z *
+ (qs0[4] + z * (qs0[5] + z))))))) * t;
}
+
return (isqrtpi * (p1 * ss + q1 * cc) / sqrt(x));
}
+
if (ix <= 0x219392ef) /* x <= 1.0e-18 */
return (-tpi / x);
+
z = x * x;
+
if (ix < 0x3fa3d70a) {
/* x < 1.28 */
u = c0[0] + z * (c0[1] + z * (c0[2] + z * c0[3]));
v = d0[0] + z * (d0[1] + z * (d0[2] + z * (d0[3] + z * d0[4])));
} else {
u = c1[11];
+
for (i = 10; i >= 0; i--)
u = u * z + c1[i];
+
v = d1[0] + z * (d1[1] + z * (d1[2] + z * (d1[3] + z * d1[4])));
}
+
return (x * (u / v) + tpi * (__k_j1f(fx) * log(x) - one / x));
}
float
y1f(float fx)
{
float f;
int ix;
+
#if defined(__i386) && !defined(__amd64)
int rp;
#endif
ix = *(int *)&fx;
+
if ((ix & ~0x80000000) > 0x7f800000) /* nan */
return (fx * fx);
+
if (ix <= 0) { /* zero or negative */
if ((ix << 1) == 0)
return (-onef / zerof);
+
return (zerof / zerof);
}
+
if (ix == 0x7f800000) /* +inf */
return (zerof);
#if defined(__i386) && !defined(__amd64)
rp = __swapRP(fp_extended);
@@ -608,67 +662,81 @@
jnf(int n, float fx)
{
double a, b, temp, x, z, w, t, q0, q1, h;
float f;
int i, ix, sgn, m, k;
+
#if defined(__i386) && !defined(__amd64)
int rp;
#endif
if (n < 0) {
n = -n;
fx = -fx;
}
+
if (n == 0)
return (j0f(fx));
+
if (n == 1)
return (j1f(fx));
ix = *(int *)&fx;
- sgn = (n & 1)? ((unsigned)ix >> 31) : 0;
+ sgn = (n & 1) ? ((unsigned)ix >> 31) : 0;
ix &= ~0x80000000;
+
if (ix >= 0x7f800000) { /* nan or inf */
if (ix > 0x7f800000)
return (fx * fx);
- return ((sgn)? -zerof : zerof);
+
+ return ((sgn) ? -zerof : zerof);
}
+
if ((ix << 1) == 0)
- return ((sgn)? -zerof : zerof);
+ return ((sgn) ? -zerof : zerof);
#if defined(__i386) && !defined(__amd64)
rp = __swapRP(fp_extended);
#endif
fx = fabsf(fx);
x = (double)fx;
+
if ((double)n <= x) {
/* safe to use J(n+1,x) = 2n/x * J(n,x) - J(n-1,x) */
a = __k_j0f(fx);
b = __k_j1f(fx);
+
for (i = 1; i < n; i++) {
temp = b;
b = b * ((double)(i + i) / x) - a;
a = temp;
}
+
f = (float)b;
#if defined(__i386) && !defined(__amd64)
if (rp != fp_extended)
(void) __swapRP(rp);
#endif
- return ((sgn)? -f : f);
+ return ((sgn) ? -f : f);
}
+
if (ix < 0x3089705f) {
/* x < 1.0e-9; use J(n,x) = 1/n! * (x / 2)^n */
if (n > 6)
n = 6; /* result underflows to zero for n >= 6 */
+
b = t = half * x;
a = one;
+
for (i = 2; i <= n; i++) {
b *= t;
a *= (double)i;
}
+
b /= a;
} else {
+ /* BEGIN CSTYLED */
/*
* Use the backward recurrence:
*
* x x^2 x^2
* J(n,x)/J(n-1,x) = ---- - ------ - ------ .....
@@ -692,40 +760,47 @@
* Q(k) = (w+k*h)*Q(k-1) - Q(k-2).
*
* Then when Q(k) > 1e4, k is large enough for single
* precision.
*/
+ /* END CSTYLED */
/* XXX NOT DONE - rework this */
w = (n + n) / x;
h = two / x;
q0 = w;
z = w + h;
q1 = w * z - one;
k = 1;
+
while (q1 < big) {
k++;
z += h;
temp = z * q1 - q0;
q0 = q1;
q1 = temp;
}
+
m = n + n;
t = zero;
+
for (i = (n + k) << 1; i >= m; i -= 2)
t = one / ((double)i / x - t);
+
a = t;
b = one;
+
/*
* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
* double 7.09782712893383973096e+02
* then recurrent value may overflow and the result is
* likely underflow to zero
*/
temp = (double)n;
temp *= log((two / x) * temp);
+
if (temp < 7.09782712893383973096e+02) {
for (i = n - 1; i > 0; i--) {
temp = b;
b = b * ((double)(i + i) / x) - a;
a = temp;
@@ -733,75 +808,91 @@
} else {
for (i = n - 1; i > 0; i--) {
temp = b;
b = b * ((double)(i + i) / x) - a;
a = temp;
+
if (b > 1.0e100) {
a /= b;
t /= b;
b = one;
}
}
}
+
b = (t * __k_j0f(fx) / b);
}
+
f = (float)b;
#if defined(__i386) && !defined(__amd64)
if (rp != fp_extended)
(void) __swapRP(rp);
#endif
- return ((sgn)? -f : f);
+ return ((sgn) ? -f : f);
}
float
ynf(int n, float fx)
{
double a, b, temp, x;
float f;
int i, sign, ix;
+
#if defined(__i386) && !defined(__amd64)
int rp;
#endif
sign = 0;
+
if (n < 0) {
n = -n;
+
if (n & 1)
sign = 1;
}
+
if (n == 0)
return (y0f(fx));
+
if (n == 1)
- return ((sign)? -y1f(fx) : y1f(fx));
+ return ((sign) ? -y1f(fx) : y1f(fx));
ix = *(int *)&fx;
+
if ((ix & ~0x80000000) > 0x7f800000) /* nan */
return (fx * fx);
+
if (ix <= 0) { /* zero or negative */
if ((ix << 1) == 0)
return (-onef / zerof);
+
return (zerof / zerof);
}
+
if (ix == 0x7f800000) /* +inf */
return (zerof);
#if defined(__i386) && !defined(__amd64)
rp = __swapRP(fp_extended);
#endif
a = __k_y0f(fx);
b = __k_y1f(fx);
x = (double)fx;
+
for (i = 1; i < n; i++) {
temp = b;
b *= (double)(i + i) / x;
+
if (b <= -DBL_MAX)
break;
+
b -= a;
a = temp;
}
+
f = (float)b;
#if defined(__i386) && !defined(__amd64)
if (rp != fp_extended)
(void) __swapRP(rp);
#endif
- return ((sign)? -f : f);
+ return ((sign) ? -f : f);
}