mirror of
https://github.com/wine-mirror/wine
synced 2024-07-09 04:16:08 +00:00
msvcrt: Use the __rem_pio2()/__rem_pio2f() implementation from the bundled musl library.
This commit is contained in:
Alexandre Julliard
parent
8252be560e
commit
b9b19ea2db
|
|
@ -218,212 +218,6 @@ float CDECL _chgsignf( float num )
|
|||
|
||||
#endif
|
||||
|
||||
/* Copied from musl: src/math/__rem_pio2_large.c */
|
||||
static int __rem_pio2_large(double *x, double *y, int e0, int nx, int prec)
|
||||
{
|
||||
static const int init_jk[] = {3, 4};
|
||||
static const INT32 ipio2[] = {
|
||||
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
|
||||
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
|
||||
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
|
||||
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
|
||||
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
|
||||
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
|
||||
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
|
||||
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
|
||||
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
|
||||
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
|
||||
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
|
||||
};
|
||||
static const double PIo2[] = {
|
||||
1.57079625129699707031e+00,
|
||||
7.54978941586159635335e-08,
|
||||
5.39030252995776476554e-15,
|
||||
3.28200341580791294123e-22,
|
||||
1.27065575308067607349e-29,
|
||||
1.22933308981111328932e-36,
|
||||
2.73370053816464559624e-44,
|
||||
2.16741683877804819444e-51,
|
||||
};
|
||||
|
||||
INT32 jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
|
||||
double z, fw, f[20], fq[20] = {0}, q[20];
|
||||
|
||||
/* 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);
|
||||
|
||||
/* 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 ? 0.0 : (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;
|
||||
}
|
||||
|
||||
jz = jk;
|
||||
recompute:
|
||||
/* distill q[] into iq[] reversingly */
|
||||
for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
|
||||
fw = (double)(INT32)(0x1p-24 * z);
|
||||
iq[i] = (INT32)(z - 0x1p24 * 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)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 == 2) {
|
||||
z = 1.0 - z;
|
||||
if (carry != 0)
|
||||
z -= scalbn(1.0, q0);
|
||||
}
|
||||
}
|
||||
|
||||
/* check if recomputation is needed */
|
||||
if (z == 0.0) {
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
/* 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 >= 0x1p24) {
|
||||
fw = (double)(INT32)(0x1p-24 * z);
|
||||
iq[jz] = (INT32)(z - 0x1p24 * fw);
|
||||
jz += 1;
|
||||
q0 += 24;
|
||||
iq[jz] = (INT32)fw;
|
||||
} else
|
||||
iq[jz] = (INT32)z;
|
||||
}
|
||||
|
||||
/* convert integer "bit" chunk to floating-point value */
|
||||
fw = scalbn(1.0, q0);
|
||||
for (i = jz; i >= 0; i--) {
|
||||
q[i] = fw * (double)iq[i];
|
||||
fw *= 0x1p-24;
|
||||
}
|
||||
|
||||
/* 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:
|
||||
fw = 0.0;
|
||||
for (i = jz; i >= 0; i--)
|
||||
fw += fq[i];
|
||||
fw = (double)fw;
|
||||
y[0] = ih==0 ? fw : -fw;
|
||||
fw = fq[0] - fw;
|
||||
for (i = 1; i <= jz; i++)
|
||||
fw += fq[i];
|
||||
y[1] = ih == 0 ? fw : -fw;
|
||||
break;
|
||||
case 3: /* painful */
|
||||
for (i = jz; i > 0; i--) {
|
||||
fw = fq[i - 1] + fq[i];
|
||||
fq[i] += fq[i - 1] - fw;
|
||||
fq[i - 1] = fw;
|
||||
}
|
||||
for (i = jz; i > 1; i--) {
|
||||
fw = fq[i - 1] + fq[i];
|
||||
fq[i] += fq[i - 1] - fw;
|
||||
fq[i - 1] = fw;
|
||||
}
|
||||
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;
|
||||
}
|
||||
|
||||
#ifndef __i386__
|
||||
/* Copied from musl: src/math/__sindf.c */
|
||||
static float __sindf(double x)
|
||||
|
|
@ -731,56 +525,7 @@ float CDECL atanf( float x )
|
|||
return sign ? -z : z;
|
||||
}
|
||||
|
||||
/* Copied from musl: src/math/__rem_pio2f.c */
|
||||
static int __rem_pio2f(float x, double *y)
|
||||
{
|
||||
static const double toint = 1.5 / DBL_EPSILON,
|
||||
pio4 = 0x1.921fb6p-1,
|
||||
invpio2 = 6.36619772367581382433e-01,
|
||||
pio2_1 = 1.57079631090164184570e+00,
|
||||
pio2_1t = 1.58932547735281966916e-08;
|
||||
|
||||
union {float f; uint32_t i;} u = {x};
|
||||
double tx[1], ty[1], fn;
|
||||
UINT32 ix;
|
||||
int n, sign, e0;
|
||||
|
||||
ix = u.i & 0x7fffffff;
|
||||
/* 25+53 bit pi is good enough for medium size */
|
||||
if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */
|
||||
/* Use a specialized rint() to get fn. */
|
||||
fn = fp_barrier(x * invpio2 + toint) - toint;
|
||||
n = (int)fn;
|
||||
*y = x - fn * pio2_1 - fn * pio2_1t;
|
||||
/* Matters with directed rounding. */
|
||||
if (*y < -pio4) {
|
||||
n--;
|
||||
fn--;
|
||||
*y = x - fn * pio2_1 - fn * pio2_1t;
|
||||
} else if (*y > pio4) {
|
||||
n++;
|
||||
fn++;
|
||||
*y = x - fn * pio2_1 - fn * pio2_1t;
|
||||
}
|
||||
return n;
|
||||
}
|
||||
if(ix >= 0x7f800000) { /* x is inf or NaN */
|
||||
*y = x - x;
|
||||
return 0;
|
||||
}
|
||||
/* scale x into [2^23, 2^24-1] */
|
||||
sign = u.i >> 31;
|
||||
e0 = (ix >> 23) - (0x7f + 23); /* e0 = ilogb(|x|)-23, positive */
|
||||
u.i = ix - (e0 << 23);
|
||||
tx[0] = u.f;
|
||||
n = __rem_pio2_large(tx, ty, e0, 1, 0);
|
||||
if (sign) {
|
||||
*y = -ty[0];
|
||||
return -n;
|
||||
}
|
||||
*y = ty[0];
|
||||
return n;
|
||||
}
|
||||
extern int __rem_pio2f(float x, double *y);
|
||||
|
||||
/*********************************************************************
|
||||
* cosf (MSVCRT.@)
|
||||
|
|
@ -1627,8 +1372,12 @@ double CDECL atan( double x )
|
|||
return sign ? -z : z;
|
||||
}
|
||||
|
||||
/* Copied from musl: src/math/rint.c */
|
||||
static double __rint(double x)
|
||||
/*********************************************************************
|
||||
* rint (MSVCR120.@)
|
||||
*
|
||||
* Copied from musl: src/math/rint.c
|
||||
*/
|
||||
double CDECL rint(double x)
|
||||
{
|
||||
static const double toint = 1 / DBL_EPSILON;
|
||||
|
||||
|
|
@ -1654,156 +1403,7 @@ static double __rint(double x)
|
|||
return y;
|
||||
}
|
||||
|
||||
/* Copied from musl: src/math/__rem_pio2.c */
|
||||
static int __rem_pio2(double x, double *y)
|
||||
{
|
||||
static const double pio4 = 0x1.921fb54442d18p-1,
|
||||
invpio2 = 6.36619772367581382433e-01,
|
||||
pio2_1 = 1.57079632673412561417e+00,
|
||||
pio2_1t = 6.07710050650619224932e-11,
|
||||
pio2_2 = 6.07710050630396597660e-11,
|
||||
pio2_2t = 2.02226624879595063154e-21,
|
||||
pio2_3 = 2.02226624871116645580e-21,
|
||||
pio2_3t = 8.47842766036889956997e-32;
|
||||
|
||||
union {double f; UINT64 i;} u = {x};
|
||||
double z, w, t, r, fn, tx[3], ty[2];
|
||||
UINT32 ix;
|
||||
int sign, n, ex, ey, i;
|
||||
|
||||
sign = u.i >> 63;
|
||||
ix = u.i >> 32 & 0x7fffffff;
|
||||
if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
|
||||
if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
|
||||
goto medium; /* cancellation -- use medium case */
|
||||
if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
|
||||
if (!sign) {
|
||||
z = x - pio2_1; /* one round good to 85 bits */
|
||||
y[0] = z - pio2_1t;
|
||||
y[1] = (z - y[0]) - pio2_1t;
|
||||
return 1;
|
||||
} else {
|
||||
z = x + pio2_1;
|
||||
y[0] = z + pio2_1t;
|
||||
y[1] = (z - y[0]) + pio2_1t;
|
||||
return -1;
|
||||
}
|
||||
} else {
|
||||
if (!sign) {
|
||||
z = x - 2 * pio2_1;
|
||||
y[0] = z - 2 * pio2_1t;
|
||||
y[1] = (z - y[0]) - 2 * pio2_1t;
|
||||
return 2;
|
||||
} else {
|
||||
z = x + 2 * pio2_1;
|
||||
y[0] = z + 2 * pio2_1t;
|
||||
y[1] = (z - y[0]) + 2 * pio2_1t;
|
||||
return -2;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */
|
||||
if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
|
||||
if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */
|
||||
goto medium;
|
||||
if (!sign) {
|
||||
z = x - 3 * pio2_1;
|
||||
y[0] = z - 3 * pio2_1t;
|
||||
y[1] = (z - y[0]) - 3 * pio2_1t;
|
||||
return 3;
|
||||
} else {
|
||||
z = x + 3 * pio2_1;
|
||||
y[0] = z + 3 * pio2_1t;
|
||||
y[1] = (z - y[0]) + 3 * pio2_1t;
|
||||
return -3;
|
||||
}
|
||||
} else {
|
||||
if (ix == 0x401921fb) /* |x| ~= 4pi/2 */
|
||||
goto medium;
|
||||
if (!sign) {
|
||||
z = x - 4 * pio2_1;
|
||||
y[0] = z - 4 * pio2_1t;
|
||||
y[1] = (z - y[0]) - 4 * pio2_1t;
|
||||
return 4;
|
||||
} else {
|
||||
z = x + 4 * pio2_1;
|
||||
y[0] = z + 4 * pio2_1t;
|
||||
y[1] = (z - y[0]) + 4 * pio2_1t;
|
||||
return -4;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
|
||||
medium:
|
||||
fn = __rint(x * invpio2);
|
||||
n = (INT32)fn;
|
||||
r = x - fn * pio2_1;
|
||||
w = fn * pio2_1t; /* 1st round, good to 85 bits */
|
||||
/* Matters with directed rounding. */
|
||||
if (r - w < -pio4) {
|
||||
n--;
|
||||
fn--;
|
||||
r = x - fn * pio2_1;
|
||||
w = fn * pio2_1t;
|
||||
} else if (r - w > pio4) {
|
||||
n++;
|
||||
fn++;
|
||||
r = x - fn * pio2_1;
|
||||
w = fn * pio2_1t;
|
||||
}
|
||||
y[0] = r - w;
|
||||
u.f = y[0];
|
||||
ey = u.i >> 52 & 0x7ff;
|
||||
ex = ix >> 20;
|
||||
if (ex - ey > 16) { /* 2nd round, good to 118 bits */
|
||||
t = r;
|
||||
w = fn * pio2_2;
|
||||
r = t - w;
|
||||
w = fn * pio2_2t - ((t - r) - w);
|
||||
y[0] = r - w;
|
||||
u.f = y[0];
|
||||
ey = u.i >> 52 & 0x7ff;
|
||||
if (ex - ey > 49) { /* 3rd round, good to 151 bits, covers all cases */
|
||||
t = r;
|
||||
w = fn * pio2_3;
|
||||
r = t - w;
|
||||
w = fn * pio2_3t - ((t - r) - w);
|
||||
y[0] = r - w;
|
||||
}
|
||||
}
|
||||
y[1] = (r - y[0]) - w;
|
||||
return n;
|
||||
}
|
||||
/*
|
||||
* all other (large) arguments
|
||||
*/
|
||||
if (ix >= 0x7ff00000) { /* x is inf or NaN */
|
||||
y[0] = y[1] = x - x;
|
||||
return 0;
|
||||
}
|
||||
/* set z = scalbn(|x|,-ilogb(x)+23) */
|
||||
u.f = x;
|
||||
u.i &= (UINT64)-1 >> 12;
|
||||
u.i |= (UINT64)(0x3ff + 23) << 52;
|
||||
z = u.f;
|
||||
for (i = 0; i < 2; i++) {
|
||||
tx[i] = (double)(INT32)z;
|
||||
z = (z - tx[i]) * 0x1p24;
|
||||
}
|
||||
tx[i] = z;
|
||||
/* skip zero terms, first term is non-zero */
|
||||
while (tx[i] == 0.0)
|
||||
i--;
|
||||
n = __rem_pio2_large(tx, ty, (int)(ix >> 20) - (0x3ff + 23), i + 1, 1);
|
||||
if (sign) {
|
||||
y[0] = -ty[0];
|
||||
y[1] = -ty[1];
|
||||
return -n;
|
||||
}
|
||||
y[0] = ty[0];
|
||||
y[1] = ty[1];
|
||||
return n;
|
||||
}
|
||||
extern int __rem_pio2(double x, double *y);
|
||||
|
||||
/* Copied from musl: src/math/__sin.c */
|
||||
static double __sin(double x, double y, int iy)
|
||||
|
|
@ -4962,14 +4562,6 @@ void __cdecl __libm_sse2_sqrt_precise(void)
|
|||
|
||||
#if _MSVCR_VER>=120
|
||||
|
||||
/*********************************************************************
|
||||
* rint (MSVCR120.@)
|
||||
*/
|
||||
double CDECL rint(double x)
|
||||
{
|
||||
return __rint(x);
|
||||
}
|
||||
|
||||
/*********************************************************************
|
||||
* rintf (MSVCR120.@)
|
||||
*
|
||||
|
|
|
|||
|
|
@ -119,7 +119,7 @@ int __rem_pio2(double x, double *y)
|
|||
if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
|
||||
medium:
|
||||
/* rint(x/(pi/2)) */
|
||||
fn = (double_t)x*invpio2 + toint - toint;
|
||||
fn = rint(x * invpio2);
|
||||
n = (int32_t)fn;
|
||||
r = x - fn*pio2_1;
|
||||
w = fn*pio2_1t; /* 1st round, good to 85 bits */
|
||||
|
|
|
|||
|
|
@ -52,7 +52,7 @@ int __rem_pio2f(float x, double *y)
|
|||
/* 25+53 bit pi is good enough for medium size */
|
||||
if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */
|
||||
/* Use a specialized rint() to get fn. */
|
||||
fn = (double_t)x*invpio2 + toint - toint;
|
||||
fn = fp_barrier(x * invpio2 + toint) - toint;
|
||||
n = (int32_t)fn;
|
||||
*y = x - fn*pio2_1 - fn*pio2_1t;
|
||||
/* Matters with directed rounding. */
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user