Format floating routines.

1 files changed, 182 insertions, 133 deletions

diff --git a/sysdeps/ieee754/dbl-64/k_rem_pio2.c b/sysdeps/ieee754/dbl-64/k_rem_pio2.c
index ec4b4cf600..047c6c2886 100644
--- a/sysdeps/ieee754/dbl-64/k_rem_pio2.c
+++ b/sysdeps/ieee754/dbl-64/k_rem_pio2.c
@@ -147,166 +147,215 @@ static const double PIo2[] = {
 };
 
 static const double
-zero   = 0.0,
-one    = 1.0,
-two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
-twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
+  zero   = 0.0,
+  one    = 1.0,
+  two24  = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+  twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
 
-int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
+int
+__kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
+                   const int32_t *ipio2)
 {
-	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
-	double z,fw,f[20],fq[20],q[20];
+  int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
+  double z, fw, f[20], fq[20], q[20];
 
-    /* initialize jk*/
-	jk = init_jk[prec];
-	jp = jk;
+  /* initialize jk*/
+  jk = init_jk[prec];
+  jp = jk;
 
-    /* determine jx,jv,q0, note that 3>q0 */
-	jx =  nx-1;
-	jv = (e0-3)/24; if(jv<0) jv=0;
-	q0 =  e0-24*(jv+1);
+  /* determine jx,jv,q0, note that 3>q0 */
+  jx = nx - 1;
+  jv = (e0 - 3) / 24; if (jv < 0)
+    jv = 0;
+  q0 = e0 - 24 * (jv + 1);
 
-    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
-	j = jv-jx; m = jx+jk;
-	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
+  /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+  j = jv - jx; m = jx + jk;
+  for (i = 0; i <= m; i++, j++)
+    f[i] = (j < 0) ? zero : (double) ipio2[j];
 
-    /* compute q[0],q[1],...q[jk] */
-	for (i=0;i<=jk;i++) {
-	    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
-	}
+  /* compute q[0],q[1],...q[jk] */
+  for (i = 0; i <= jk; i++)
+    {
+      for (j = 0, fw = 0.0; j <= jx; j++)
+	fw += x[j] * f[jx + i - j];
+      q[i] = fw;
+    }
 
-	jz = jk;
+  jz = jk;
 recompute:
-    /* distill q[] into iq[] reversingly */
-	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
-	    fw    =  (double)((int32_t)(twon24* z));
-	    iq[i] =  (int32_t)(z-two24*fw);
-	    z     =  q[j-1]+fw;
-	}
+  /* distill q[] into iq[] reversingly */
+  for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
+    {
+      fw = (double) ((int32_t) (twon24 * z));
+      iq[i] = (int32_t) (z - two24 * fw);
+      z = q[j - 1] + fw;
+    }
 
-    /* compute n */
-	z  = __scalbn(z,q0);		/* actual value of z */
-	z -= 8.0*__floor(z*0.125);		/* trim off integer >= 8 */
-	n  = (int32_t) z;
-	z -= (double)n;
-	ih = 0;
-	if(q0>0) {	/* need iq[jz-1] to determine n */
-	    i  = (iq[jz-1]>>(24-q0)); n += i;
-	    iq[jz-1] -= i<<(24-q0);
-	    ih = iq[jz-1]>>(23-q0);
-	}
-	else if(q0==0) ih = iq[jz-1]>>23;
-	else if(z>=0.5) ih=2;
+  /* compute n */
+  z = __scalbn (z, q0);                 /* actual value of z */
+  z -= 8.0 * __floor (z * 0.125);               /* trim off integer >= 8 */
+  n = (int32_t) z;
+  z -= (double) n;
+  ih = 0;
+  if (q0 > 0)           /* need iq[jz-1] to determine n */
+    {
+      i = (iq[jz - 1] >> (24 - q0)); n += i;
+      iq[jz - 1] -= i << (24 - q0);
+      ih = iq[jz - 1] >> (23 - q0);
+    }
+  else if (q0 == 0)
+    ih = iq[jz - 1] >> 23;
+  else if (z >= 0.5)
+    ih = 2;
 
-	if(ih>0) {	/* q > 0.5 */
-	    n += 1; carry = 0;
-	    for(i=0;i<jz ;i++) {	/* compute 1-q */
-		j = iq[i];
-		if(carry==0) {
-		    if(j!=0) {
-			carry = 1; iq[i] = 0x1000000- j;
-		    }
-		} else  iq[i] = 0xffffff - j;
-	    }
-	    if(q0>0) {		/* rare case: chance is 1 in 12 */
-	        switch(q0) {
-	        case 1:
-		   iq[jz-1] &= 0x7fffff; break;
-		case 2:
-		   iq[jz-1] &= 0x3fffff; break;
-	        }
+  if (ih > 0)           /* q > 0.5 */
+    {
+      n += 1; carry = 0;
+      for (i = 0; i < jz; i++)          /* compute 1-q */
+	{
+	  j = iq[i];
+	  if (carry == 0)
+	    {
+	      if (j != 0)
+		{
+		  carry = 1; iq[i] = 0x1000000 - j;
+		}
 	    }
-	    if(ih==2) {
-		z = one - z;
-		if(carry!=0) z -= __scalbn(one,q0);
+	  else
+	    iq[i] = 0xffffff - j;
+	}
+      if (q0 > 0)               /* rare case: chance is 1 in 12 */
+	{
+	  switch (q0)
+	    {
+	    case 1:
+	      iq[jz - 1] &= 0x7fffff; break;
+	    case 2:
+	      iq[jz - 1] &= 0x3fffff; break;
 	    }
 	}
+      if (ih == 2)
+	{
+	  z = one - z;
+	  if (carry != 0)
+	    z -= __scalbn (one, q0);
+	}
+    }
 
-    /* check if recomputation is needed */
-	if(z==zero) {
-	    j = 0;
-	    for (i=jz-1;i>=jk;i--) j |= iq[i];
-	    if(j==0) { /* need recomputation */
-		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
+  /* check if recomputation is needed */
+  if (z == zero)
+    {
+      j = 0;
+      for (i = jz - 1; i >= jk; i--)
+	j |= iq[i];
+      if (j == 0)      /* need recomputation */
+	{
+	  for (k = 1; iq[jk - k] == 0; k++)
+	    ;                               /* k = no. of terms needed */
 
-		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
-		    f[jx+i] = (double) ipio2[jv+i];
-		    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
-		    q[i] = fw;
-		}
-		jz += k;
-		goto recompute;
+	  for (i = jz + 1; i <= jz + k; i++) /* add q[jz+1] to q[jz+k] */
+	    {
+	      f[jx + i] = (double) ipio2[jv + i];
+	      for (j = 0, fw = 0.0; j <= jx; j++)
+		fw += x[j] * f[jx + i - j];
+	      q[i] = fw;
 	    }
+	  jz += k;
+	  goto recompute;
 	}
+    }
 
-    /* chop off zero terms */
-	if(z==0.0) {
-	    jz -= 1; q0 -= 24;
-	    while(iq[jz]==0) { jz--; q0-=24;}
-	} else { /* break z into 24-bit if necessary */
-	    z = __scalbn(z,-q0);
-	    if(z>=two24) {
-		fw = (double)((int32_t)(twon24*z));
-		iq[jz] = (int32_t)(z-two24*fw);
-		jz += 1; q0 += 24;
-		iq[jz] = (int32_t) fw;
-	    } else iq[jz] = (int32_t) z ;
+  /* chop off zero terms */
+  if (z == 0.0)
+    {
+      jz -= 1; q0 -= 24;
+      while (iq[jz] == 0)
+	{
+	  jz--; q0 -= 24;
 	}
-
-    /* convert integer "bit" chunk to floating-point value */
-	fw = __scalbn(one,q0);
-	for(i=jz;i>=0;i--) {
-	    q[i] = fw*(double)iq[i]; fw*=twon24;
+    }
+  else           /* break z into 24-bit if necessary */
+    {
+      z = __scalbn (z, -q0);
+      if (z >= two24)
+	{
+	  fw = (double) ((int32_t) (twon24 * z));
+	  iq[jz] = (int32_t) (z - two24 * fw);
+	  jz += 1; q0 += 24;
+	  iq[jz] = (int32_t) fw;
 	}
+      else
+	iq[jz] = (int32_t) z;
+    }
 
-    /* compute PIo2[0,...,jp]*q[jz,...,0] */
-	for(i=jz;i>=0;i--) {
-	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
-	    fq[jz-i] = fw;
-	}
+  /* convert integer "bit" chunk to floating-point value */
+  fw = __scalbn (one, q0);
+  for (i = jz; i >= 0; i--)
+    {
+      q[i] = fw * (double) iq[i]; fw *= twon24;
+    }
 
-    /* compress fq[] into y[] */
-	switch(prec) {
-	    case 0:
-		fw = 0.0;
-		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw;
-		break;
-	    case 1:
-	    case 2:;
+  /* compute PIo2[0,...,jp]*q[jz,...,0] */
+  for (i = jz; i >= 0; i--)
+    {
+      for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
+	fw += PIo2[k] * q[i + k];
+      fq[jz - i] = fw;
+    }
+
+  /* compress fq[] into y[] */
+  switch (prec)
+    {
+    case 0:
+      fw = 0.0;
+      for (i = jz; i >= 0; i--)
+	fw += fq[i];
+      y[0] = (ih == 0) ? fw : -fw;
+      break;
+    case 1:
+    case 2:;
 #if __FLT_EVAL_METHOD__ != 0
-		volatile
+      volatile
 #endif
-		double fv = 0.0;
-		for (i=jz;i>=0;i--) fv += fq[i];
-		y[0] = (ih==0)? fv: -fv;
-		fv = fq[0]-fv;
-		for (i=1;i<=jz;i++) fv += fq[i];
-		y[1] = (ih==0)? fv: -fv;
-		break;
-	    case 3:	/* painful */
-		for (i=jz;i>0;i--) {
+      double fv = 0.0;
+      for (i = jz; i >= 0; i--)
+	fv += fq[i];
+      y[0] = (ih == 0) ? fv : -fv;
+      fv = fq[0] - fv;
+      for (i = 1; i <= jz; i++)
+	fv += fq[i];
+      y[1] = (ih == 0) ? fv : -fv;
+      break;
+    case 3:             /* painful */
+      for (i = jz; i > 0; i--)
+	{
 #if __FLT_EVAL_METHOD__ != 0
-		    volatile
+	  volatile
 #endif
-		    double fv = (double)(fq[i-1]+fq[i]);
-		    fq[i]  += fq[i-1]-fv;
-		    fq[i-1] = fv;
-		}
-		for (i=jz;i>1;i--) {
+	  double fv = (double) (fq[i - 1] + fq[i]);
+	  fq[i] += fq[i - 1] - fv;
+	  fq[i - 1] = fv;
+	}
+      for (i = jz; i > 1; i--)
+	{
 #if __FLT_EVAL_METHOD__ != 0
-		    volatile
+	  volatile
 #endif
-		    double fv = (double)(fq[i-1]+fq[i]);
-		    fq[i]  += fq[i-1]-fv;
-		    fq[i-1] = fv;
-		}
-		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
-		if(ih==0) {
-		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
-		} else {
-		    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
-		}
+	  double fv = (double) (fq[i - 1] + fq[i]);
+	  fq[i] += fq[i - 1] - fv;
+	  fq[i - 1] = fv;
+	}
+      for (fw = 0.0, i = jz; i >= 2; i--)
+	fw += fq[i];
+      if (ih == 0)
+	{
+	  y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
+	}
+      else
+	{
+	  y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
 	}
-	return n&7;
+    }
+  return n & 7;
 }