pallas_math/
math.rs

1/*!
2# Cardano Math functions
3 */
4
5use std::fmt::{Debug, Display};
6use std::ops::{Div, Mul, Neg, Sub};
7use std::sync::LazyLock;
8use thiserror::Error;
9
10pub type FixedDecimal = crate::math_dashu::Decimal;
11
12pub static ZERO: LazyLock<FixedDecimal> = LazyLock::new(|| FixedDecimal::from(0u64));
13pub static MINUS_ONE: LazyLock<FixedDecimal> = LazyLock::new(|| FixedDecimal::from(-1i64));
14pub static ONE: LazyLock<FixedDecimal> = LazyLock::new(|| FixedDecimal::from(1u64));
15
16#[derive(Debug, Error)]
17pub enum Error {
18    #[error("error in regex")]
19    RegexFailure(#[from] regex::Error),
20
21    #[error("string contained a nul byte")]
22    NulFailure(#[from] std::ffi::NulError),
23}
24
25pub const DEFAULT_PRECISION: u64 = 34;
26
27pub trait FixedPrecision:
28    Neg + Mul + Div + Sub + Display + Clone + PartialEq + PartialOrd + Debug + From<u64> + From<i64>
29{
30    /// Creates a new fixed point number with the given precision
31    fn new(precision: u64) -> Self;
32
33    /// Creates a new fixed point number from an integer string. Precision tells us how many decimals
34    fn from_str(s: &str, precision: u64) -> Result<Self, Error>;
35
36    /// Returns the precision of the fixed point number
37    fn precision(&self) -> u64;
38
39    /// Performs the 'exp' approximation. First does the scaling of 'x' to \[0,1\]
40    /// and then calls the continued fraction approximation function.
41    fn exp(&self) -> Self;
42
43    /// Entry point for 'ln' approximation. First does the necessary scaling, and
44    /// then calls the continued fraction calculation. For any value outside the
45    /// domain, i.e., 'x in (-inf,0]', the function panics.
46    fn ln(&self) -> Self;
47
48    /// Entry point for 'pow' function. x^y = exp(y * ln x)
49    fn pow(&self, y: &Self) -> Self;
50
51    /// Entry point for bounded iterations for comparing two exp values.
52    fn exp_cmp(&self, max_n: u64, bound_self: i64, compare: &Self) -> ExpCmpOrdering;
53
54    /// Round to the nearest integer number
55    #[must_use]
56    fn round(&self) -> Self;
57
58    /// Round down to the nearest integer number
59    #[must_use]
60    fn floor(&self) -> Self;
61
62    /// Round up to the nearest integer number
63    #[must_use]
64    fn ceil(&self) -> Self;
65
66    /// Truncate to the nearest integer number
67    #[must_use]
68    fn trunc(&self) -> Self;
69}
70
71#[derive(Debug, Clone, PartialEq)]
72pub enum ExpOrdering {
73    GT,
74    LT,
75    UNKNOWN,
76}
77
78impl From<&str> for ExpOrdering {
79    fn from(s: &str) -> Self {
80        match s {
81            "GT" => ExpOrdering::GT,
82            "LT" => ExpOrdering::LT,
83            _ => ExpOrdering::UNKNOWN,
84        }
85    }
86}
87
88#[derive(Debug, Clone, PartialEq)]
89pub struct ExpCmpOrdering {
90    pub iterations: u64,
91    pub estimation: ExpOrdering,
92    pub approx: FixedDecimal,
93}
94
95#[cfg(test)]
96mod tests {
97    use super::*;
98    use dashu_base::Abs;
99    use proptest::prelude::Strategy;
100    use proptest::proptest;
101    use std::fs::File;
102    use std::io::BufRead;
103    use std::path::PathBuf;
104
105    #[test]
106    fn test_fixed_precision() {
107        let fp: FixedDecimal = FixedDecimal::new(34);
108        assert_eq!(fp.precision(), 34);
109        assert_eq!(fp.to_string(), "0.0000000000000000000000000000000000");
110    }
111
112    #[test]
113    fn test_fixed_precision_eq() {
114        let fp1: FixedDecimal = FixedDecimal::new(34);
115        let fp2: FixedDecimal = FixedDecimal::new(34);
116        assert_eq!(fp1, fp2);
117    }
118
119    #[test]
120    fn test_fixed_precision_from_str() {
121        let fp: FixedDecimal =
122            FixedDecimal::from_str("1234567890123456789012345678901234", 34).unwrap();
123        assert_eq!(fp.precision(), 34);
124        assert_eq!(fp.to_string(), "0.1234567890123456789012345678901234");
125
126        let fp: FixedDecimal =
127            FixedDecimal::from_str("-1234567890123456789012345678901234", 30).unwrap();
128        assert_eq!(fp.precision(), 30);
129        assert_eq!(fp.to_string(), "-1234.567890123456789012345678901234");
130
131        let fp: FixedDecimal =
132            FixedDecimal::from_str("-1234567890123456789012345678901234", 34).unwrap();
133        assert_eq!(fp.precision(), 34);
134        assert_eq!(fp.to_string(), "-0.1234567890123456789012345678901234");
135    }
136
137    #[test]
138    fn test_fixed_precision_exp() {
139        let fp: FixedDecimal = FixedDecimal::from(1u64);
140        assert_eq!(fp.to_string(), "1.0000000000000000000000000000000000");
141        let exp_fp = fp.exp();
142        assert_eq!(exp_fp.to_string(), "2.7182818284590452353602874043083282");
143    }
144
145    #[test]
146    fn test_fixed_precision_mul() {
147        let fp1: FixedDecimal =
148            FixedDecimal::from_str("52500000000000000000000000000000000", 34).unwrap();
149        let fp2: FixedDecimal =
150            FixedDecimal::from_str("43000000000000000000000000000000000", 34).unwrap();
151        let fp3 = &fp1 * &fp2;
152        assert_eq!(fp3.to_string(), "22.5750000000000000000000000000000000");
153        let fp4 = fp1 * fp2;
154        assert_eq!(fp4.to_string(), "22.5750000000000000000000000000000000");
155    }
156
157    #[test]
158    fn test_fixed_precision_div() {
159        let fp1: FixedDecimal = FixedDecimal::from_str("1", 34).unwrap();
160        let fp2: FixedDecimal = FixedDecimal::from_str("10", 34).unwrap();
161        let fp3 = &fp1 / &fp2;
162        assert_eq!(fp3.to_string(), "0.1000000000000000000000000000000000");
163        let fp4 = fp1 / fp2;
164        assert_eq!(fp4.to_string(), "0.1000000000000000000000000000000000");
165    }
166
167    #[test]
168    fn test_fixed_precision_sub() {
169        let fp1: FixedDecimal = FixedDecimal::from_str("1", 34).unwrap();
170        assert_eq!(fp1.to_string(), "0.0000000000000000000000000000000001");
171        let fp2: FixedDecimal = FixedDecimal::from_str("10", 34).unwrap();
172        assert_eq!(fp2.to_string(), "0.0000000000000000000000000000000010");
173        let fp3 = &fp1 - &fp2;
174        assert_eq!(fp3.to_string(), "-0.0000000000000000000000000000000009");
175        let fp4 = fp1 - fp2;
176        assert_eq!(fp4.to_string(), "-0.0000000000000000000000000000000009");
177    }
178
179    #[test]
180    fn test_fixed_precision_round() {
181        let fp1: FixedDecimal =
182            FixedDecimal::from_str("11234567890123456789012345678901234", 34).unwrap();
183        assert_eq!(
184            fp1.round().to_string(),
185            "1.0000000000000000000000000000000000"
186        );
187        let fp2: FixedDecimal =
188            FixedDecimal::from_str("14999999999999999999999999999999999", 34).unwrap();
189        assert_eq!(
190            fp2.round().to_string(),
191            "1.0000000000000000000000000000000000"
192        );
193        let fp3: FixedDecimal =
194            FixedDecimal::from_str("15000000000000000000000000000000000", 34).unwrap();
195        assert_eq!(
196            fp3.round().to_string(),
197            "2.0000000000000000000000000000000000"
198        );
199        let fp4: FixedDecimal = FixedDecimal::from_str("1500", 3).unwrap();
200        assert_eq!(fp4.round().to_string(), "2.000");
201        let fp5: FixedDecimal = FixedDecimal::from_str("1499", 3).unwrap();
202        assert_eq!(fp5.round().to_string(), "1.000");
203        let fp6: FixedDecimal =
204            FixedDecimal::from_str("-11234567890123456789012345678901234", 34).unwrap();
205        assert_eq!(
206            fp6.round().to_string(),
207            "-1.0000000000000000000000000000000000"
208        );
209        let fp2: FixedDecimal =
210            FixedDecimal::from_str("-14999999999999999999999999999999999", 34).unwrap();
211        assert_eq!(
212            fp2.round().to_string(),
213            "-1.0000000000000000000000000000000000"
214        );
215        let fp3: FixedDecimal =
216            FixedDecimal::from_str("-15000000000000000000000000000000000", 34).unwrap();
217        assert_eq!(
218            fp3.round().to_string(),
219            "-2.0000000000000000000000000000000000"
220        );
221        let fp4: FixedDecimal = FixedDecimal::from_str("-1500", 3).unwrap();
222        assert_eq!(fp4.round().to_string(), "-2.000");
223        let fp5: FixedDecimal = FixedDecimal::from_str("-1499", 3).unwrap();
224        assert_eq!(fp5.round().to_string(), "-1.000");
225        let fp6: FixedDecimal = FixedDecimal::from_str("1000", 3).unwrap();
226        assert_eq!(fp6.round().to_string(), "1.000");
227        let fp7: FixedDecimal = FixedDecimal::from_str("-1000", 3).unwrap();
228        assert_eq!(fp7.round().to_string(), "-1.000");
229    }
230
231    #[test]
232    fn test_fixed_precision_floor() {
233        let fp1: FixedDecimal =
234            FixedDecimal::from_str("11234567890123456789012345678901234", 34).unwrap();
235        assert_eq!(
236            fp1.floor().to_string(),
237            "1.0000000000000000000000000000000000"
238        );
239        let fp2: FixedDecimal =
240            FixedDecimal::from_str("14999999999999999999999999999999999", 34).unwrap();
241        assert_eq!(
242            fp2.floor().to_string(),
243            "1.0000000000000000000000000000000000"
244        );
245        let fp3: FixedDecimal =
246            FixedDecimal::from_str("15000000000000000000000000000000000", 34).unwrap();
247        assert_eq!(
248            fp3.floor().to_string(),
249            "1.0000000000000000000000000000000000"
250        );
251        let fp4: FixedDecimal = FixedDecimal::from_str("1500", 3).unwrap();
252        assert_eq!(fp4.floor().to_string(), "1.000");
253        let fp5: FixedDecimal = FixedDecimal::from_str("1499", 3).unwrap();
254        assert_eq!(fp5.floor().to_string(), "1.000");
255        let fp6: FixedDecimal =
256            FixedDecimal::from_str("-11234567890123456789012345678901234", 34).unwrap();
257        assert_eq!(
258            fp6.floor().to_string(),
259            "-2.0000000000000000000000000000000000"
260        );
261        let fp2: FixedDecimal =
262            FixedDecimal::from_str("-14999999999999999999999999999999999", 34).unwrap();
263        assert_eq!(
264            fp2.floor().to_string(),
265            "-2.0000000000000000000000000000000000"
266        );
267        let fp3: FixedDecimal =
268            FixedDecimal::from_str("-15000000000000000000000000000000000", 34).unwrap();
269        assert_eq!(
270            fp3.floor().to_string(),
271            "-2.0000000000000000000000000000000000"
272        );
273        let fp4: FixedDecimal = FixedDecimal::from_str("-1500", 3).unwrap();
274        assert_eq!(fp4.floor().to_string(), "-2.000");
275        let fp5: FixedDecimal = FixedDecimal::from_str("-1499", 3).unwrap();
276        assert_eq!(fp5.floor().to_string(), "-2.000");
277        let fp6: FixedDecimal = FixedDecimal::from_str("1000", 3).unwrap();
278        assert_eq!(fp6.floor().to_string(), "1.000");
279        let fp7: FixedDecimal = FixedDecimal::from_str("-1000", 3).unwrap();
280        assert_eq!(fp7.floor().to_string(), "-1.000");
281    }
282
283    #[test]
284    fn test_fixed_precision_ceil() {
285        let fp1: FixedDecimal =
286            FixedDecimal::from_str("11234567890123456789012345678901234", 34).unwrap();
287        assert_eq!(
288            fp1.ceil().to_string(),
289            "2.0000000000000000000000000000000000"
290        );
291        let fp2: FixedDecimal =
292            FixedDecimal::from_str("14999999999999999999999999999999999", 34).unwrap();
293        assert_eq!(
294            fp2.ceil().to_string(),
295            "2.0000000000000000000000000000000000"
296        );
297        let fp3: FixedDecimal =
298            FixedDecimal::from_str("15000000000000000000000000000000000", 34).unwrap();
299        assert_eq!(
300            fp3.ceil().to_string(),
301            "2.0000000000000000000000000000000000"
302        );
303        let fp4: FixedDecimal = FixedDecimal::from_str("1500", 3).unwrap();
304        assert_eq!(fp4.ceil().to_string(), "2.000");
305        let fp5: FixedDecimal = FixedDecimal::from_str("1499", 3).unwrap();
306        assert_eq!(fp5.ceil().to_string(), "2.000");
307        let fp6: FixedDecimal =
308            FixedDecimal::from_str("-11234567890123456789012345678901234", 34).unwrap();
309        assert_eq!(
310            fp6.ceil().to_string(),
311            "-1.0000000000000000000000000000000000"
312        );
313        let fp2: FixedDecimal =
314            FixedDecimal::from_str("-14999999999999999999999999999999999", 34).unwrap();
315        assert_eq!(
316            fp2.ceil().to_string(),
317            "-1.0000000000000000000000000000000000"
318        );
319        let fp3: FixedDecimal =
320            FixedDecimal::from_str("-15000000000000000000000000000000000", 34).unwrap();
321        assert_eq!(
322            fp3.ceil().to_string(),
323            "-1.0000000000000000000000000000000000"
324        );
325        let fp4: FixedDecimal = FixedDecimal::from_str("-1500", 3).unwrap();
326        assert_eq!(fp4.ceil().to_string(), "-1.000");
327        let fp5: FixedDecimal = FixedDecimal::from_str("-1499", 3).unwrap();
328        assert_eq!(fp5.ceil().to_string(), "-1.000");
329        let fp6: FixedDecimal = FixedDecimal::from_str("1000", 3).unwrap();
330        assert_eq!(fp6.ceil().to_string(), "1.000");
331        let fp7: FixedDecimal = FixedDecimal::from_str("-1000", 3).unwrap();
332        assert_eq!(fp7.ceil().to_string(), "-1.000");
333    }
334
335    #[test]
336    fn test_fixed_precision_trunc() {
337        let fp1: FixedDecimal =
338            FixedDecimal::from_str("11234567890123456789012345678901234", 34).unwrap();
339        assert_eq!(
340            fp1.trunc().to_string(),
341            "1.0000000000000000000000000000000000"
342        );
343        let fp2: FixedDecimal =
344            FixedDecimal::from_str("14999999999999999999999999999999999", 34).unwrap();
345        assert_eq!(
346            fp2.trunc().to_string(),
347            "1.0000000000000000000000000000000000"
348        );
349        let fp3: FixedDecimal =
350            FixedDecimal::from_str("15000000000000000000000000000000000", 34).unwrap();
351        assert_eq!(
352            fp3.trunc().to_string(),
353            "1.0000000000000000000000000000000000"
354        );
355        let fp4: FixedDecimal = FixedDecimal::from_str("1500", 3).unwrap();
356        assert_eq!(fp4.trunc().to_string(), "1.000");
357        let fp5: FixedDecimal = FixedDecimal::from_str("1499", 3).unwrap();
358        assert_eq!(fp5.trunc().to_string(), "1.000");
359        let fp6: FixedDecimal =
360            FixedDecimal::from_str("-11234567890123456789012345678901234", 34).unwrap();
361        assert_eq!(
362            fp6.trunc().to_string(),
363            "-1.0000000000000000000000000000000000"
364        );
365        let fp2: FixedDecimal =
366            FixedDecimal::from_str("-14999999999999999999999999999999999", 34).unwrap();
367        assert_eq!(
368            fp2.trunc().to_string(),
369            "-1.0000000000000000000000000000000000"
370        );
371        let fp3: FixedDecimal =
372            FixedDecimal::from_str("-15000000000000000000000000000000000", 34).unwrap();
373        assert_eq!(
374            fp3.trunc().to_string(),
375            "-1.0000000000000000000000000000000000"
376        );
377        let fp4: FixedDecimal = FixedDecimal::from_str("-1500", 3).unwrap();
378        assert_eq!(fp4.trunc().to_string(), "-1.000");
379        let fp5: FixedDecimal = FixedDecimal::from_str("-1499", 3).unwrap();
380        assert_eq!(fp5.trunc().to_string(), "-1.000");
381        let fp6: FixedDecimal = FixedDecimal::from_str("1000", 3).unwrap();
382        assert_eq!(fp6.trunc().to_string(), "1.000");
383        let fp7: FixedDecimal = FixedDecimal::from_str("-1000", 3).unwrap();
384        assert_eq!(fp7.trunc().to_string(), "-1.000");
385    }
386
387    #[test]
388    fn golden_tests() {
389        let mut data_path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
390        data_path.push("tests/data/golden_tests.txt");
391
392        // read each line of golden_tests.txt
393        let file = File::open(data_path).expect("golden_tests.txt: file not found");
394        let reader = std::io::BufReader::new(file);
395
396        // read each line of golden_tests_result.txt
397        let mut data_path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
398        data_path.push("tests/data/golden_tests_result.txt");
399        let file = File::open(data_path).expect("golden_tests_result.txt: file not found");
400        let result_reader = std::io::BufReader::new(file);
401
402        let one: FixedDecimal = FixedDecimal::from(1u64);
403        let ten: FixedDecimal = FixedDecimal::from(10u64);
404        let f: FixedDecimal = &one / &ten;
405        assert_eq!(f.to_string(), "0.1000000000000000000000000000000000");
406
407        for (test_line, result_line) in reader.lines().zip(result_reader.lines()) {
408            let test_line = test_line.expect("failed to read line");
409            // println!("test_line: {}", test_line);
410            let mut parts = test_line.split_whitespace();
411            let x = FixedDecimal::from_str(parts.next().unwrap(), DEFAULT_PRECISION)
412                .expect("failed to parse x");
413            let a = FixedDecimal::from_str(parts.next().unwrap(), DEFAULT_PRECISION)
414                .expect("failed to parse a");
415            let b = FixedDecimal::from_str(parts.next().unwrap(), DEFAULT_PRECISION)
416                .expect("failed to parse b");
417            let result_line = result_line.expect("failed to read line");
418            // println!("result_line: {}", result_line);
419            let mut result_parts = result_line.split_whitespace();
420            let expected_exp_x = result_parts.next().expect("expected_exp_x not found");
421            let expected_ln_a = result_parts.next().expect("expected_ln_a not found");
422            let expected_threshold_b = result_parts.next().expect("expected_threshold_b not found");
423            let expected_approx_exp = result_parts.next().expect("expected_approx_exp not found");
424            let expected_estimation =
425                ExpOrdering::from(result_parts.next().expect("expected_estimation not found"));
426            let expected_iterations = result_parts.next().expect("expected_iterations not found");
427
428            // calculate exp' x
429            let exp_x = x.exp();
430            assert_eq!(exp_x.to_string(), expected_exp_x);
431
432            // calculate ln' a, print -ln' a
433            let ln_a = a.ln();
434            assert_eq!((-ln_a).to_string(), expected_ln_a);
435
436            // calculate (1 - f) *** b
437            let c = &one - &f;
438            assert_eq!(c.to_string(), "0.9000000000000000000000000000000000");
439            let threshold_b = c.pow(&b);
440            assert_eq!(
441                (&one - &threshold_b).to_string(),
442                expected_threshold_b,
443                "(1 - f) *** b failed to match! - (1 - f)={}, b={}",
444                &c,
445                &b
446            );
447
448            // do Taylor approximation for
449            //  a < 1 - (1 - f) *** b <=> 1/(1-a) < exp(-b * ln' (1 - f))
450            // using Lagrange error term calculation
451            let c = &one - &f;
452            let temp = c.ln();
453            let alpha = &b * &temp;
454            let alpha = -alpha;
455            let q_ = &one - &a;
456            let q = &one / &q_;
457            let res = alpha.exp_cmp(1000, 3, &q);
458
459            // println!("alpha: {}", alpha);
460            // println!("q: {}", q);
461            // println!("res.approx: {}", res.approx);
462            // println!("res.estimation: {:?}", res.estimation);
463            // println!("res.iterations: {}", res.iterations);
464
465            // we compare 1/(1-p) < e^-(1-(1-f)^sigma)
466            let threshold = &one - &threshold_b;
467            if a < threshold && res.estimation != ExpOrdering::LT {
468                panic!(
469                    "wrong result should be leader {} should be more like {}",
470                    &temp, threshold
471                );
472            }
473
474            if a >= threshold && res.estimation != ExpOrdering::GT {
475                panic!(
476                    "wrong result should not be leader {} should be more like {}",
477                    &temp, threshold
478                );
479            }
480
481            assert_eq!(res.approx.to_string(), expected_approx_exp);
482            assert_eq!(res.estimation, expected_estimation);
483            assert_eq!(res.iterations.to_string(), expected_iterations);
484        }
485    }
486
487    #[test]
488    #[should_panic(expected = "ln of a value in (-inf,0] is undefined")]
489    fn ln_of_0_should_be_undefined() {
490        ZERO.ln();
491    }
492
493    #[test]
494    #[should_panic(expected = "ln of a value in (-inf,0] is undefined")]
495    fn ln_of_negative_should_be_undefined() {
496        MINUS_ONE.ln();
497    }
498
499    #[test]
500    fn pow_of_zero_to_any_positive_power_should_be_zero() {
501        proptest!(|(y in 1u64..=u64::MAX)| {
502            assert_eq!(ZERO.pow(&FixedDecimal::from(y)), *ZERO);
503        });
504    }
505
506    #[test]
507    #[should_panic(expected = "zero to a negative power is undefined")]
508    fn pow_of_zero_to_neg_power_should_be_undefined() {
509        let y = FixedDecimal::from(-1i64);
510        ZERO.pow(&y);
511    }
512
513    #[test]
514    fn pow_of_any_to_power_0_should_be_1() {
515        proptest!(|(x in i64::MIN..=i64::MAX)| {
516            assert_eq!(FixedDecimal::from(x).pow(&*ZERO), *ONE);
517        });
518    }
519
520    #[test]
521    fn pow_of_any_to_power_1_should_be_same() {
522        proptest!(|(x in i64::MIN..=i64::MAX)| {
523            assert_eq!(FixedDecimal::from(x).pow(&*ONE), FixedDecimal::from(x));
524        });
525    }
526
527    #[test]
528    fn pow_to_positive_times_pow_to_negative_should_be_1() {
529        let epsilon = FixedDecimal::from_str("1000000000000000000", 34).unwrap();
530        proptest!(|(x in (-5i64..=5i64).prop_filter("Exclude zero", |&x| x != 0), y in 1i64..=25i64)| {
531            let x = FixedDecimal::from(x);
532            let y = FixedDecimal::from(y);
533            let minus_y = -&y;
534            let x_to_y = x.pow(&y);
535            let x_to_minus_y = x.pow(&minus_y);
536            let result = &x_to_y * &x_to_minus_y;
537            let diff = (&result - &*ONE).abs();
538            // println!("x: {}, y: {}, x^y: {}, x^-y: {}, x^y * x^-y: {}, diff: {}", x, y, x_to_y, x_to_minus_y, result, diff);
539            assert!(diff <= epsilon);
540        });
541    }
542}
pallas_math/math.rs

pallas_math/
math.rs
