Skip to content

Commit cd84aa7

Browse files
OkladnojOkladnoj
committed
feat: Replace hardware-dependent trig with fdlibm
Replace all platform-specific math in WWMath (x87 FPU asm on Win32, system libm on macOS/64-bit) with Sun's fdlibm 5.3 library. This ensures bit-identical IEEE 754 floating-point results across all platforms, eliminating the root cause of multiplayer desyncs. Functions replaced: Sin, Cos, Sqrt, Acos, Asin, Atan, Atan2, Inv_Sqrt. Functions preserved: Fast_Sin, Fast_Cos (LUT-based).
1 parent 6266009 commit cd84aa7

14 files changed

Lines changed: 1970 additions & 137 deletions

File tree

‎Core/Libraries/Source/WWVegas/WWMath/CMakeLists.txt‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,19 @@ set(WWMATH_SRC
7878
wwmath.cpp
7979
wwmath.h
8080
wwmathids.h
81+
82+
fdlibm/fdlibm_det.h
83+
fdlibm/fdlibm.h
84+
fdlibm/fdlibm_det.c
85+
fdlibm/e_acos.c
86+
fdlibm/e_asin.c
87+
fdlibm/e_atan2.c
88+
fdlibm/e_rem_pio2.c
89+
fdlibm/e_sqrt.c
90+
fdlibm/k_cos.c
91+
fdlibm/k_rem_pio2.c
92+
fdlibm/k_sin.c
93+
fdlibm/k_tan.c
8194
)
8295

8396
add_library(core_wwmath STATIC)

‎Core/Libraries/Source/WWVegas/WWMath/fdlibm/e_acos.c‎

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
2+
/* @(#)e_acos.c 1.3 95/01/18 */
3+
/*
4+
* ====================================================
5+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6+
*
7+
* Developed at SunSoft, a Sun Microsystems, Inc. business.
8+
* Permission to use, copy, modify, and distribute this
9+
* software is freely granted, provided that this notice
10+
* is preserved.
11+
* ====================================================
12+
*/
13+
14+
/* __ieee754_acos(x)
15+
* Method :
16+
* acos(x) = pi/2 - asin(x)
17+
* acos(-x) = pi/2 + asin(x)
18+
* For |x|<=0.5
19+
* acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
20+
* For x>0.5
21+
* acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
22+
* = 2asin(sqrt((1-x)/2))
23+
* = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
24+
* = 2f + (2c + 2s*z*R(z))
25+
* where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
26+
* for f so that f+c ~ sqrt(z).
27+
* For x<-0.5
28+
* acos(x) = pi - 2asin(sqrt((1-|x|)/2))
29+
* = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
30+
*
31+
* Special cases:
32+
* if x is NaN, return x itself;
33+
* if |x|>1, return NaN with invalid signal.
34+
*
35+
* Function needed: sqrt
36+
*/
37+
38+
#include "fdlibm.h"
39+
40+
#ifdef __STDC__
41+
static const double
42+
#else
43+
static double
44+
#endif
45+
one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
46+
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
47+
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
48+
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
49+
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
50+
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
51+
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
52+
pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
53+
pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
54+
pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
55+
qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
56+
qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
57+
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
58+
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
59+
60+
#ifdef __STDC__
61+
double __ieee754_acos(double x)
62+
#else
63+
double __ieee754_acos(x)
64+
double x;
65+
#endif
66+
{
67+
double z,p,q,r,w,s,c,df;
68+
int hx,ix;
69+
hx = __HI(x);
70+
ix = hx&0x7fffffff;
71+
if(ix>=0x3ff00000) { /* |x| >= 1 */
72+
if(((ix-0x3ff00000)|__LO(x))==0) { /* |x|==1 */
73+
if(hx>0) return 0.0; /* acos(1) = 0 */
74+
else return pi+2.0*pio2_lo; /* acos(-1)= pi */
75+
}
76+
return (x-x)/(x-x); /* acos(|x|>1) is NaN */
77+
}
78+
if(ix<0x3fe00000) { /* |x| < 0.5 */
79+
if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
80+
z = x*x;
81+
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
82+
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
83+
r = p/q;
84+
return pio2_hi - (x - (pio2_lo-x*r));
85+
} else if (hx<0) { /* x < -0.5 */
86+
z = (one+x)*0.5;
87+
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
88+
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
89+
s = sqrt(z);
90+
r = p/q;
91+
w = r*s-pio2_lo;
92+
return pi - 2.0*(s+w);
93+
} else { /* x > 0.5 */
94+
z = (one-x)*0.5;
95+
s = sqrt(z);
96+
df = s;
97+
__LO(df) = 0;
98+
c = (z-df*df)/(s+df);
99+
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
100+
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
101+
r = p/q;
102+
w = r*s+c;
103+
return 2.0*(df+w);
104+
}
105+
}

‎Core/Libraries/Source/WWVegas/WWMath/fdlibm/e_asin.c‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
2+
/* @(#)e_asin.c 1.3 95/01/18 */
3+
/*
4+
* ====================================================
5+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6+
*
7+
* Developed at SunSoft, a Sun Microsystems, Inc. business.
8+
* Permission to use, copy, modify, and distribute this
9+
* software is freely granted, provided that this notice
10+
* is preserved.
11+
* ====================================================
12+
*/
13+
14+
/* __ieee754_asin(x)
15+
* Method :
16+
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
17+
* we approximate asin(x) on [0,0.5] by
18+
* asin(x) = x + x*x^2*R(x^2)
19+
* where
20+
* R(x^2) is a rational approximation of (asin(x)-x)/x^3
21+
* and its remez error is bounded by
22+
* |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
23+
*
24+
* For x in [0.5,1]
25+
* asin(x) = pi/2-2*asin(sqrt((1-x)/2))
26+
* Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
27+
* then for x>0.98
28+
* asin(x) = pi/2 - 2*(s+s*z*R(z))
29+
* = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
30+
* For x<=0.98, let pio4_hi = pio2_hi/2, then
31+
* f = hi part of s;
32+
* c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
33+
* and
34+
* asin(x) = pi/2 - 2*(s+s*z*R(z))
35+
* = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
36+
* = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
37+
*
38+
* Special cases:
39+
* if x is NaN, return x itself;
40+
* if |x|>1, return NaN with invalid signal.
41+
*
42+
*/
43+
44+
45+
#include "fdlibm.h"
46+
47+
#ifdef __STDC__
48+
static const double
49+
#else
50+
static double
51+
#endif
52+
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
53+
huge = 1.000e+300,
54+
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
55+
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
56+
pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
57+
/* coefficient for R(x^2) */
58+
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
59+
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
60+
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
61+
pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
62+
pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
63+
pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
64+
qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
65+
qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
66+
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
67+
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
68+
69+
#ifdef __STDC__
70+
double __ieee754_asin(double x)
71+
#else
72+
double __ieee754_asin(x)
73+
double x;
74+
#endif
75+
{
76+
double t,w,p,q,c,r,s;
77+
int hx,ix;
78+
hx = __HI(x);
79+
ix = hx&0x7fffffff;
80+
if(ix>= 0x3ff00000) { /* |x|>= 1 */
81+
if(((ix-0x3ff00000)|__LO(x))==0)
82+
/* asin(1)=+-pi/2 with inexact */
83+
return x*pio2_hi+x*pio2_lo;
84+
return (x-x)/(x-x); /* asin(|x|>1) is NaN */
85+
} else if (ix<0x3fe00000) { /* |x|<0.5 */
86+
if(ix<0x3e400000) { /* if |x| < 2**-27 */
87+
if(huge+x>one) return x;/* return x with inexact if x!=0*/
88+
} else
89+
t = x*x;
90+
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
91+
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
92+
w = p/q;
93+
return x+x*w;
94+
}
95+
/* 1> |x|>= 0.5 */
96+
w = one-fabs(x);
97+
t = w*0.5;
98+
p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
99+
q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
100+
s = sqrt(t);
101+
if(ix>=0x3FEF3333) { /* if |x| > 0.975 */
102+
w = p/q;
103+
t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
104+
} else {
105+
w = s;
106+
__LO(w) = 0;
107+
c = (t-w*w)/(s+w);
108+
r = p/q;
109+
p = 2.0*s*r-(pio2_lo-2.0*c);
110+
q = pio4_hi-2.0*w;
111+
t = pio4_hi-(p-q);
112+
}
113+
if(hx>0) return t; else return -t;
114+
}

‎Core/Libraries/Source/WWVegas/WWMath/fdlibm/e_atan2.c‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
2+
/* @(#)e_atan2.c 1.3 95/01/18 */
3+
/*
4+
* ====================================================
5+
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6+
*
7+
* Developed at SunSoft, a Sun Microsystems, Inc. business.
8+
* Permission to use, copy, modify, and distribute this
9+
* software is freely granted, provided that this notice
10+
* is preserved.
11+
* ====================================================
12+
*
13+
*/
14+
15+
/* __ieee754_atan2(y,x)
16+
* Method :
17+
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
18+
* 2. Reduce x to positive by (if x and y are unexceptional):
19+
* ARG (x+iy) = arctan(y/x) ... if x > 0,
20+
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
21+
*
22+
* Special cases:
23+
*
24+
* ATAN2((anything), NaN ) is NaN;
25+
* ATAN2(NAN , (anything) ) is NaN;
26+
* ATAN2(+-0, +(anything but NaN)) is +-0 ;
27+
* ATAN2(+-0, -(anything but NaN)) is +-pi ;
28+
* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
29+
* ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
30+
* ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
31+
* ATAN2(+-INF,+INF ) is +-pi/4 ;
32+
* ATAN2(+-INF,-INF ) is +-3pi/4;
33+
* ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
34+
*
35+
* Constants:
36+
* The hexadecimal values are the intended ones for the following
37+
* constants. The decimal values may be used, provided that the
38+
* compiler will convert from decimal to binary accurately enough
39+
* to produce the hexadecimal values shown.
40+
*/
41+
42+
#include "fdlibm.h"
43+
44+
#ifdef __STDC__
45+
static const double
46+
#else
47+
static double
48+
#endif
49+
tiny = 1.0e-300,
50+
zero = 0.0,
51+
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
52+
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
53+
pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
54+
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
55+
56+
#ifdef __STDC__
57+
double __ieee754_atan2(double y, double x)
58+
#else
59+
double __ieee754_atan2(y,x)
60+
double y,x;
61+
#endif
62+
{
63+
double z;
64+
int k,m,hx,hy,ix,iy;
65+
unsigned lx,ly;
66+
67+
hx = __HI(x); ix = hx&0x7fffffff;
68+
lx = __LO(x);
69+
hy = __HI(y); iy = hy&0x7fffffff;
70+
ly = __LO(y);
71+
if(((ix|((lx|-lx)>>31))>0x7ff00000)||
72+
((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
73+
return x+y;
74+
if((hx-0x3ff00000|lx)==0) return atan(y); /* x=1.0 */
75+
m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
76+
77+
/* when y = 0 */
78+
if((iy|ly)==0) {
79+
switch(m) {
80+
case 0:
81+
case 1: return y; /* atan(+-0,+anything)=+-0 */
82+
case 2: return pi+tiny;/* atan(+0,-anything) = pi */
83+
case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
84+
}
85+
}
86+
/* when x = 0 */
87+
if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
88+
89+
/* when x is INF */
90+
if(ix==0x7ff00000) {
91+
if(iy==0x7ff00000) {
92+
switch(m) {
93+
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
94+
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
95+
case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
96+
case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
97+
}
98+
} else {
99+
switch(m) {
100+
case 0: return zero ; /* atan(+...,+INF) */
101+
case 1: return -zero ; /* atan(-...,+INF) */
102+
case 2: return pi+tiny ; /* atan(+...,-INF) */
103+
case 3: return -pi-tiny ; /* atan(-...,-INF) */
104+
}
105+
}
106+
}
107+
/* when y is INF */
108+
if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
109+
110+
/* compute y/x */
111+
k = (iy-ix)>>20;
112+
if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */
113+
else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
114+
else z=atan(fabs(y/x)); /* safe to do y/x */
115+
switch (m) {
116+
case 0: return z ; /* atan(+,+) */
117+
case 1: __HI(z) ^= 0x80000000;
118+
return z ; /* atan(-,+) */
119+
case 2: return pi-(z-pi_lo);/* atan(+,-) */
120+
default: /* case 3 */
121+
return (z-pi_lo)-pi;/* atan(-,-) */
122+
}
123+
}

0 commit comments

Comments
 (0)