2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
};
-static const double
-zero = 0.0,
-one = 1.0,
-two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
-twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
-
int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
{
int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
/* 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];
+ f[i] = j<0 ? 0.0 : (double)ipio2[j];
/* compute q[0],q[1],...q[jk] */
for (i=0; i<=jk; i++) {
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);
+ fw = (double)(int32_t)(0x1p-24*z);
+ iq[i] = (int32_t)(z - 0x1p24*fw);
z = q[j-1]+fw;
}
if (carry == 0) {
if (j != 0) {
carry = 1;
- iq[i] = 0x1000000- j;
+ iq[i] = 0x1000000 - j;
}
} else
iq[i] = 0xffffff - j;
}
}
if (ih == 2) {
- z = one - z;
+ z = 1.0 - z;
if (carry != 0)
- z -= scalbn(one,q0);
+ z -= scalbn(1.0,q0);
}
}
/* check if recomputation is needed */
- if (z == zero) {
+ if (z == 0.0) {
j = 0;
for (i=jz-1; i>=jk; i--) j |= iq[i];
if (j == 0) { /* need recomputation */
}
} 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);
+ if (z >= 0x1p24) {
+ fw = (double)(int32_t)(0x1p-24*z);
+ iq[jz] = (int32_t)(z - 0x1p24*fw);
jz += 1;
q0 += 24;
iq[jz] = (int32_t)fw;
}
/* convert integer "bit" chunk to floating-point value */
- fw = scalbn(one,q0);
+ fw = scalbn(1.0,q0);
for (i=jz; i>=0; i--) {
q[i] = fw*(double)iq[i];
- fw *= twon24;
+ fw *= 0x1p-24;
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */
fw = 0.0;
for (i=jz; i>=0; i--)
fw += fq[i];
- STRICT_ASSIGN(double,fw,fw);
+ // TODO: drop excess precision here once double_t is used
+ fw = (double)fw;
y[0] = ih==0 ? fw : -fw;
fw = fq[0]-fw;
for (i=1; i<=jz; i++)