about summary refs log tree commit diff
path: root/sysdeps/ieee754/dbl-64/s_log1p.c
diff options
context:
space:
mode:
authorOndřej Bílka <neleai@seznam.cz>2013-10-17 16:03:24 +0200
committerOndřej Bílka <neleai@seznam.cz>2013-10-17 16:03:24 +0200
commitc5d5d574cbfa96d0f6c1db24d1e072c472627e41 (patch)
tree83b97e29ee65636dfe1247ea8d2344ca3f0b04b4 /sysdeps/ieee754/dbl-64/s_log1p.c
parente5c2c2d0c0315ca24cc9cd638cdb1a2d8dcc4b0d (diff)
downloadglibc-c5d5d574cbfa96d0f6c1db24d1e072c472627e41.tar.gz
glibc-c5d5d574cbfa96d0f6c1db24d1e072c472627e41.tar.xz
glibc-c5d5d574cbfa96d0f6c1db24d1e072c472627e41.zip
Format floating routines.
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_log1p.c')
-rw-r--r--sysdeps/ieee754/dbl-64/s_log1p.c167
1 files changed, 96 insertions, 71 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c
index 0c6f76371a..ea1dc6ca76 100644
--- a/sysdeps/ieee754/dbl-64/s_log1p.c
+++ b/sysdeps/ieee754/dbl-64/s_log1p.c
@@ -11,7 +11,7 @@
  */
 /* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
    for performance improvement on pipelined processors.
-*/
+ */
 
 /* double log1p(double x)
  *
@@ -82,87 +82,112 @@
 #include <math_private.h>
 
 static const double
-ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
-ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
-two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
-Lp[] = {0.0, 6.666666666666735130e-01,  /* 3FE55555 55555593 */
- 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
- 2.857142874366239149e-01,  /* 3FD24924 94229359 */
- 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
- 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
- 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
- 1.479819860511658591e-01};  /* 3FC2F112 DF3E5244 */
+  ln2_hi = 6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */
+  ln2_lo = 1.90821492927058770002e-10,  /* 3dea39ef 35793c76 */
+  two54 = 1.80143985094819840000e+16,   /* 43500000 00000000 */
+  Lp[] = { 0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */
+	   3.999999999940941908e-01, /* 3FD99999 9997FA04 */
+	   2.857142874366239149e-01, /* 3FD24924 94229359 */
+	   2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
+	   1.818357216161805012e-01, /* 3FC74664 96CB03DE */
+	   1.531383769920937332e-01, /* 3FC39A09 D078C69F */
+	   1.479819860511658591e-01 }; /* 3FC2F112 DF3E5244 */
 
 static const double zero = 0.0;
 
 double
-__log1p(double x)
+__log1p (double x)
 {
-	double hfsq,f,c,s,z,R,u,z2,z4,z6,R1,R2,R3,R4;
-	int32_t k,hx,hu,ax;
+  double hfsq, f, c, s, z, R, u, z2, z4, z6, R1, R2, R3, R4;
+  int32_t k, hx, hu, ax;
 
-	GET_HIGH_WORD(hx,x);
-	ax = hx&0x7fffffff;
+  GET_HIGH_WORD (hx, x);
+  ax = hx & 0x7fffffff;
 
-	k = 1;
-	if (hx < 0x3FDA827A) {			/* x < 0.41422  */
-	    if(__builtin_expect(ax>=0x3ff00000, 0)) { /* x <= -1.0 */
-		if(x==-1.0) return -two54/(x-x);/* log1p(-1)=+inf */
-		else return (x-x)/(x-x);	/* log1p(x<-1)=NaN */
-	    }
-	    if(__builtin_expect(ax<0x3e200000, 0)) { /* |x| < 2**-29 */
-		math_force_eval(two54+x);	/* raise inexact */
-		if (ax<0x3c900000)		/* |x| < 2**-54 */
-		    return x;
-		else
-		    return x - x*x*0.5;
-	    }
-	    if(hx>0||hx<=((int32_t)0xbfd2bec3)) {
-		k=0;f=x;hu=1;}	/* -0.2929<x<0.41422 */
+  k = 1;
+  if (hx < 0x3FDA827A)                          /* x < 0.41422  */
+    {
+      if (__builtin_expect (ax >= 0x3ff00000, 0))     /* x <= -1.0 */
+	{
+	  if (x == -1.0)
+	    return -two54 / (x - x);            /* log1p(-1)=+inf */
+	  else
+	    return (x - x) / (x - x);           /* log1p(x<-1)=NaN */
 	}
-	else if (__builtin_expect(hx >= 0x7ff00000, 0)) return x+x;
-	if(k!=0) {
-	    if(hx<0x43400000) {
-		u  = 1.0+x;
-		GET_HIGH_WORD(hu,u);
-		k  = (hu>>20)-1023;
-		c  = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
-		c /= u;
-	    } else {
-		u  = x;
-		GET_HIGH_WORD(hu,u);
-		k  = (hu>>20)-1023;
-		c  = 0;
-	    }
-	    hu &= 0x000fffff;
-	    if(hu<0x6a09e) {
-		SET_HIGH_WORD(u,hu|0x3ff00000);	/* normalize u */
-	    } else {
-		k += 1;
-		SET_HIGH_WORD(u,hu|0x3fe00000);	/* normalize u/2 */
-		hu = (0x00100000-hu)>>2;
-	    }
-	    f = u-1.0;
+      if (__builtin_expect (ax < 0x3e200000, 0))     /* |x| < 2**-29 */
+	{
+	  math_force_eval (two54 + x);          /* raise inexact */
+	  if (ax < 0x3c900000)                  /* |x| < 2**-54 */
+	    return x;
+	  else
+	    return x - x * x * 0.5;
+	}
+      if (hx > 0 || hx <= ((int32_t) 0xbfd2bec3))
+	{
+	  k = 0; f = x; hu = 1;
+	}                       /* -0.2929<x<0.41422 */
+    }
+  else if (__builtin_expect (hx >= 0x7ff00000, 0))
+    return x + x;
+  if (k != 0)
+    {
+      if (hx < 0x43400000)
+	{
+	  u = 1.0 + x;
+	  GET_HIGH_WORD (hu, u);
+	  k = (hu >> 20) - 1023;
+	  c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0); /* correction term */
+	  c /= u;
+	}
+      else
+	{
+	  u = x;
+	  GET_HIGH_WORD (hu, u);
+	  k = (hu >> 20) - 1023;
+	  c = 0;
+	}
+      hu &= 0x000fffff;
+      if (hu < 0x6a09e)
+	{
+	  SET_HIGH_WORD (u, hu | 0x3ff00000);   /* normalize u */
+	}
+      else
+	{
+	  k += 1;
+	  SET_HIGH_WORD (u, hu | 0x3fe00000);   /* normalize u/2 */
+	  hu = (0x00100000 - hu) >> 2;
 	}
-	hfsq=0.5*f*f;
-	if(hu==0) {	/* |f| < 2**-20 */
-	    if(f==zero) {
-	      if(k==0) return zero;
-			else {c += k*ln2_lo; return k*ln2_hi+c;}
+      f = u - 1.0;
+    }
+  hfsq = 0.5 * f * f;
+  if (hu == 0)          /* |f| < 2**-20 */
+    {
+      if (f == zero)
+	{
+	  if (k == 0)
+	    return zero;
+	  else
+	    {
+	      c += k * ln2_lo; return k * ln2_hi + c;
 	    }
-	    R = hfsq*(1.0-0.66666666666666666*f);
-	    if(k==0) return f-R; else
-		     return k*ln2_hi-((R-(k*ln2_lo+c))-f);
 	}
-	s = f/(2.0+f);
-	z = s*s;
-	R1 = z*Lp[1]; z2=z*z;
-	R2 = Lp[2]+z*Lp[3]; z4=z2*z2;
-	R3 = Lp[4]+z*Lp[5]; z6=z4*z2;
-	R4 = Lp[6]+z*Lp[7];
-	R = R1 + z2*R2 + z4*R3 + z6*R4;
-	if(k==0) return f-(hfsq-s*(hfsq+R)); else
-		 return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
+      R = hfsq * (1.0 - 0.66666666666666666 * f);
+      if (k == 0)
+	return f - R;
+      else
+	return k * ln2_hi - ((R - (k * ln2_lo + c)) - f);
+    }
+  s = f / (2.0 + f);
+  z = s * s;
+  R1 = z * Lp[1]; z2 = z * z;
+  R2 = Lp[2] + z * Lp[3]; z4 = z2 * z2;
+  R3 = Lp[4] + z * Lp[5]; z6 = z4 * z2;
+  R4 = Lp[6] + z * Lp[7];
+  R = R1 + z2 * R2 + z4 * R3 + z6 * R4;
+  if (k == 0)
+    return f - (hfsq - s * (hfsq + R));
+  else
+    return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f);
 }
 weak_alias (__log1p, log1p)
 #ifdef NO_LONG_DOUBLE