brine_fp/
exp.rs

1use super::consts::*;
2use super::unsigned::UnsignedNumeric;
3use super::signed::SignedNumeric;
4
5// Modified from the original to support precise numbers instead of floats
6// origin: FreeBSD /usr/src/lib/msun/src/e_exp.c */
7// 
8// ====================================================
9// Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
10// 
11// Permission to use, copy, modify, and distribute this
12// software is freely granted, provided that this notice
13// is preserved.
14// ====================================================
15// 
16// exp(x)
17// Returns the exponential of x.
18// 
19// Method
20//   1. Argument reduction:
21//      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
22//      Given x, find r and integer k such that
23// 
24//               x = k*ln2 + r,  |r| <= 0.5*ln2.
25// 
26//      Here r will be represented as r = hi-lo for better
27//      accuracy.
28// 
29//   2. Approximation of exp(r) by a special rational function on
30//      the interval [0,0.34658]:
31//      Write
32//          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
33//      We use a special Remez algorithm on [0,0.34658] to generate
34//      a polynomial of degree 5 to approximate R. The maximum error
35//      of this polynomial approximation is bounded by 2**-59. In
36//      other words,
37//          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
38//      (where z=r*r, and the values of P1 to P5 are listed below)
39//      and
40//          |                  5          |     -59
41//          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
42//          |                             |
43//      The computation of exp(r) thus becomes
44//                              2*r
45//              exp(r) = 1 + ----------
46//                            R(r) - r
47//                                 r*c(r)
48//                     = 1 + r + ----------- (for better accuracy)
49//                                2 - c(r)
50//      where
51//                              2       4             10
52//              c(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
53// 
54//   3. Scale back to obtain exp(x):
55//      From step 1, we have
56//         exp(x) = 2^k * exp(r)
57// 
58// Special cases:
59//      exp(INF) is INF, exp(NaN) is NaN;
60//      exp(-INF) is 0, and
61//      for finite argument, only exp(0)=1 is exact.
62// 
63// Accuracy:
64//      according to an error analysis, the error is always less than
65//      1 ulp (unit in the last place).
66// 
67// Misc. info.
68//      For IEEE double
69//          if x >  709.782712893383973096 then exp(x) overflows
70//          if x < -745.133219101941108420 then exp(x) underflows
71
72impl SignedNumeric {
73
74    /// Calculate the exponential of `x`, that is, *e* raised to the power `x`
75    /// (where *e* is the base of the natural system of logarithms, approximately 2.71828).
76    /// Note that precision can start to get inaccurate for larger numbers (> 20).
77    pub fn exp(&self) -> Option<UnsignedNumeric> {
78        let hi: Self;
79        let lo: Self;
80        let k: Self;
81        let x: Self;
82
83        // argument reduction
84        // if |x| > 0.5 ln2
85        if self.value.greater_than(&HALFLN2) {
86            // if |x| >= 1.5 ln2
87            if self.value.greater_than_or_equal(&THREEHALFLN2) {
88                k = INVLN2
89                    .signed()
90                    .checked_mul(self)?
91                    .checked_add(&Self {
92                        value: HALF,
93                        is_negative: self.is_negative,
94                    })?
95                    .floor()?;
96
97                // A K larger than this value will cause less than 9 decimals of precision
98                // if k.value.to_imprecise()? > 29 {
99                //   return None
100                // }
101            } else {
102                k = Self {
103                    value: UnsignedNumeric::one(),
104                    is_negative: self.is_negative,
105                }
106            }
107            hi = self.checked_sub(
108                &k.checked_mul(&LN2HI.signed())?
109                    .checked_div(&LN2HI_SCALE.signed())?,
110            )?;
111
112            lo = k
113                .checked_mul(&LN2LO.signed())?
114                .checked_div(&LN2LO_SCALE.signed())?;
115            x = hi.checked_sub(&lo)?
116        } else {
117            x = self.clone();
118            k = UnsignedNumeric::zero().signed();
119            hi = self.clone();
120            lo = UnsignedNumeric::zero().signed()
121        }
122
123        // x is now in primary range
124        let xx = x.checked_mul(&x)?;
125        // c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5))));
126        let p4p5 = P4.checked_add(&xx.checked_mul(&P5)?)?;
127        let p3p4p5 = P3.checked_add(&xx.checked_mul(&p4p5)?)?;
128        let p2p3p4p5 = P2.checked_add(&xx.checked_mul(&p3p4p5)?)?;
129        let p1p2p3p4p5 = P1.checked_add(&xx.checked_mul(&p2p3p4p5)?)?;
130        let c = x.checked_sub(&p1p2p3p4p5.checked_mul(&xx)?)?;
131
132        // y = 1. + (x * c / (2. - c) - lo + hi);
133        let y = ONE_PREC.signed().checked_add(
134            &x.checked_mul(&c)?
135                .checked_div(&TWO_PREC.signed().checked_sub(&c)?)?
136                .checked_sub(&lo)?
137                .checked_add(&hi)?,
138        )?;
139
140        if k.value.eq(&UnsignedNumeric::zero()) {
141            Some(y.value)
142        } else {
143            let bits = k.value.to_imprecise()?;
144
145            if k.is_negative {
146                Some(UnsignedNumeric {
147                    value: y.value.value >> bits,
148                })
149            } else {
150                Some(UnsignedNumeric {
151                    value: y.value.value << bits,
152                })
153            }
154        }
155    }
156
157    /// Returns the square root of `self`, returning `None` if the number is negative.
158    pub fn sqrt(&self) -> Option<Self> {
159        if self.is_negative {
160            return None;
161        }
162        self.value.sqrt().map(|v| Self { value: v, is_negative: false })
163    }
164}
165
166impl UnsignedNumeric {
167
168    /// Frexp breaks f into a normalized fraction and an integral power of two. It returns frac and
169    /// exp satisfying f == frac × 2**exp, with the absolute value of frac in the interval [½, 1).
170    ///
171    /// Special cases are:
172    ///	Frexp(±0) = ±0, 0
173    ///	Frexp(±Inf) = ±Inf, 0
174    ///	Frexp(NaN) = NaN, 0
175    pub fn frexp(&self) -> Option<(Self, i64)> {
176        if self.eq(&ZERO_PREC) {
177            Some((ZERO_PREC.clone(), 0))
178        } else if self.less_than(&ONE_PREC) {
179            let first_leading = self.value.0[0].leading_zeros();
180            let one_leading = ONE_PREC.value.0[0].leading_zeros();
181            let bits = i64::from(first_leading.checked_sub(one_leading).unwrap());
182            let frac = UnsignedNumeric {
183                value: self.value << bits,
184            };
185            if frac.less_than(&HALF) {
186                Some((frac.checked_mul(&TWO_PREC).unwrap(), -bits - 1))
187            } else {
188                Some((frac, -bits))
189            }
190        } else {
191            let bits = 128_i64.checked_sub(i64::from(self.to_imprecise()?.leading_zeros()))?;
192            let frac = UnsignedNumeric {
193                value: self.value >> bits,
194            };
195            if frac.less_than(&HALF) {
196                Some((frac.checked_mul(&TWO_PREC).unwrap(), bits - 1))
197            } else {
198                Some((frac, bits))
199            }
200        }
201    }
202
203    /// Raises `self` to the power of `exp`, returning the result as a new `UnsignedNumeric`.
204    /// Returns `None` if the operation would overflow or if `self` is zero.
205    ///
206    /// b = pow/frac
207    /// y = a^b
208    /// ln (y) = bln (a)
209    /// y = e^(b ln (a))
210    pub fn pow(&self, exp: &Self) -> Option<Self> {
211        if self.eq(&ZERO_PREC) {
212            return Some(ZERO_PREC.clone());
213        }
214
215        let lg = self.log()?;
216        let x = exp.signed().checked_mul(&lg)?;
217        x.exp()
218    }
219
220    /// Returns the square root of `self`
221    pub fn sqrt(&self) -> Option<Self> {
222        self.pow(&HALF)
223    }
224}
225
226#[cfg(test)]
227mod tests {
228    use super::*;
229    use crate::InnerUint;
230
231    fn frexp_recombine(frac: UnsignedNumeric, exp: i64) -> UnsignedNumeric {
232        let shifted = if exp >= 0 {
233            UnsignedNumeric {
234                value: frac.value << (exp as usize),
235            }
236        } else {
237            UnsignedNumeric {
238                value: frac.value >> ((-exp) as usize),
239            }
240        };
241        shifted
242    }
243
244    #[test]
245    fn test_signed_exp() {
246        let precision = InnerUint::from(1_000_000_000_u128); // correct to at least 9 decimal places
247
248        let half = UnsignedNumeric { value: half() }.signed();
249        assert!(half.exp().unwrap().almost_eq(
250            &UnsignedNumeric::new(16487212707001282)
251                .checked_div(&UnsignedNumeric::new(10000000000000000))
252                .unwrap(),
253            precision
254        ));
255
256        let three_half = UnsignedNumeric::new(15)
257            .checked_div(&UnsignedNumeric::new(10))
258            .unwrap()
259            .signed();
260        assert!(three_half.exp().unwrap().almost_eq(
261            &UnsignedNumeric::new(44816890703380645)
262                .checked_div(&UnsignedNumeric::new(10000000000000000))
263                .unwrap(),
264            precision
265        ));
266
267        let point_one = UnsignedNumeric::new(1)
268            .checked_div(&UnsignedNumeric::new(10))
269            .unwrap()
270            .signed();
271        assert!(point_one.exp().unwrap().almost_eq(
272            &UnsignedNumeric::new(11051709180756477)
273                .checked_div(&UnsignedNumeric::new(10000000000000000))
274                .unwrap(),
275            precision
276        ));
277
278        let negative = UnsignedNumeric::new(55)
279            .checked_div(&UnsignedNumeric::new(100))
280            .unwrap()
281            .signed()
282            .negate();
283        assert!(negative.exp().unwrap().almost_eq(
284            &UnsignedNumeric::new(5769498103804866)
285                .checked_div(&UnsignedNumeric::new(10000000000000000))
286                .unwrap(),
287            precision
288        ));
289
290        let test = UnsignedNumeric::new(19).signed();
291        assert!(test.exp().unwrap().almost_eq(
292            &UnsignedNumeric::new(178482300963187260)
293                .checked_div(&UnsignedNumeric::new(1000000000))
294                .unwrap(),
295            precision
296        ));
297    }
298
299    #[test]
300    fn test_pow() {
301        let precision = InnerUint::from(5_000_000_000_000_u128); // correct to at least 12 decimal places
302        let test = UnsignedNumeric::new(8);
303        let sqrt = test.pow(&HALF).unwrap();
304        let expected = UnsignedNumeric::new(28284271247461903)
305            .checked_div(&UnsignedNumeric::new(10000000000000000))
306            .unwrap();
307        assert!(sqrt.almost_eq(&expected, precision));
308
309        let test2 = UnsignedNumeric::new(55)
310            .checked_div(&UnsignedNumeric::new(100))
311            .unwrap();
312        let squared = test2.pow(&TWO_PREC).unwrap();
313        let expected = UnsignedNumeric::new(3025)
314            .checked_div(&UnsignedNumeric::new(10000))
315            .unwrap();
316        assert!(squared.almost_eq(&expected, precision));
317    }
318
319    #[test]
320    fn test_sqrt() {
321        let precision = InnerUint::from(5_000_000_000_000_u128); // correct to at least 12 decimal places
322        let test = UnsignedNumeric::new(12);
323        let sqrt = test.sqrt().unwrap();
324        let expected = UnsignedNumeric::new(34641016151377544)
325            .checked_div(&UnsignedNumeric::new(10000000000000000))
326            .unwrap();
327        assert!(sqrt.almost_eq(&expected, precision));
328    }
329
330    #[test]
331    pub fn test_signed_sqrt() {
332        let precision = InnerUint::from(5_000_000_000_000_u128); // correct to at least 12 decimal places
333        let test = SignedNumeric {
334            value: UnsignedNumeric::new(8),
335            is_negative: false,
336        };
337        let sqrt = test.sqrt().unwrap();
338        let expected = UnsignedNumeric::new(28284271247461903)
339            .checked_div(&UnsignedNumeric::new(10000000000000000))
340            .unwrap();
341        assert!(sqrt.value.almost_eq(&expected, precision));
342
343        let neg_test = SignedNumeric {
344            value: UnsignedNumeric::new(8),
345            is_negative: true,
346        };
347        assert!(neg_test.sqrt().is_none());
348    }
349
350    #[test]
351    fn test_exp_zero() {
352        let zero = UnsignedNumeric::zero().signed();
353        let result = zero.exp().unwrap();
354        assert_eq!(result, UnsignedNumeric::one());
355    }
356
357    #[test]
358    fn test_exp_small_negative() {
359        let x = UnsignedNumeric::new(1)
360            .checked_div(&UnsignedNumeric::new(1000))
361            .unwrap()
362            .signed()
363            .negate();
364        let result = x.exp().unwrap();
365
366        // e^-0.001 ≈ 0.999000499833375 (scaled: 999000499833375000)
367        let expected = UnsignedNumeric::from_scaled_u128(999_000_499_833_375_000);
368        assert!(
369            result.almost_eq(&expected, InnerUint::from(10_000_000)),
370            "got: {} expected: {}",
371            result.to_string(),
372            expected.to_string()
373        );
374    }
375
376    #[test]
377    fn test_exp_large_positive() {
378        let x = UnsignedNumeric::new(10).signed(); // e^10 ≈ 22026.4657948067
379        let result = x.exp().unwrap();
380
381        // Correctly scaled expected value: e^10 * 1e18
382        let expected = UnsignedNumeric::from_scaled_u128(22_026_465_794_806_700_000_000);
383        assert!(
384            result.almost_eq(&expected, InnerUint::from(1_000_000_000_000_u64)),
385            "got: {} expected: {}",
386            result.to_string(),
387            expected.to_string()
388        );
389    }
390
391    #[test]
392    fn test_exp_large_negative() {
393        let x = UnsignedNumeric::new(10).signed().negate(); // e^-10 ≈ 0.00004539992
394        let result = x.exp().unwrap();
395
396        let expected = UnsignedNumeric::from_scaled_u128(45_399_920_000_000); // scaled by 10^18
397        assert!(
398            result.almost_eq(&expected, InnerUint::from(1_000_000_000)),
399            "got: {} expected: {}",
400            result.to_string(),
401            expected.to_string()
402        );
403    }
404
405    #[test]
406    fn test_frexp_zero() {
407        let zero = UnsignedNumeric::zero();
408        let (frac, exp) = zero.frexp().unwrap();
409        assert_eq!(frac, zero);
410        assert_eq!(exp, 0);
411    }
412
413    #[test]
414    fn test_frexp_one() {
415        let one = UnsignedNumeric::one();
416        let (frac, exp) = one.frexp().unwrap();
417        // Should return exactly 0.5 ≤ frac < 1, and 1.0 == frac × 2^exp
418        let recombined = frexp_recombine(frac.clone(), exp);
419        assert!(
420            recombined.almost_eq(&one, InnerUint::from(1_000_000_000)),
421            "Expected: {}, Got: {} × 2^{} = {}",
422            one.to_string(),
423            frac.to_string(),
424            exp,
425            recombined.to_string()
426        );
427        assert!(frac.greater_than_or_equal(&HALF));
428        assert!(frac.less_than(&ONE_PREC));
429    }
430
431    #[test]
432    fn test_frexp_two() {
433        let two = UnsignedNumeric::new(2);
434        let (frac, exp) = two.frexp().unwrap();
435        let recombined = frexp_recombine(frac.clone(), exp);
436        assert!(
437            recombined.almost_eq(&two, InnerUint::from(1_000_000_000)),
438            "Expected: {}, Got: {} × 2^{} = {}",
439            two.to_string(),
440            frac.to_string(),
441            exp,
442            recombined.to_string()
443        );
444        assert!(frac.greater_than_or_equal(&HALF));
445        assert!(frac.less_than(&ONE_PREC));
446    }
447
448    #[test]
449    fn test_frexp_fractional() {
450        let val = UnsignedNumeric::new(3)
451            .checked_div(&UnsignedNumeric::new(4)) // 0.75
452            .unwrap();
453        let (frac, exp) = val.frexp().unwrap();
454        let recombined = frexp_recombine(frac.clone(), exp);
455
456        assert!(
457            recombined.almost_eq(&val, InnerUint::from(1_000_000_000)),
458            "Expected: {}, Got: {} × 2^{} = {}",
459            val.to_string(),
460            frac.to_string(),
461            exp,
462            recombined.to_string()
463        );
464        assert!(frac.greater_than_or_equal(&HALF));
465        assert!(frac.less_than(&ONE_PREC));
466    }
467
468    #[test]
469    fn test_frexp_large_value() {
470        let val = UnsignedNumeric {
471            value: InnerUint([0, 0, 1]), // 2^128 scaled fixed-point
472        };
473        let (frac, exp) = val.frexp().unwrap();
474        let recombined = frexp_recombine(frac.clone(), exp);
475
476        assert!(
477            recombined.almost_eq(&val, InnerUint::from(10_000_000_000_u64)),
478            "Expected: {}, Got: {} × 2^{} = {}",
479            val.to_string(),
480            frac.to_string(),
481            exp,
482            recombined.to_string()
483        );
484    }
485}
brine_fp/exp.rs

brine_fp/
exp.rs
