diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/k_rem_pio2.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/k_rem_pio2.c | 315 |
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; } |