Print this page
11210 libm should be cstyle(1ONBLD) clean
@@ -20,10 +20,11 @@
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
+
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
@@ -54,72 +55,80 @@
static long double neg(long double, int *);
static long double poly(long double, const long double *, int);
static long double polytail(long double);
static long double primary(long double);
-static const long double
-c0 = 0.0L,
-ch = 0.5L,
-c1 = 1.0L,
-c2 = 2.0L,
-c3 = 3.0L,
-c4 = 4.0L,
-c5 = 5.0L,
-c6 = 6.0L,
-pi = 3.1415926535897932384626433832795028841971L,
-tiny = 1.0e-40L;
+static const long double c0 = 0.0L,
+ ch = 0.5L,
+ c1 = 1.0L,
+ c2 = 2.0L,
+ c3 = 3.0L,
+ c4 = 4.0L,
+ c5 = 5.0L,
+ c6 = 6.0L,
+ pi = 3.1415926535897932384626433832795028841971L,
+ tiny = 1.0e-40L;
long double
-__k_lgammal(long double x, int *signgamlp) {
+__k_lgammal(long double x, int *signgamlp)
+{
long double t, y;
int i;
/* purge off +-inf, NaN and negative arguments */
if (!finitel(x))
- return (x*x);
+ return (x * x);
+
*signgamlp = 1;
+
if (signbitl(x))
return (neg(x, signgamlp));
/* for x < 8.0 */
if (x < 8.0L) {
y = anintl(x);
- i = (int) y;
+ i = (int)y;
+
switch (i) {
case 0:
+
if (x < 1.0e-40L)
return (-logl(x));
else
- return (primary(x)-log1pl(x))-logl(x);
+ return ((primary(x) - log1pl(x)) - logl(x));
+
case 1:
- return (primary(x-y)-logl(x));
+ return (primary(x - y) - logl(x));
case 2:
- return (primary(x-y));
+ return (primary(x - y));
case 3:
- return (primary(x-y)+logl(x-c1));
+ return (primary(x - y) + logl(x - c1));
case 4:
- return (primary(x-y)+logl((x-c1)*(x-c2)));
+ return (primary(x - y) + logl((x - c1) * (x - c2)));
case 5:
- return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)));
+ return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
+ c3)));
case 6:
- return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)));
+ return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
+ c3) * (x - c4)));
case 7:
- return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)));
+ return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
+ c3) * (x - c4) * (x - c5)));
case 8:
- return primary(x-y)+
- logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)*(x-c6));
+ return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
+ c3) * (x - c4) * (x - c5) * (x - c6)));
}
}
/* 8.0 <= x < 1.0e40 */
if (x < 1.0e40L) {
t = logl(x);
- return (x*(t-c1)-(ch*t-polytail(c1/x)));
+ return (x * (t - c1) - (ch * t - polytail(c1 / x)));
}
/* 1.0e40 <= x <= inf */
- return (x*(logl(x)-c1));
+ return (x * (logl(x) - c1));
}
static const long double an1[] = { /* 20 terms */
-0.0772156649015328606065120900824024309741L,
3.224670334241132182362075833230130289059e-0001L,
@@ -300,34 +309,50 @@
2.961350193868935325526962209019387821584e-0008L,
-2.834602215195368130104649234505033159842e-0009L,
};
static long double
-primary(long double s) { /* assume |s|<=0.5 */
+primary(long double s) /* assume |s|<=0.5 */
+{
int i;
- i = (int) (8.0L * (s + 0.5L));
+ i = (int)(8.0L * (s + 0.5L));
+
switch (i) {
- case 0: return ch*s+s*poly(s, an4, 21);
- case 1: return ch*s+s*poly(s, an3, 20);
- case 2: return ch*s+s*poly(s, an2, 20);
- case 3: return ch*s+s*poly(s, an1, 20);
- case 4: return ch*s+s*poly(s, ap1, 19);
- case 5: return ch*s+s*poly(s, ap2, 19);
- case 6: return ch*s+s*poly(s, ap3, 19);
- case 7: return ch*s+s*poly(s, ap4, 19);
+ case 0:
+ return (ch * s + s * poly(s, an4, 21));
+ case 1:
+ return (ch * s + s * poly(s, an3, 20));
+ case 2:
+ return (ch * s + s * poly(s, an2, 20));
+ case 3:
+ return (ch * s + s * poly(s, an1, 20));
+ case 4:
+ return (ch * s + s * poly(s, ap1, 19));
+ case 5:
+ return (ch * s + s * poly(s, ap2, 19));
+ case 6:
+ return (ch * s + s * poly(s, ap3, 19));
+ case 7:
+ return (ch * s + s * poly(s, ap4, 19));
}
+
/* NOTREACHED */
return (0.0L);
}
static long double
-poly(long double s, const long double *p, int n) {
+poly(long double s, const long double *p, int n)
+{
long double y;
int i;
- y = p[n-1];
- for (i = n-2; i >= 0; i--) y = p[i]+s*y;
+
+ y = p[n - 1];
+
+ for (i = n - 2; i >= 0; i--)
+ y = p[i] + s * y;
+
return (y);
}
static const long double pt[] = {
9.189385332046727417803297364056176804663e-0001L,
@@ -350,23 +375,30 @@
7.316526934569686459641438882340322673357e+0007L,
-3.806450279064900548836571789284896711473e+0008L,
};
static long double
-polytail(long double s) {
+polytail(long double s)
+{
long double t, z;
int i;
- z = s*s;
+
+ z = s * s;
t = pt[18];
- for (i = 17; i >= 1; i--) t = pt[i]+z*t;
- return (pt[0]+s*t);
+
+ for (i = 17; i >= 1; i--)
+ t = pt[i] + z * t;
+
+ return (pt[0] + s * t);
}
static long double
-neg(long double z, int *signgamlp) {
+neg(long double z, int *signgamlp)
+{
long double t, p;
+ /* BEGIN CSTYLED */
/*
* written by K.C. Ng, Feb 2, 1989.
*
* Since
* -z*G(-z)*G(z) = pi/sin(pi*z),
@@ -384,18 +416,24 @@
* return -log(z);
* else
* return log(pi/(t*z))-lgamma(z);
*
*/
+ /* END CSTYLED */
t = sinpil(z); /* t := sin(pi*z) */
+
if (t == c0) /* return 1.0/0.0 = +INF */
- return (c1/c0);
+ return (c1 / c0);
z = -z;
+
if (z <= tiny)
p = -logl(z);
else
- p = logl(pi/(fabsl(t)*z)) - __k_lgammal(z, signgamlp);
- if (t < c0) *signgamlp = -1;
+ p = logl(pi / (fabsl(t) * z)) - __k_lgammal(z, signgamlp);
+
+ if (t < c0)
+ *signgamlp = -1;
+
return (p);
}