Print this page
11210 libm should be cstyle(1ONBLD) clean
*** 20,29 ****
--- 20,30 ----
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
*** 43,74 ****
#include "libm_protos.h"
#include <math.h>
#include <values.h>
#define GENERIC double
- static const GENERIC
- zero = 0.0,
- small = 1.0e-5,
- tiny = 1.0e-18,
- one = 1.0,
- eight = 8.0,
- invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
- tpi = 0.636619772367581343075535053490057448;
! static GENERIC pzero(GENERIC), qzero(GENERIC);
static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
-2.500000000000003622131880894830476755537e-0001,
1.095597547334830263234433855932375353303e-0002,
-1.819734750463320921799187258987098087697e-0004,
9.977001946806131657544212501069893930846e-0007,
};
static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
1.0,
1.867609810662950169966782360588199673741e-0002,
1.590389206181565490878430827706972074208e-0004,
6.520867386742583632375520147714499522721e-0007,
};
static const GENERIC r1[9] = { /* [1.28,8] */
9.999999999999999942156495584397047660949e-0001,
-2.389887722731319130476839836908143731281e-0001,
1.293359476138939027791270393439493640570e-0002,
-2.770985642343140122168852400228563364082e-0004,
--- 44,79 ----
#include "libm_protos.h"
#include <math.h>
#include <values.h>
#define GENERIC double
! static const GENERIC zero = 0.0,
! small = 1.0e-5,
! tiny = 1.0e-18,
! one = 1.0,
! eight = 8.0,
! invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
! tpi = 0.636619772367581343075535053490057448;
!
! static GENERIC pzero(GENERIC);
! static GENERIC qzero(GENERIC);
!
static const GENERIC r0[4] = { /* [1.e-5, 1.28] */
-2.500000000000003622131880894830476755537e-0001,
1.095597547334830263234433855932375353303e-0002,
-1.819734750463320921799187258987098087697e-0004,
9.977001946806131657544212501069893930846e-0007,
};
+
static const GENERIC s0[4] = { /* [1.e-5, 1.28] */
1.0,
1.867609810662950169966782360588199673741e-0002,
1.590389206181565490878430827706972074208e-0004,
6.520867386742583632375520147714499522721e-0007,
};
+
static const GENERIC r1[9] = { /* [1.28,8] */
9.999999999999999942156495584397047660949e-0001,
-2.389887722731319130476839836908143731281e-0001,
1.293359476138939027791270393439493640570e-0002,
-2.770985642343140122168852400228563364082e-0004,
*** 76,85 ****
--- 81,91 ----
-1.636846356264052597969042009265043251279e-0008,
5.072306160724884775085431059052611737827e-0011,
-8.187060730684066824228914775146536139112e-0014,
5.422219326959949863954297860723723423842e-0017,
};
+
static const GENERIC s1[9] = { /* [1.28,8] */
1.0,
1.101122772686807702762104741932076228349e-0002,
6.140169310641649223411427764669143978228e-0005,
2.292035877515152097976946119293215705250e-0007,
*** 89,111 ****
2.801559820648939665270492520004836611187e-0018,
2.073101088320349159764410261466350732968e-0021,
};
GENERIC
! j0(GENERIC x) {
GENERIC z, s, c, ss, cc, r, u, v, ox;
int i;
if (isnan(x))
! return (x*x); /* + -> * for Cheetah */
ox = x;
x = fabs(x);
if (x > 8.0) {
if (!finite(x))
return (zero);
s = sin(x);
c = cos(x);
/*
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
--- 95,123 ----
2.801559820648939665270492520004836611187e-0018,
2.073101088320349159764410261466350732968e-0021,
};
GENERIC
! j0(GENERIC x)
! {
GENERIC z, s, c, ss, cc, r, u, v, ox;
int i;
if (isnan(x))
! return (x * x); /* + -> * for Cheetah */
!
ox = x;
x = fabs(x);
+
if (x > 8.0) {
if (!finite(x))
return (zero);
+
s = sin(x);
c = cos(x);
+
+ /* BEGIN CSTYLED */
/*
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
*** 114,165 ****
* = 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.
*/
if (x > 8.9e307) { /* x+x may overflow */
! ss = s-c;
! cc = s+c;
} else if (signbit(s) != signbit(c)) {
ss = s - c;
! cc = -cos(x+x)/ss;
} else {
cc = s + c;
! ss = -cos(x+x)/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 (x > 1.0e40) z = (invsqrtpi*cc)/sqrt(x);
! else {
! u = pzero(x); v = qzero(x);
! z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
}
/* force to pass SVR4 even the result is wrong (sign) */
if (x > X_TLOSS)
return (_SVID_libm_err(ox, z, 34));
else
return (z);
}
if (x <= small) {
if (x <= tiny)
! return (one-x);
else
! return (one-x*x*0.25);
}
! z = x*x;
if (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));
} else {
for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
! r = r*z + r1[i];
! s = s*z + s1[i];
}
! return (r/s);
}
}
static const GENERIC u0[13] = {
-7.380429510868722526754723020704317641941e-0002,
--- 126,186 ----
* = 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.
*/
+ /* END CSTYLED */
if (x > 8.9e307) { /* x+x may overflow */
! ss = s - c;
! cc = s + c;
} else if (signbit(s) != signbit(c)) {
ss = s - c;
! cc = -cos(x + x) / ss;
} else {
cc = s + c;
! ss = -cos(x + x) / 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 (x > 1.0e40) {
! z = (invsqrtpi * cc) / sqrt(x);
! } else {
! u = pzero(x);
! v = qzero(x);
! z = invsqrtpi * (u * cc - v * ss) / sqrt(x);
}
+
/* force to pass SVR4 even the result is wrong (sign) */
if (x > X_TLOSS)
return (_SVID_libm_err(ox, z, 34));
else
return (z);
}
+
if (x <= small) {
if (x <= tiny)
! return (one - x);
else
! return (one - x * x * 0.25);
}
!
! z = x * x;
!
if (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));
} else {
for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
! r = r * z + r1[i];
! s = s * z + s1[i];
}
!
! return (r / s);
}
}
static const GENERIC u0[13] = {
-7.380429510868722526754723020704317641941e-0002,
*** 174,211 ****
7.182052899726138381739945881914874579696e-0018,
-9.060556574619677567323741194079797987200e-0021,
7.124435467408860515265552217131230511455e-0024,
-2.709726774636397615328813121715432044771e-0027,
};
static const GENERIC v0[5] = {
1.0,
4.678678931512549002587702477349214886475e-0003,
9.486828955529948534822800829497565178985e-0006,
1.001495929158861646659010844136682454906e-0008,
4.725338116256021660204443235685358593611e-0012,
};
GENERIC
! y0(GENERIC x) {
GENERIC z, /* d, */ s, c, ss, cc, u, v;
int i;
if (isnan(x))
! return (x*x); /* + -> * for Cheetah */
if (x <= zero) {
if (x == zero)
/* d= -one/(x-x); */
return (_SVID_libm_err(x, x, 8));
else
/* d = zero/(x-x); */
return (_SVID_libm_err(x, x, 9));
}
if (x > 8.0) {
if (!finite(x))
return (zero);
s = sin(x);
c = cos(x);
/*
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
--- 195,239 ----
7.182052899726138381739945881914874579696e-0018,
-9.060556574619677567323741194079797987200e-0021,
7.124435467408860515265552217131230511455e-0024,
-2.709726774636397615328813121715432044771e-0027,
};
+
static const GENERIC v0[5] = {
1.0,
4.678678931512549002587702477349214886475e-0003,
9.486828955529948534822800829497565178985e-0006,
1.001495929158861646659010844136682454906e-0008,
4.725338116256021660204443235685358593611e-0012,
};
GENERIC
! y0(GENERIC x)
! {
GENERIC z, /* d, */ s, c, ss, cc, u, v;
int i;
if (isnan(x))
! return (x * x); /* + -> * for Cheetah */
!
if (x <= zero) {
if (x == zero)
/* d= -one/(x-x); */
return (_SVID_libm_err(x, x, 8));
else
/* d = zero/(x-x); */
return (_SVID_libm_err(x, x, 9));
}
+
if (x > 8.0) {
if (!finite(x))
return (zero);
+
s = sin(x);
c = cos(x);
+
+ /* BEGIN CSTYLED */
/*
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4
* Better formula:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
*** 214,321 ****
* = 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.
*/
if (x > 8.9e307) { /* x+x may overflow */
! ss = s-c;
! cc = s+c;
} else if (signbit(s) != signbit(c)) {
ss = s - c;
! cc = -cos(x+x)/ss;
} else {
cc = s + c;
! ss = -cos(x+x)/cc;
}
/*
* j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
* y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
*/
if (x > 1.0e40)
! z = (invsqrtpi*ss)/sqrt(x);
else
! z = invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x);
if (x > X_TLOSS)
return (_SVID_libm_err(x, z, 35));
else
return (z);
-
- }
- if (x <= tiny) {
- return (u0[0] + tpi*log(x));
}
! z = x*x;
! for (u = u0[12], 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*(j0(x)*log(x)));
}
static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
! .4861344183386052721391238447e5,
! .1377662549407112278133438945e6,
! .1222466364088289731869114004e6,
! .4107070084315176135583353374e5,
! .5026073801860637125889039915e4,
! .1783193659125479654541542419e3,
.88010344055383421691677564e0,
};
static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
! .4861344183386052721414037058e5,
! .1378196632630384670477582699e6,
! .1223967185341006542748936787e6,
! .4120150243795353639995862617e5,
! .5068271181053546392490184353e4,
! .1829817905472769960535671664e3,
1.0,
};
- static const GENERIC huge = 1.0e10;
static GENERIC
! pzero(GENERIC x) {
GENERIC s, r, t, z;
int i;
if (x > huge)
return (one);
! t = eight/x; z = t*t;
! r = pr[5]+z*pr[6];
! s = ps[5]+z;
for (i = 4; i >= 0; i--) {
! r = r*z + pr[i];
! s = s*z + ps[i];
}
! return (r/s);
}
static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
! -.1731210995701068539185611951e3,
! -.5522559165936166961235240613e3,
! -.5604935606637346590614529613e3,
! -.2200430300226009379477365011e3,
! -.323869355375648849771296746e2,
! -.14294979207907956223499258e1,
-.834690374102384988158918e-2,
};
static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
! .1107975037248683865326709645e5,
! .3544581680627082674651471873e5,
! .3619118937918394132179019059e5,
! .1439895563565398007471485822e5,
! .2190277023344363955930226234e4,
! .106695157020407986137501682e3,
1.0,
};
static GENERIC
! qzero(GENERIC x) {
GENERIC s, r, t, z;
int i;
if (x > huge)
! return (-0.125/x);
! t = eight/x; z = t*t;
! r = qr[5]+z*qr[6];
! s = qs[5]+z;
for (i = 4; i >= 0; i--) {
! r = r*z + qr[i];
! s = s*z + qs[i];
}
! return (t*(r/s));
}
--- 242,358 ----
* = 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.
*/
+ /* END CSTYLED */
if (x > 8.9e307) { /* x+x may overflow */
! ss = s - c;
! cc = s + c;
} else if (signbit(s) != signbit(c)) {
ss = s - c;
! cc = -cos(x + x) / ss;
} else {
cc = s + c;
! ss = -cos(x + x) / cc;
}
+
/*
* j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
* y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
*/
if (x > 1.0e40)
! z = (invsqrtpi * ss) / sqrt(x);
else
! z = invsqrtpi * (pzero(x) * ss + qzero(x) * cc) /
! sqrt(x);
!
if (x > X_TLOSS)
return (_SVID_libm_err(x, z, 35));
else
return (z);
}
!
! if (x <= tiny)
! return (u0[0] + tpi * log(x));
!
! z = x * x;
!
! for (u = u0[12], 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 * (j0(x) * log(x)));
}
static const GENERIC pr[7] = { /* [8 -- inf] pzero 6550 */
! .4861344183386052721391238447e5, .1377662549407112278133438945e6,
! .1222466364088289731869114004e6, .4107070084315176135583353374e5,
! .5026073801860637125889039915e4, .1783193659125479654541542419e3,
.88010344055383421691677564e0,
};
+
static const GENERIC ps[7] = { /* [8 -- inf] pzero 6550 */
! .4861344183386052721414037058e5, .1378196632630384670477582699e6,
! .1223967185341006542748936787e6, .4120150243795353639995862617e5,
! .5068271181053546392490184353e4, .1829817905472769960535671664e3,
1.0,
};
+ static const GENERIC huge = 1.0e10;
static GENERIC
! pzero(GENERIC x)
! {
GENERIC s, r, t, z;
int i;
+
if (x > huge)
return (one);
!
! t = eight / x;
! z = t * t;
! r = pr[5] + z * pr[6];
! s = ps[5] + z;
!
for (i = 4; i >= 0; i--) {
! r = r * z + pr[i];
! s = s * z + ps[i];
}
!
! return (r / s);
}
static const GENERIC qr[7] = { /* [8 -- inf] qzero 6950 */
! -.1731210995701068539185611951e3, -.5522559165936166961235240613e3,
! -.5604935606637346590614529613e3, -.2200430300226009379477365011e3,
! -.323869355375648849771296746e2, -.14294979207907956223499258e1,
-.834690374102384988158918e-2,
};
+
static const GENERIC qs[7] = { /* [8 -- inf] qzero 6950 */
! .1107975037248683865326709645e5, .3544581680627082674651471873e5,
! .3619118937918394132179019059e5, .1439895563565398007471485822e5,
! .2190277023344363955930226234e4, .106695157020407986137501682e3,
1.0,
};
static GENERIC
! qzero(GENERIC x)
! {
GENERIC s, r, t, z;
int i;
+
if (x > huge)
! return (-0.125 / x);
!
! t = eight / x;
! z = t * t;
! r = qr[5] + z * qr[6];
! s = qs[5] + z;
!
for (i = 4; i >= 0; i--) {
! r = r * z + qr[i];
! s = s * z + qs[i];
}
!
! return (t * (r / s));
}