rt/portable/
math.rs

1// Portable replacement for the native x87 entry points. Decimal scaling emulates the
2// x87 64-bit significand after each multiply/divide. Transcendentals use explicit scalar
3// binary64 series; native x87 remains unchanged and is the comparison reference.
4
5use super::fx;
6
7const PI: f64 = core::f64::consts::PI;
8const HALF_PI: f64 = core::f64::consts::FRAC_PI_2;
9const LN2: f64 = core::f64::consts::LN_2;
10const POW10: [f64; 23] = [
11    1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16,
12    1e17, 1e18, 1e19, 1e20, 1e21, 1e22,
13];
14
15#[derive(Clone, Copy)]
16struct Extended {
17    m: u64,
18    e: i32,
19}
20
21fn round_shift(v: u128, shift: u32) -> u128 {
22    if shift == 0 {
23        return v;
24    }
25    if shift > 128 {
26        return 0;
27    }
28    if shift == 128 {
29        return (v > 1u128 << 127) as u128;
30    }
31    let q = v >> shift;
32    let remainder = v & ((1u128 << shift) - 1);
33    let half = 1u128 << (shift - 1);
34    q + (remainder > half || (remainder == half && q & 1 != 0)) as u128
35}
36
37impl Extended {
38    fn integer(v: u64) -> Self {
39        if v == 0 {
40            return Self { m: 0, e: 0 };
41        }
42        let sh = v.leading_zeros();
43        Self {
44            m: v << sh,
45            e: -(sh as i32),
46        }
47    }
48    fn from_f64(v: f64) -> Self {
49        let bits = v.to_bits();
50        let exp = ((bits >> 52) & 2047) as i32;
51        let m = bits & ((1u64 << 52) - 1);
52        if exp == 0 {
53            let mut n = Self::integer(m);
54            n.e -= 1074;
55            n
56        } else {
57            Self {
58                m: ((1u64 << 52) | m) << 11,
59                e: exp - 1023 - 63,
60            }
61        }
62    }
63    fn mul(self, b: Self) -> Self {
64        if self.m == 0 || b.m == 0 {
65            return Self::integer(0);
66        }
67        let product = self.m as u128 * b.m as u128;
68        let sh = 64 - product.leading_zeros();
69        let mut m = round_shift(product, sh);
70        let mut e = self.e + b.e + sh as i32;
71        if m > u64::MAX as u128 {
72            m >>= 1;
73            e += 1;
74        }
75        Self { m: m as u64, e }
76    }
77    fn div(self, b: Self) -> Self {
78        if self.m == 0 {
79            return self;
80        }
81        let sh = if self.m >= b.m { 63 } else { 64 };
82        let numerator = (self.m as u128) << sh;
83        let mut m = numerator / b.m as u128;
84        let remainder = numerator % b.m as u128;
85        if remainder * 2 > b.m as u128 || (remainder * 2 == b.m as u128 && m & 1 != 0) {
86            m += 1;
87        }
88        let mut e = self.e - b.e - sh;
89        if m > u64::MAX as u128 {
90            m >>= 1;
91            e += 1;
92        }
93        Self { m: m as u64, e }
94    }
95    fn double(self) -> f64 {
96        if self.m == 0 {
97            return 0.0;
98        }
99        let mut exponent = self.e + 63;
100        if exponent > 1023 {
101            return f64::INFINITY;
102        }
103        if exponent < -1022 {
104            let sh = -(self.e + 1074);
105            return f64::from_bits(if sh > 128 {
106                0
107            } else {
108                round_shift(self.m as u128, sh.max(0) as u32) as u64
109            });
110        }
111        let mut m = round_shift(self.m as u128, 11) as u64;
112        if m == 1u64 << 53 {
113            m >>= 1;
114            exponent += 1;
115        }
116        if exponent > 1023 {
117            return f64::INFINITY;
118        }
119        f64::from_bits(((exponent + 1023) as u64) << 52 | (m & ((1u64 << 52) - 1)))
120    }
121    fn rounded_integer(self, negative: bool) -> i64 {
122        let m = if self.e >= 0 {
123            if self.e >= 64 {
124                return i64::MIN;
125            }
126            (self.m as u128) << self.e
127        } else {
128            round_shift(self.m as u128, self.e.unsigned_abs())
129        };
130        if negative {
131            if m > 1u128 << 63 {
132                i64::MIN
133            } else {
134                (m as i64).wrapping_neg()
135            }
136        } else if m > i64::MAX as u128 {
137            i64::MIN
138        } else {
139            m as i64
140        }
141    }
142}
143
144fn power10(mut n: u32) -> Extended {
145    let mut p = Extended::integer(1);
146    while n > 22 {
147        p = p.mul(Extended::from_f64(POW10[22]));
148        n -= 22;
149    }
150    p.mul(Extended::from_f64(POW10[n as usize]))
151}
152
153/// Decimal scaling with the same 64-bit significand intermediate rounding as x87.
154pub fn scale(m: u64, e: i32) -> f64 {
155    let signed = m as i64;
156    if signed == 0 {
157        return 0.0;
158    }
159    if e > 5000 {
160        return if signed < 0 {
161            f64::NEG_INFINITY
162        } else {
163            f64::INFINITY
164        };
165    }
166    if e < -5000 {
167        return if signed < 0 { -0.0 } else { 0.0 };
168    }
169    let a = Extended::integer(signed.unsigned_abs());
170    let p = power10(e.unsigned_abs());
171    let out = if e < 0 { a.div(p) } else { a.mul(p) }.double();
172    if signed < 0 { -out } else { out }
173}
174/// Nearest-even integer after decimal scaling, including native indefinite overflow.
175pub fn scale_round(v: f64, k: i32) -> i64 {
176    if !v.is_finite() || k > 5000 {
177        return i64::MIN;
178    }
179    if k < -5000 || v == 0.0 {
180        return 0;
181    }
182    let a = Extended::from_f64(v.abs());
183    let p = power10(k.unsigned_abs());
184    if k < 0 { a.div(p) } else { a.mul(p) }.rounded_integer(v.is_sign_negative())
185}
186/// Nearest integer decimal logarithm; nonpositive inputs produce the SSE indefinite value.
187pub fn log10_round(v: f64) -> i32 {
188    fx::cvtd(log2(v) * core::f64::consts::LOG10_2)
189}
190
191/// Binary logarithm using a normalized mantissa and a convergent atanh series.
192pub fn log2(x: f64) -> f64 {
193    if x.is_nan() {
194        return x + x;
195    }
196    if x == 0.0 {
197        return f64::NEG_INFINITY;
198    }
199    if x < 0.0 {
200        return f64::from_bits(0xfff8000000000000);
201    }
202    if x == f64::INFINITY {
203        return x;
204    }
205    let ext = Extended::from_f64(x);
206    let mut m = ext.m as f64 / 9223372036854775808.0;
207    let mut e = ext.e + 63;
208    if m > core::f64::consts::SQRT_2 {
209        m *= 0.5;
210        e += 1;
211    }
212    let z = (m - 1.0) / (m + 1.0);
213    let zz = z * z;
214    let mut term = z;
215    let mut sum = z;
216    let mut i = 3;
217    while i <= 39 {
218        term *= zz;
219        sum += term / i as f64;
220        i += 2;
221    }
222    e as f64 + (2.0 * sum) / LN2
223}
224/// Platform-independent base-two exponent; finite underflow rounds into subnormals.
225pub fn exp2(x: f64) -> f64 {
226    if x.is_nan() {
227        return x + x;
228    }
229    if x >= 1024.0 {
230        return f64::INFINITY;
231    }
232    if x <= -1075.0 {
233        return 0.0;
234    }
235    let n = fx::cvtd(x);
236    let f = (x - n as f64) * LN2;
237    let mut sum = 1.0;
238    let mut term = 1.0;
239    let mut i = 1;
240    while i <= 20 {
241        term *= f / i as f64;
242        sum += term;
243        i += 1;
244    }
245    if n > 1023 {
246        (sum * 2.0) * f64::from_bits(2046u64 << 52)
247    } else if n < -1022 {
248        (sum * f64::from_bits(((n + 1077) as u64) << 52)) * f64::from_bits(969u64 << 52)
249    } else {
250        sum * f64::from_bits(((n + 1023) as u64) << 52)
251    }
252}
253/// Existing rotation/zoom base-two logarithm API.
254pub fn rz_log2(x: f64) -> f64 {
255    log2(x)
256}
257/// Existing rotation/zoom base-two exponent API.
258pub fn rz_exp2(x: f64) -> f64 {
259    exp2(x)
260}
261/// Effect power contract: y=0 yields 1, nonpositive bases yield 0, NaNs propagate canonically.
262pub fn pow(x: f64, y: f64) -> f64 {
263    if y.is_nan() {
264        return f64::NAN;
265    }
266    if y == 0.0 {
267        return 1.0;
268    }
269    if x.is_nan() {
270        return f64::NAN;
271    }
272    if x <= 0.0 {
273        return 0.0;
274    }
275    if y == 1.0 {
276        return x;
277    }
278    exp2(y * log2(x))
279}
280
281fn atan_reduced(mut x: f64) -> f64 {
282    let inverse = x > 1.0;
283    if inverse {
284        x = 1.0 / x;
285    }
286    let transformed = x > 0.41421356237309503;
287    if transformed {
288        x = (x - 1.0) / (x + 1.0);
289    }
290    let xx = -x * x;
291    let mut term = x;
292    let mut sum = x;
293    let mut i = 3;
294    while i <= 61 {
295        term *= xx;
296        sum += term / i as f64;
297        i += 2;
298    }
299    if transformed {
300        sum += core::f64::consts::FRAC_PI_4;
301    }
302    if inverse { HALF_PI - sum } else { sum }
303}
304
305fn remainder(x: f64, divisor: f64) -> f64 {
306    if x.abs() < divisor {
307        return x;
308    }
309    let a = Extended::from_f64(x.abs());
310    let b = Extended::from_f64(divisor);
311    let mut m = a.m as u128;
312    let mut e = a.e;
313    while e > b.e {
314        if m >= b.m as u128 {
315            m -= b.m as u128;
316        }
317        m <<= 1;
318        e -= 1;
319    }
320    if m >= b.m as u128 {
321        m -= b.m as u128;
322    }
323    let mut result = Extended::integer(m as u64);
324    result.e += b.e;
325    let v = result.double();
326    if x.is_sign_negative() { -v } else { v }
327}
328/// Quadrant-aware arctangent, including signed zeros and infinite operands.
329pub fn atan2(y: f64, x: f64) -> f64 {
330    if y.is_nan() || x.is_nan() {
331        return x + y;
332    }
333    let ax = x.abs();
334    let ay = y.abs();
335    let mut a = if ay == 0.0 {
336        0.0
337    } else if ax == 0.0 {
338        HALF_PI
339    } else if ax.is_infinite() && ay.is_infinite() {
340        core::f64::consts::FRAC_PI_4
341    } else {
342        atan_reduced(ay / ax)
343    };
344    if x.is_sign_negative() {
345        a = PI - a;
346    }
347    if y.is_sign_negative() { -a } else { a }
348}
349
350/// Cosine and sine with shared quadrant reduction; return order matches the x87 API.
351pub fn sincos(x: f64) -> (f64, f64) {
352    if !x.is_finite() {
353        return (f64::NAN, f64::NAN);
354    }
355    if x == 0.0 {
356        return (1.0, x);
357    }
358    let reduced = remainder(x, 2.0 * PI);
359    let quadrant = fx::cvtd(reduced / HALF_PI);
360    // Split pi/2 minimizes cancellation for the normal UI/effect angle range.
361    let r =
362        (reduced - quadrant as f64 * 1.5707963267948966) - quadrant as f64 * 6.123233995736766e-17;
363    let rr = -r * r;
364    let mut s = r;
365    let mut st = r;
366    let mut c = 1.0;
367    let mut ct = 1.0;
368    let mut i = 1;
369    while i <= 10 {
370        st *= rr / ((2 * i) * (2 * i + 1)) as f64;
371        ct *= rr / ((2 * i - 1) * (2 * i)) as f64;
372        s += st;
373        c += ct;
374        i += 1;
375    }
376    match quadrant.rem_euclid(4) {
377        0 => (c, s),
378        1 => (-s, c),
379        2 => (-c, -s),
380        _ => (s, -c),
381    }
382}
383/// Scalar cosine using the common angle reduction.
384pub fn cos(x: f64) -> f64 {
385    sincos(x).0
386}
387/// Scalar sine using the common angle reduction.
388pub fn sin(x: f64) -> f64 {
389    sincos(x).1
390}
rt/portable/math.rs

rt/portable/
math.rs
