summaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/k_rem_pio2.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/k_rem_pio2.c')
-rw-r--r--sysdeps/ieee754/dbl-64/k_rem_pio2.c315
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;
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-26 00:16:55 +0000