num_bigfloat/
lib.rs

1//! Increased precision floating point numbers implemented purely in Rust.
2//! Number has fixed-size mantissa and exponent, but increased precision compared to f32 or f64 values.
3//!
4//! Number characteristics:
5//!
6//! | Name                          | Value  |
7//! |:------------------------------|-------:|
8//! | Base                          |     10 |
9//! | Decimal positions in mantissa |     40 |
10//! | Exponent minimum value        |   -128 |
11//! | Exponent maximum value        |    127 |
12//!
13//! ## Examples
14//!
15//! ```
16//! use num_bigfloat::BigFloat;
17//! use num_bigfloat::ONE;
18//! use num_bigfloat::PI;
19//!
20//! // compute pi: pi = 6*arctan(1/sqrt(3))
21//! # #[cfg(feature = "std")] {
22//! let six: BigFloat = 6.0.into();
23//! let three: BigFloat = BigFloat::parse("3.0").unwrap();
24//! let pi = six * (ONE / three.sqrt()).atan();
25//! let epsilon = 1.0e-38.into();
26//!
27//! assert!((pi - PI).abs() < epsilon);
28//!
29//! println!("{}", pi);
30//! # }
31//! // output: 3.141592653589793238462643383279502884196e-39
32//! ```
33//!
34//! The same example using functions:
35//!
36//! ```
37//! use num_bigfloat::BigFloat;
38//! use num_bigfloat::ONE;
39//! use num_bigfloat::PI;
40//!
41//! // compute pi: pi = 6*arctan(1/sqrt(3))
42//! let six: BigFloat = BigFloat::from_u8(6);
43//! let three: BigFloat = BigFloat::from_u8(3);
44//! let pi = six.mul(&ONE.div(&three.sqrt()).atan());
45//! let epsilon = BigFloat::from_f64(1.0e-38);  // note: conversion from f64,f32 are not loss-less for `no_std`.
46//!
47//! assert!(pi.sub(&PI).abs().cmp(&epsilon).unwrap() < 0);
48//! ```
49//!
50//! ## no_std
51//!
52//! Library can be used without the standard Rust library. This can be achieved by turning off `std` feature.
53//!
54//! ## Other features
55//!
56//! The library depends on [rand](https://crates.io/crates/rand) which is used by `BigFloat::random_normal`. This dependecy can be excluded by turning off the `rand` feature.
57//! `rand` feature requires `std` feature.
58//!
59//! The library depends on [serde](https://crates.io/crates/serde) which is also optional and can be eliminated by turning off the `serde` feature.
60//!
61//! In addition, the library implements [num-traits](https://crates.io/crates/num_traits). Dependency on `num-traits` can be excluded by turning off `num-traits` feature.
62
63#![deny(missing_docs)]
64#![deny(clippy::suspicious)]
65#![cfg_attr(not(feature = "std"), no_std)]
66
67mod defs;
68mod ext;
69mod inc;
70mod ops;
71mod parser;
72mod util;
73
74#[cfg(feature = "num-traits")]
75mod num_traits;
76
77#[cfg(feature = "serde")]
78pub mod serde;
79
80pub use crate::defs::Error;
81pub use crate::defs::RoundingMode;
82pub use crate::ext::BigFloat;
83pub use crate::ext::E;
84pub use crate::ext::EPSILON;
85pub use crate::ext::FRAC_1_PI;
86pub use crate::ext::FRAC_1_SQRT_2;
87pub use crate::ext::FRAC_2_PI;
88pub use crate::ext::FRAC_2_SQRT_PI;
89pub use crate::ext::FRAC_PI_3;
90pub use crate::ext::FRAC_PI_4;
91pub use crate::ext::FRAC_PI_6;
92pub use crate::ext::FRAC_PI_8;
93pub use crate::ext::HALF_PI;
94pub use crate::ext::INF_NEG;
95pub use crate::ext::INF_POS;
96pub use crate::ext::LN_10;
97pub use crate::ext::LN_2;
98pub use crate::ext::LOG10_E;
99pub use crate::ext::LOG2_E;
100pub use crate::ext::MAX;
101pub use crate::ext::MAX_EXP;
102pub use crate::ext::MIN;
103pub use crate::ext::MIN_EXP;
104pub use crate::ext::MIN_POSITIVE;
105pub use crate::ext::MIN_POSITIVE_NORMAL;
106pub use crate::ext::NAN;
107pub use crate::ext::ONE;
108pub use crate::ext::PI;
109pub use crate::ext::RADIX;
110pub use crate::ext::SQRT_2;
111pub use crate::ext::TWO;
112pub use crate::ext::ZERO;
113
114#[cfg(test)]
115mod tests {
116
117    use crate::defs::{
118        DECIMAL_BASE, DECIMAL_MAX_EXPONENT, DECIMAL_MIN_EXPONENT, DECIMAL_PARTS, DECIMAL_POSITIONS,
119        DECIMAL_SIGN_NEG, DECIMAL_SIGN_POS,
120    };
121    use crate::{
122        BigFloat, HALF_PI, INF_NEG, INF_POS, MAX, MIN, MIN_POSITIVE, NAN, ONE, PI, TWO, ZERO,
123    };
124    use rand::random;
125
126    #[test]
127    fn test_lib_basic() {
128        let mut d1;
129        let mut d2;
130        let mut d3;
131        let mut ref_num;
132
133        // creation & deconstruction
134
135        // regular buf
136        let bytes1: [u8; 20] =
137            [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 112, 13, 14, 15, 16, 17, 18, 19, 20];
138        let expected1: [u8; 30] = [
139            1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
140            0,
141        ];
142        let exp1 = 123;
143        let d4 = BigFloat::from_bytes(&bytes1, 1, exp1);
144
145        let mut mantissa_buf1 = [0; 30];
146        d4.get_mantissa_bytes(&mut mantissa_buf1);
147        assert!(mantissa_buf1 == expected1);
148        assert!(d4.get_mantissa_len() == bytes1.len());
149        assert!(d4.get_sign() == 1);
150        assert!(d4.get_exponent() == exp1);
151
152        // too long buf
153        let bytes2: [u8; 45] = [
154            1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 112, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5,
155            6, 7, 8, 9, 10, 11, 112, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 3, 4, 5,
156        ];
157        let expected2: [u8; 42] = [
158            1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
159            0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 0,
160        ];
161        let exp2 = -128;
162        let d4 = BigFloat::from_bytes(&bytes2, -2, exp2);
163
164        let mut mantissa_buf2 = [0; 42];
165        d4.get_mantissa_bytes(&mut mantissa_buf2);
166        assert!(mantissa_buf2 == expected2);
167        assert!(d4.get_mantissa_len() == 40);
168        assert!(d4.get_sign() == -1);
169        assert!(d4.get_exponent() == exp2);
170
171        // conversions
172
173        // inf
174        assert!(BigFloat::from_f64(f64::INFINITY).cmp(&INF_POS) == Some(0));
175        assert!(BigFloat::from_f64(f64::NEG_INFINITY).cmp(&INF_NEG) == Some(0));
176
177        // nan
178        assert!(BigFloat::from_f64(f64::NAN).is_nan());
179
180        // 0.0
181        assert!(BigFloat::from_f64(0.0).to_f64() == 0.0);
182
183        // conversions
184        for _ in 0..10000 {
185            let f: f64 = random_f64_exp(50, 25);
186            if f.is_finite() && f != 0.0 {
187                d1 = BigFloat::from_f64(f);
188                assert!(d1.to_f64() == f);
189                if (f as f32).is_finite() && (f as f32) != 0.0 {
190                    d1 = BigFloat::from_f32(f as f32);
191                    assert!(d1.to_f32() == f as f32);
192                }
193            }
194        }
195
196        // 0 * 0
197        d1 = BigFloat::new();
198        d2 = BigFloat::new();
199        ref_num = BigFloat::new();
200        d3 = d1.mul(&d2);
201        assert!(d3.cmp(&ref_num) == Some(0));
202
203        // 0.99 * 0
204        d1 = BigFloat::from_f64(0.99);
205        d3 = d1.mul(&d2);
206        assert!(d3.cmp(&ref_num) == Some(0));
207
208        // 0 * 12349999
209        d1 = BigFloat::new();
210        d2 = BigFloat::from_u32(12349999);
211        d3 = d1.mul(&d2);
212        assert!(d3.cmp(&ref_num) == Some(0));
213
214        // 1 * 1
215        d1 = BigFloat::from_u8(1);
216        d2 = BigFloat::from_u8(1);
217        d3 = d1.mul(&d2);
218        assert!(d3.cmp(&d1) == Some(0));
219
220        // 1 * -1
221        d1 = BigFloat::from_u8(1);
222        d2 = BigFloat::from_u8(1).inv_sign();
223        d3 = d1.mul(&d2);
224        assert!(d3.cmp(&d2) == Some(0));
225
226        // -1 * 1
227        d3 = d2.mul(&d1);
228        assert!(d3.cmp(&d2) == Some(0));
229
230        // -1 * -1
231        d1 = d1.inv_sign();
232        d3 = d1.mul(&d2);
233        ref_num = BigFloat::from_u8(1);
234        assert!(d3.cmp(&ref_num) == Some(0));
235
236        // 0 / 0
237        d1 = BigFloat::new();
238        d2 = BigFloat::new();
239        assert!(d1.div(&d2).is_nan());
240
241        // d2 / 0
242        d2 = BigFloat::from_u8(123);
243        assert!(d2.div(&d1).is_inf_pos());
244
245        // 0 / d2
246        d3 = d1.div(&d2);
247        ref_num = BigFloat::new();
248        assert!(d3.cmp(&ref_num) == Some(0));
249
250        // 0 / -d2
251        d2 = d2.inv_sign();
252        d3 = d1.div(&d2);
253        assert!(d3.cmp(&ref_num) == Some(0));
254    }
255
256    #[test]
257    fn test_lib_add_sub() {
258        let mut d1;
259        let mut d2;
260        let mut ref_num;
261
262        // add & sub
263        for _ in 0..10000 {
264            // avoid subnormal numbers
265            d1 = random_normal_float(4, 30);
266            d2 = random_normal_float(4, 34);
267            let n1 = d1.add(&d2);
268            let n2 = n1.sub(&d2);
269            assert_small(&n2.sub(&d1), 2, &d1);
270        }
271
272        // mul & div
273        for _ in 0..10000 {
274            // avoid subnormal numbers
275            d1 = random_normal_float(40, 40);
276            d2 = random_normal_float(40, 40);
277            if d2.cmp(&ZERO).unwrap() != 0 {
278                let n1 = d1.div(&d2);
279                let n2 = n1.mul(&d2);
280                assert_small(&n2.sub(&d1), 2, &d1);
281            }
282        }
283
284        // subnormal numbers
285        d1 = MIN_POSITIVE;
286        d2 = MIN_POSITIVE;
287        ref_num = MIN_POSITIVE.mul(&TWO);
288
289        // min_positive + min_positive = 2*min_positive
290        assert!(d1.add(&d2).cmp(&ref_num) == Some(0));
291        assert!(d1.add(&d2).cmp(&d1).unwrap() > 0);
292        assert!(d1.cmp(&d1.add(&d2)).unwrap() < 0);
293
294        // min_positive - min_positive = 0
295        ref_num = BigFloat::new();
296        assert!(d1.sub(&d2).cmp(&ref_num) == Some(0));
297
298        // 1 * min_positive = min_positive
299        assert!(ONE.mul(&d2).cmp(&d2) == Some(0));
300
301        // min_positive / 1 = min_positive
302        assert!(d2.div(&ONE).cmp(&d2) == Some(0));
303
304        // min_positive / 1 = min_positive
305        assert!(d2.div(&ONE).cmp(&d2) == Some(0));
306
307        // normal -> subnormal -> normal
308        d1 = ONE;
309        d1.set_exponent(DECIMAL_MIN_EXPONENT);
310        d2 = MIN_POSITIVE;
311        assert!(!d1.is_subnormal());
312        assert!(d1.sub(&d2).cmp(&d1).unwrap() < 0);
313        assert!(d1.cmp(&d1.sub(&d2)).unwrap() > 0);
314        d1 = d1.sub(&d2);
315        assert!(d1.is_subnormal());
316        d1 = d1.add(&d2);
317        assert!(!d1.is_subnormal());
318
319        // overflow
320        d1 = ONE;
321        d1.set_exponent(DECIMAL_MAX_EXPONENT - (DECIMAL_POSITIONS as i8 - 1));
322        assert!(MAX.add(&d1).is_inf_pos());
323        assert!(MIN.sub(&d1).is_inf_neg());
324        assert!(MAX.mul(&MAX).is_inf_pos());
325        d1 = ONE;
326        d1.set_exponent(DECIMAL_MIN_EXPONENT);
327        assert!(MAX.div(&d1).is_inf_pos());
328    }
329
330    #[test]
331    fn test_lib_fract_int_abs() {
332        let mut d1;
333
334        // fract & int
335        let f1 = 12345.6789;
336        d1 = BigFloat::from_f64(f1);
337        assert!((d1.frac().to_f64() - f1.fract()).abs() < 10000.0 * f64::EPSILON);
338        assert!((d1.int().to_f64() - (f1 as u64) as f64).abs() < 10000.0 * f64::EPSILON);
339
340        let f1 = -0.006789;
341        d1 = BigFloat::from_f64(f1);
342        assert!(d1.frac().cmp(&d1) == Some(0));
343        assert!(d1.int().is_zero());
344
345        d1 = BigFloat::from_bytes(&[2, 2, 2, 2, 2, 0, 0, 0], DECIMAL_SIGN_POS, -2);
346        assert!(d1.frac().is_zero());
347        assert!(d1.int().cmp(&d1) == Some(0));
348
349        assert!(MIN_POSITIVE.frac().cmp(&MIN_POSITIVE) == Some(0));
350        assert!(MIN_POSITIVE.int().is_zero());
351
352        d1 = BigFloat::new();
353        assert!(d1.frac().is_zero());
354        assert!(d1.int().is_zero());
355
356        // ceil & floor
357        d1 = BigFloat::from_f64(12.3);
358        assert!(d1.floor().to_f64() == 12.0);
359        assert!(d1.ceil().to_f64() == 13.0);
360        d1 = BigFloat::from_f64(12.0);
361        assert!(d1.floor().to_f64() == 12.0);
362        assert!(d1.ceil().to_f64() == 12.0);
363
364        d1 = BigFloat::from_f64(-12.3);
365        assert!(d1.floor().to_f64() == -13.0);
366        assert!(d1.ceil().to_f64() == -12.0);
367        d1 = BigFloat::from_f64(-12.0);
368        assert!(d1.floor().to_f64() == -12.0);
369        assert!(d1.ceil().to_f64() == -12.0);
370
371        // abs
372        d1 = BigFloat::from_f64(12.3);
373        assert!(d1.abs().to_f64() == 12.3);
374        d1 = BigFloat::from_f64(-12.3);
375        assert!(d1.abs().to_f64() == 12.3);
376    }
377
378    #[test]
379    fn test_lib_sqrt() {
380        // sqrt
381        for _ in 0..10000 {
382            let num = random_normal_float(256, 128).abs();
383            let sq = num.sqrt();
384            let ret = sq.mul(&sq);
385            assert_small(&num.sub(&ret), 2, &num);
386        }
387
388        // sqrt of max
389        let n = MAX.sqrt();
390        assert!(n.cmp(&ZERO).unwrap() > 0 && MAX.cmp(&n).unwrap() > 0);
391        let n = n.mul(&n);
392        assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 2, &MAX);
393
394        // sqrt of zero
395        let n = ZERO.sqrt();
396        assert!(n.cmp(&ZERO).unwrap() == 0);
397
398        // sqrt of min positive
399        let n = MIN_POSITIVE.sqrt();
400        assert!(n.cmp(&ZERO).unwrap() > 0 && MIN_POSITIVE.cmp(&n).unwrap() < 0);
401        let n = n.mul(&n);
402        assert_small_if_not(n.is_zero(), &MIN_POSITIVE.sub(&n), 2, &MIN_POSITIVE);
403
404        // sqrt of negative
405        let n = MIN_POSITIVE.inv_sign().sqrt();
406        assert!(n.is_nan());
407    }
408
409    #[test]
410    fn test_lib_cbrt() {
411        // cbrt
412        for _ in 0..10000 {
413            let num = random_normal_float(256, 128);
414            let sq = num.cbrt();
415            let ret = sq.mul(&sq).mul(&sq);
416            assert_small(&num.sub(&ret), 2, &num);
417        }
418
419        // cbrt of max
420        let n = MAX.cbrt();
421        assert!(n.cmp(&ZERO).unwrap() > 0 && MAX.cmp(&n).unwrap() > 0);
422        let n = n.mul(&n).mul(&n);
423        assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 2, &MAX);
424
425        // cbrt of zero
426        let n = ZERO.cbrt();
427        assert!(n.cmp(&ZERO).unwrap() == 0);
428
429        // cbrt of min positive
430        let n = MIN_POSITIVE.cbrt();
431        assert!(n.cmp(&ZERO).unwrap() > 0 && MIN_POSITIVE.cmp(&n).unwrap() < 0);
432        let n = n.mul(&n).mul(&n);
433        assert_small_if_not(n.is_zero(), &MIN_POSITIVE.sub(&n), 2, &MIN_POSITIVE);
434
435        // cbrt of negative MAX
436        let n = MAX.inv_sign().cbrt();
437        assert!(n.cmp(&ZERO).unwrap() < 0 && MAX.inv_sign().cmp(&n).unwrap() < 0);
438        let n = n.mul(&n).mul(&n);
439        assert_small_if_not(n.is_inf_neg(), &MAX.inv_sign().sub(&n), 2, &MAX.inv_sign());
440
441        // cbrt of negative MIN_POSITIVE
442        let n = MIN_POSITIVE.inv_sign().cbrt();
443        assert!(n.cmp(&ZERO).unwrap() < 0 && MIN_POSITIVE.inv_sign().cmp(&n).unwrap() > 0);
444        let n = n.mul(&n).mul(&n);
445        assert_small_if_not(
446            n.is_zero(),
447            &MIN_POSITIVE.inv_sign().sub(&n),
448            2,
449            &MIN_POSITIVE.inv_sign(),
450        );
451    }
452
453    #[test]
454    fn test_lib_pow() {
455        // pow
456        for _ in 0..10000 {
457            let a = random_normal_float(4, 40);
458            let n = random_normal_float(4, 40);
459            let inv = ONE.div(&n);
460            let p = a.pow(&n);
461            if !p.is_inf() && p.get_mantissa_len() >= DECIMAL_POSITIONS - 1 {
462                let ret = p.pow(&inv);
463                assert_small(&a.sub(&ret), 3, &a);
464            }
465        }
466
467        // one^any = one
468        assert!(ONE.pow(&MIN_POSITIVE).cmp(&ONE).unwrap() == 0);
469        assert!(ONE.pow(&MIN_POSITIVE.inv_sign()).cmp(&ONE).unwrap() == 0);
470        assert!(ONE.pow(&MAX).cmp(&ONE).unwrap() == 0);
471        assert!(ONE.pow(&MIN).cmp(&ONE).unwrap() == 0);
472        assert!(ONE.pow(&TWO).cmp(&ONE).unwrap() == 0);
473        assert!(ONE.pow(&ONE.div(&TWO)).cmp(&ONE).unwrap() == 0);
474        assert!(ONE.pow(&ZERO).cmp(&ONE).unwrap() == 0);
475
476        // (> 1)^n
477        assert!(TWO.pow(&MIN_POSITIVE).cmp(&ONE).unwrap() == 0);
478        assert!(TWO.pow(&MIN_POSITIVE.inv_sign()).cmp(&ONE).unwrap() == 0);
479        assert!(TWO.pow(&MAX).is_inf_pos());
480        assert!(TWO.pow(&MIN).is_zero());
481        assert!(TWO.pow(&ONE).cmp(&TWO).unwrap() == 0);
482        assert!(TWO.pow(&ZERO).cmp(&ONE).unwrap() == 0);
483
484        // (< 1)^n
485        let n = ONE.div(&TWO);
486        assert!(n.pow(&MIN_POSITIVE).cmp(&ONE).unwrap() == 0);
487        assert!(n.pow(&MIN_POSITIVE.inv_sign()).cmp(&ONE).unwrap() == 0);
488        assert!(n.pow(&MAX).is_zero());
489        assert!(n.pow(&MIN).is_inf_pos());
490        assert!(n.pow(&ONE).cmp(&n).unwrap() == 0);
491        assert!(n.pow(&ZERO).cmp(&ONE).unwrap() == 0);
492
493        // 0^n
494        assert!(ZERO.pow(&MIN_POSITIVE).is_zero());
495        assert!(ZERO.pow(&MIN_POSITIVE.inv_sign()).is_inf_pos());
496        assert!(ZERO.pow(&MAX).is_zero());
497        assert!(ZERO.pow(&MIN).is_inf_pos());
498        assert!(ZERO.pow(&ONE).is_zero());
499        assert!(ZERO.pow(&ZERO).is_nan());
500
501        // -(> 1)^n
502        let n = TWO.inv_sign();
503        assert!(n.pow(&MIN_POSITIVE).cmp(&ONE.inv_sign()).unwrap() == 0);
504        assert!(
505            n.pow(&MIN_POSITIVE.inv_sign())
506                .cmp(&ONE.inv_sign())
507                .unwrap()
508                == 0
509        );
510        assert!(n.pow(&MAX).is_inf_neg());
511        assert!(n.pow(&MIN).is_zero());
512        assert!(n.pow(&ONE).cmp(&TWO.inv_sign()).unwrap() == 0);
513        assert!(n.pow(&TWO).is_positive());
514        assert!(n.pow(&ZERO).cmp(&ONE).unwrap() == 0);
515
516        // -(< 1)^n
517        let n = ONE.div(&TWO).inv_sign();
518        assert!(n.pow(&MIN_POSITIVE).cmp(&ONE.inv_sign()).unwrap() == 0);
519        assert!(
520            n.pow(&MIN_POSITIVE.inv_sign())
521                .cmp(&ONE.inv_sign())
522                .unwrap()
523                == 0
524        );
525        assert!(n.pow(&MAX).is_zero());
526        assert!(n.pow(&MIN).is_inf_neg());
527        assert!(n.pow(&ONE).cmp(&n).unwrap() == 0);
528        assert!(n.pow(&TWO).is_positive());
529        assert!(n.pow(&ZERO).cmp(&ONE).unwrap() == 0);
530
531        // MAX^n
532        assert!(MAX.pow(&MIN_POSITIVE).cmp(&ONE).unwrap() == 0);
533        assert!(MAX.pow(&MIN_POSITIVE.inv_sign()).cmp(&ONE).unwrap() == 0);
534        assert!(MAX.pow(&MAX).is_inf_pos());
535        assert!(MAX.pow(&MIN).is_zero());
536        assert!(MAX.pow(&ONE).cmp(&MAX).unwrap() == 0);
537        assert!(MAX.pow(&ZERO).cmp(&ONE).unwrap() == 0);
538
539        // MIN^n
540        assert!(MIN.pow(&MIN_POSITIVE).cmp(&ONE.inv_sign()).unwrap() == 0);
541        assert!(
542            MIN.pow(&MIN_POSITIVE.inv_sign())
543                .cmp(&ONE.inv_sign())
544                .unwrap()
545                == 0
546        );
547        assert!(MIN.pow(&MAX).is_inf_neg());
548        assert!(MIN.pow(&MIN).is_zero());
549        assert!(MIN.pow(&ONE).cmp(&MIN).unwrap() == 0);
550        assert!(MIN.pow(&ZERO).cmp(&ONE).unwrap() == 0);
551
552        // MIN_POSITIVE^n
553        assert!(MIN_POSITIVE.pow(&MIN_POSITIVE).cmp(&ONE).unwrap() == 0);
554        assert!(
555            MIN_POSITIVE
556                .pow(&MIN_POSITIVE.inv_sign())
557                .cmp(&ONE)
558                .unwrap()
559                == 0
560        );
561        assert!(MIN_POSITIVE.pow(&MAX).is_zero());
562        assert!(MIN_POSITIVE.pow(&MIN).is_inf_pos());
563        assert!(MIN_POSITIVE.pow(&ONE).cmp(&MIN_POSITIVE).unwrap() == 0);
564        assert!(MIN_POSITIVE.pow(&ZERO).cmp(&ONE).unwrap() == 0);
565
566        // (-MIN_POSITIVE)^n
567        let n = MIN_POSITIVE.inv_sign();
568        assert!(n.pow(&MIN_POSITIVE).cmp(&ONE.inv_sign()).unwrap() == 0);
569        assert!(
570            n.pow(&MIN_POSITIVE.inv_sign())
571                .cmp(&ONE.inv_sign())
572                .unwrap()
573                == 0
574        );
575        assert!(n.pow(&MAX).is_zero());
576        assert!(n.pow(&MIN).is_inf_neg());
577        assert!(n.pow(&ONE).cmp(&n).unwrap() == 0);
578        assert!(n.pow(&TWO).is_positive());
579        assert!(n.pow(&ZERO).cmp(&ONE).unwrap() == 0);
580    }
581
582    #[test]
583    fn test_lib_ln_log2_log10() {
584        // ln
585        let ten = BigFloat::from_u8(10);
586        for _ in 0..10000 {
587            let num = random_normal_float(256, 127).abs();
588            if num.is_zero() {
589                assert!(num.ln().is_nan());
590                assert!(num.log2().is_nan());
591                assert!(num.log10().is_nan());
592            } else {
593                let l = num.ln();
594                let e = l.exp();
595                assert_small(&num.sub(&e), 4, &num);
596
597                let l = num.log2();
598                let e = TWO.pow(&l);
599                assert_small(&num.sub(&e), 4, &num);
600
601                let l = num.log10();
602                let e = ten.pow(&l);
603                assert_small(&num.sub(&e), 5, &num);
604            }
605        }
606
607        let ops = [BigFloat::ln, BigFloat::log2, BigFloat::log10];
608        let bases = [crate::E, TWO, ten];
609        for i in 0..ops.len() {
610            let op = ops[i];
611            let base = bases[i];
612
613            // crossing x axis at x = 1
614            let n = op(&ONE);
615            assert_small_if_not(n.is_zero(), &ZERO.sub(&n), 2, &n);
616
617            // ln of max
618            let n = op(&MAX);
619            assert!(n.cmp(&ZERO).unwrap() > 0 && MAX.cmp(&n).unwrap() > 0);
620            let n = base.pow(&n);
621            assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 4, &MAX);
622
623            // ln of min positive
624            let n = op(&MIN_POSITIVE);
625            assert!(n.cmp(&ZERO).unwrap() < 0 && MIN_POSITIVE.cmp(&n.abs()).unwrap() < 0);
626            let n = base.pow(&n);
627            assert_small_if_not(n.is_inf_neg(), &MIN_POSITIVE.sub(&n), 2, &n);
628
629            // ln of negative
630            let n = op(&MIN_POSITIVE.inv_sign());
631            assert!(n.is_nan());
632
633            // ln of zero
634            let n = op(&ZERO.inv_sign());
635            assert!(n.is_nan());
636        }
637    }
638
639    #[test]
640    fn test_lib_log() {
641        // log
642        for _ in 0..10000 {
643            let base = random_normal_float(256, 127).abs();
644            let num = random_normal_float(256, 127).abs();
645            if num.is_zero() || base.is_zero() {
646                assert!(num.log(&base).is_nan());
647            } else {
648                let l = num.log(&base);
649                let e = base.pow(&l);
650                assert_small(&num.sub(&e), 5, &num);
651            }
652        }
653
654        // log crossing x axis at x = 1
655        let base1 = BigFloat::from_f64(5.4321);
656        let base2 = BigFloat::from_f64(0.4321);
657        let n = ONE.log(&base1);
658        assert_small_if_not(n.is_zero(), &ZERO.sub(&n), 2, &n);
659        let n = ONE.log(&base2);
660        assert_small_if_not(n.is_zero(), &ZERO.sub(&n), 2, &n);
661
662        // log of max
663        let n = MAX.log(&base1);
664        assert!(n.cmp(&ZERO).unwrap() > 0 && MAX.cmp(&n).unwrap() > 0);
665        let n = base1.pow(&n);
666        assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 4, &MAX);
667        let n = MAX.log(&base2);
668        assert!(n.cmp(&ZERO).unwrap() < 0 && MAX.cmp(&n.inv_sign()).unwrap() > 0);
669        let n = base2.pow(&n);
670        assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 2, &MAX);
671
672        // log of min positive
673        let n = MIN_POSITIVE.log(&base1);
674        assert!(n.cmp(&ZERO).unwrap() < 0 && MIN_POSITIVE.cmp(&n.abs()).unwrap() < 0);
675        let n = base1.pow(&n);
676        assert_small_if_not(n.is_inf_neg(), &MIN_POSITIVE.sub(&n), 2, &n);
677        let n = MIN_POSITIVE.log(&base2);
678        assert!(n.cmp(&ZERO).unwrap() > 0 && MIN_POSITIVE.cmp(&n.abs()).unwrap() < 0);
679        let n = base2.pow(&n);
680        assert_small_if_not(n.is_inf_neg(), &MIN_POSITIVE.sub(&n), 2, &n);
681
682        // log of negative
683        let n = MIN_POSITIVE.inv_sign().log(&base1);
684        assert!(n.is_nan());
685        let n = MIN_POSITIVE.inv_sign().log(&base2);
686        assert!(n.is_nan());
687
688        // negative base
689        assert!(TWO.log(&base1.inv_sign()).is_nan());
690
691        // log of zero
692        let n = ZERO.inv_sign().log(&base1);
693        assert!(n.is_nan());
694        let n = ZERO.inv_sign().log(&base2);
695        assert!(n.is_nan());
696
697        // zero base
698        assert!(TWO.log(&ZERO).is_nan());
699    }
700
701    #[test]
702    fn test_lib_sin_asin() {
703        let eps = BigFloat::from_f64(1.0e-7);
704        // sin, asin
705        for _ in 0..10000 {
706            let num = random_normal_float(3, 40);
707            let s = num.sin();
708            let a = s.asin();
709            assert!(HALF_PI.cmp(&a).unwrap() >= 0 && HALF_PI.inv_sign().cmp(&a).unwrap() <= 0);
710            if ONE.sub(&s.abs()).cmp(&eps).unwrap() >= 0 {
711                // avoid values of sin close to 1;
712                // although computation of sin and asin is precise
713                // this test does not work for them.
714                if num.abs().cmp(&HALF_PI).unwrap() <= 0 {
715                    assert_small(&num.sub(&a), 4, &num);
716                } else {
717                    let mut sub1 = num.add(&a).abs().div(&PI).frac();
718                    let mut sub2 = num.sub(&a).abs().div(&PI).frac();
719                    if sub1.get_mantissa_len() > 3 {
720                        sub1 = ONE.sub(&sub1);
721                    }
722                    if sub2.get_mantissa_len() > 3 {
723                        sub2 = ONE.sub(&sub2);
724                    }
725                    assert_small_any(&sub1, &sub2, 3, &num);
726                }
727            }
728        }
729
730        // x axis crossing points: 0, PI, 2*PI
731        let n = ZERO.sin();
732        assert!(n.is_zero() || ZERO.sub(&n).get_mantissa_len() < 2);
733
734        let n = PI.sin();
735        assert!(n.is_zero() || ZERO.sub(&n).get_exponent() < -(DECIMAL_POSITIONS as i8) - 38);
736
737        let n = PI.add(&PI).sin();
738        assert!(n.is_zero() || ZERO.sub(&n).get_exponent() < -(DECIMAL_POSITIONS as i8) - 38);
739
740        // sin near 0
741        let n = MIN_POSITIVE.sin();
742        assert!(MIN_POSITIVE.sub(&n).get_mantissa_len() < 2);
743        let n = MIN_POSITIVE.inv_sign().sin();
744        assert!(MIN_POSITIVE.inv_sign().sub(&n).get_mantissa_len() < 2);
745
746        // sin extremums PI/2, 3*PI/2
747        let eps = BigFloat::from_f64(1.0e-39);
748        let exp_err = -(DECIMAL_POSITIONS as i8) - 18;
749        test_extremum(BigFloat::sin, &HALF_PI, &ONE, 3, 2, 2, &eps, exp_err);
750        test_extremum(
751            BigFloat::sin,
752            &HALF_PI.add(&PI),
753            &ONE.inv_sign(),
754            3,
755            2,
756            2,
757            &eps,
758            exp_err,
759        );
760
761        // asin extremums: 1, -1
762        test_extremum(BigFloat::asin, &ONE, &HALF_PI, 2, 2, 2, &eps, exp_err);
763        test_extremum(
764            BigFloat::asin,
765            &ONE.inv_sign(),
766            &HALF_PI.inv_sign(),
767            1,
768            2,
769            2,
770            &eps,
771            exp_err,
772        );
773
774        // asin near 0
775        let n = MIN_POSITIVE.asin();
776        assert!(MIN_POSITIVE.sub(&n).get_mantissa_len() < 2);
777        let n = MIN_POSITIVE.inv_sign().asin();
778        assert!(MIN_POSITIVE.inv_sign().sub(&n).get_mantissa_len() < 2);
779
780        // asin resulting to NAN
781        let n = TWO.asin();
782        assert!(n.is_nan());
783        let n = TWO.inv_sign().asin();
784        assert!(n.is_nan());
785    }
786
787    #[test]
788    fn test_lib_cos_acos() {
789        let eps = BigFloat::from_f64(1.0e-8);
790        for _ in 0..10000 {
791            let num = random_normal_float(3, 40);
792            let c = num.cos();
793            let a = c.acos();
794            assert!(PI.cmp(&a).unwrap() >= 0 && ZERO.inv_sign().cmp(&a).unwrap() <= 0);
795            if ONE.sub(&c.abs()).cmp(&eps).unwrap() >= 0 {
796                // avoid values of cos close to 1;
797                // although computation of cos and acos is precise
798                // this test does not work for them.
799                if num.abs().cmp(&PI).unwrap() <= 0 {
800                    assert_small(&num.abs().sub(&a), 4, &num);
801                } else {
802                    let mut sub1 = num.add(&a).abs().div(&PI).frac();
803                    let mut sub2 = num.sub(&a).abs().div(&PI).frac();
804                    if sub1.get_mantissa_len() > 3 {
805                        sub1 = ONE.sub(&sub1);
806                    }
807                    if sub2.get_mantissa_len() > 3 {
808                        sub2 = ONE.sub(&sub2);
809                    }
810                    assert_small_any(&sub1, &sub2, 3, &num);
811                }
812            }
813        }
814
815        // x axis crossing points: PI/2, 3*PI/2
816        let n = HALF_PI.cos();
817        assert!(n.is_zero() || ZERO.sub(&n).get_exponent() < -(DECIMAL_POSITIONS as i8) - 38);
818
819        let n = HALF_PI.add(&PI).cos();
820        assert!(n.is_zero() || ZERO.sub(&n).get_exponent() < -(DECIMAL_POSITIONS as i8) - 38);
821
822        // cos extremums: 0, PI
823        let eps = BigFloat::from_f64(1.0e-39);
824        let exp_err = -(DECIMAL_POSITIONS as i8) - 18;
825        test_extremum(BigFloat::cos, &ZERO, &ONE, 3, 2, 2, &eps, exp_err);
826        test_extremum(BigFloat::cos, &PI, &ONE.inv_sign(), 3, 2, 2, &eps, exp_err);
827
828        // acos extremums: 1, -1
829        test_extremum(BigFloat::acos, &ONE, &ZERO, 2, 2, 2, &eps, exp_err);
830        test_extremum(BigFloat::acos, &ONE.inv_sign(), &PI, 1, 2, 2, &eps, exp_err);
831
832        // acos resulting to NAN
833        let n = TWO.acos();
834        assert!(n.is_nan());
835        let n = TWO.inv_sign().acos();
836        assert!(n.is_nan());
837    }
838
839    #[test]
840    fn test_lib_tan_atan() {
841        // tan, atan
842        for _ in 0..10000 {
843            let num = random_normal_float(3, 40);
844            let t = num.tan();
845            let a = t.atan();
846            assert!(HALF_PI.cmp(&a).unwrap() >= 0 && HALF_PI.inv_sign().cmp(&a).unwrap() <= 0);
847            if num.abs().cmp(&HALF_PI).unwrap() <= 0 {
848                assert_small(&num.sub(&a), 2, &num);
849            } else {
850                let mut sub1 = num.add(&a).abs().div(&PI).frac();
851                let mut sub2 = num.sub(&a).abs().div(&PI).frac();
852                if sub1.get_mantissa_len() > 2 {
853                    sub1 = ONE.sub(&sub1);
854                }
855                if sub2.get_mantissa_len() > 2 {
856                    sub2 = ONE.sub(&sub2);
857                }
858                assert_small_any(&sub1, &sub2, 2, &num);
859            }
860        }
861
862        // near pi/2, -pi/2
863        let eps = BigFloat::from_f64(1.0e-39);
864        let n = HALF_PI.tan();
865        assert!(n.get_mantissa_len() as i8 + n.get_exponent() > 39);
866        let n = HALF_PI.sub(&eps).tan();
867        assert!(n.get_mantissa_len() as i8 + n.get_exponent() > 39 && n.is_positive());
868        let n = HALF_PI.inv_sign().tan();
869        assert!(n.get_mantissa_len() as i8 + n.get_exponent() > 39);
870        let n = HALF_PI.sub(&eps).inv_sign().tan();
871        assert!(n.get_mantissa_len() as i8 + n.get_exponent() > 39 && n.is_negative());
872
873        // atan for large negative and large positive.
874        let n = MAX.atan();
875        assert_small(&HALF_PI.sub(&n), 2, &HALF_PI);
876        let n = MIN.atan();
877        assert_small(&HALF_PI.inv_sign().sub(&n), 2, &HALF_PI.inv_sign());
878
879        let mut n = MAX;
880        n.set_exponent(n.get_exponent() - (DECIMAL_POSITIONS as i8));
881        let n = n.atan();
882        assert_small(&HALF_PI.sub(&n), 2, &HALF_PI);
883
884        let mut n = MIN;
885        n.set_exponent(n.get_exponent() - (DECIMAL_POSITIONS as i8));
886        let n = n.atan();
887        assert_small(&HALF_PI.inv_sign().sub(&n), 2, &HALF_PI.inv_sign());
888    }
889
890    #[test]
891    fn test_lib_sinh_asinh() {
892        // sinh, asinh
893        for _ in 0..10000 {
894            let num = random_normal_float(91, 127);
895            let s = num.sinh();
896            let a = s.asinh();
897            if !a.is_inf() && !s.is_inf() {
898                assert_small(&num.sub(&a), 2, &num);
899            }
900        }
901
902        // sinh of MAX
903        let n = MAX.sinh();
904        assert!(n.is_inf_pos());
905        let n = MIN.sinh();
906        assert!(n.is_inf_neg());
907
908        // asinh of MAX
909        let n = MAX.asinh();
910        assert!(n.cmp(&MAX).unwrap() < 0);
911        let n = n.sinh();
912        assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 2, &MAX);
913        let n = MIN.asinh();
914        assert!(n.cmp(&MIN).unwrap() > 0);
915        let n = n.sinh();
916        assert_small_if_not(n.is_inf_neg(), &MIN.sub(&n), 2, &MIN);
917    }
918
919    #[test]
920    fn test_lib_cosh_acosh() {
921        // cosh, acosh
922        for _ in 0..10000 {
923            let num = random_normal_float(4, 40);
924            let s = num.cosh();
925            let a = s.acosh();
926            if !a.is_inf() && !s.is_inf() {
927                assert_small(&num.abs().sub(&a), 3, &num);
928            }
929        }
930
931        // cosh extremums at 0
932        let exp_err = -(DECIMAL_POSITIONS as i8) - 37;
933        let eps = BigFloat::from_f64(1.0e-19);
934        test_extremum(BigFloat::cosh, &ZERO, &ONE, 3, 2, 2, &eps, exp_err);
935
936        // cosh of MAX
937        let n = MAX.cosh();
938        assert!(n.is_inf_pos());
939        let n = MIN.cosh();
940        assert!(n.is_inf_pos());
941
942        // acosh of MAX
943        let n = MAX.acosh();
944        assert!(n.cmp(&MAX).unwrap() < 0);
945        let n = n.cosh();
946        assert_small_if_not(n.is_inf_pos(), &MAX.sub(&n), 2, &MAX);
947
948        // acosh extrema at 1
949        let exp_err = -(DECIMAL_POSITIONS as i8) - 18;
950        let eps = BigFloat::from_f64(1.0e-39);
951        test_extremum(BigFloat::acosh, &ONE, &ZERO, 1, 2, 2, &eps, exp_err);
952
953        // acosh resulting to NAN
954        let n = ZERO.acosh();
955        assert!(n.is_nan());
956    }
957
958    #[test]
959    fn test_lib_tanh_atanh() {
960        // tanh, atanh
961        for _ in 0..10000 {
962            let num = random_normal_float(88, 127);
963            let s = num.tanh();
964            let a = s.atanh();
965            assert_small(&num.sub(&a), 2, &num);
966        }
967
968        // tanh of MAX
969        let n = MAX.tanh();
970        assert!(n.cmp(&ONE).unwrap() == 0);
971        let n = MIN.tanh();
972        assert!(n.cmp(&ONE.inv_sign()).unwrap() == 0);
973
974        // atanh of 1, -1
975        let n = ONE.atanh();
976        assert!(n.is_inf_pos());
977        let n = ONE.inv_sign().atanh();
978        assert!(n.is_inf_neg());
979
980        // atanh of MAX
981        let n = MAX.atanh();
982        assert!(n.sub(&ZERO).get_mantissa_len() < 2);
983        let n = MIN.atanh();
984        assert!(n.sub(&ZERO).get_mantissa_len() < 2);
985
986        // atanh resulting to NAN
987        let n = TWO.atanh();
988        assert!(n.is_nan());
989        let n = TWO.inv_sign().atanh();
990        assert!(n.is_nan());
991    }
992
993    #[test]
994    fn test_lib_util() {
995        // min, max
996        for _ in 0..1000 {
997            let num1 = random_normal_float(256, 127);
998            let num2 = random_normal_float(256, 127);
999            let n1 = num1.max(&num2);
1000            let n2 = num1.min(&num2);
1001            if num1.cmp(&num2).unwrap() > 0 {
1002                assert!(n1.cmp(&num1).unwrap() == 0);
1003                assert!(n2.cmp(&num2).unwrap() == 0);
1004            } else {
1005                assert!(n1.cmp(&num2).unwrap() == 0);
1006                assert!(n2.cmp(&num1).unwrap() == 0);
1007            }
1008        }
1009
1010        assert!(ONE.max(&INF_POS).cmp(&INF_POS).unwrap() == 0);
1011        assert!(INF_POS.max(&ONE).cmp(&INF_POS).unwrap() == 0);
1012        assert!(ONE.max(&INF_NEG).cmp(&ONE).unwrap() == 0);
1013        assert!(INF_NEG.max(&ONE).cmp(&ONE).unwrap() == 0);
1014        assert!(INF_POS.max(&INF_NEG).cmp(&INF_POS).unwrap() == 0);
1015        assert!(INF_NEG.max(&INF_POS).cmp(&INF_POS).unwrap() == 0);
1016        assert!(NAN.max(&ONE).is_nan());
1017        assert!(NAN.max(&NAN).is_nan());
1018        assert!(ONE.max(&NAN).is_nan());
1019        assert!(INF_POS.max(&NAN).is_nan());
1020        assert!(INF_NEG.max(&NAN).is_nan());
1021        assert!(NAN.max(&INF_POS).is_nan());
1022        assert!(NAN.max(&INF_NEG).is_nan());
1023
1024        assert!(ONE.min(&INF_POS).cmp(&ONE).unwrap() == 0);
1025        assert!(INF_POS.min(&ONE).cmp(&ONE).unwrap() == 0);
1026        assert!(ONE.min(&INF_NEG).cmp(&INF_NEG).unwrap() == 0);
1027        assert!(INF_NEG.min(&ONE).cmp(&INF_NEG).unwrap() == 0);
1028        assert!(INF_POS.min(&INF_NEG).cmp(&INF_NEG).unwrap() == 0);
1029        assert!(INF_NEG.min(&INF_POS).cmp(&INF_NEG).unwrap() == 0);
1030        assert!(NAN.min(&ONE).is_nan());
1031        assert!(NAN.min(&NAN).is_nan());
1032        assert!(ONE.min(&NAN).is_nan());
1033        assert!(INF_POS.min(&NAN).is_nan());
1034        assert!(INF_NEG.min(&NAN).is_nan());
1035        assert!(NAN.min(&INF_POS).is_nan());
1036        assert!(NAN.min(&INF_NEG).is_nan());
1037
1038        // clamp
1039        for _ in 0..1000 {
1040            let num1 = random_normal_float(256, 127);
1041            let num2 = random_normal_float(256, 127);
1042            let num3 = random_normal_float(256, 127);
1043            let upper = num1.max(&num2);
1044            let lower = num1.min(&num2);
1045            let n = num3.clamp(&lower, &upper);
1046            assert!(n.cmp(&upper).unwrap() <= 0);
1047            assert!(n.cmp(&lower).unwrap() >= 0);
1048        }
1049
1050        assert!(TWO.clamp(&ONE, &ONE).cmp(&ONE).unwrap() == 0);
1051        assert!(INF_POS.clamp(&ONE, &TWO).cmp(&TWO).unwrap() == 0);
1052        assert!(INF_NEG.clamp(&ONE, &TWO).cmp(&ONE).unwrap() == 0);
1053
1054        assert!(ONE.clamp(&INF_NEG, &INF_POS).cmp(&ONE).unwrap() == 0);
1055        assert!(INF_POS.clamp(&INF_NEG, &INF_POS).cmp(&INF_POS).unwrap() == 0);
1056        assert!(INF_NEG.clamp(&INF_NEG, &INF_POS).cmp(&INF_NEG).unwrap() == 0);
1057
1058        assert!(ONE.clamp(&TWO, &INF_POS).cmp(&TWO).unwrap() == 0);
1059        assert!(INF_POS.clamp(&TWO, &INF_POS).cmp(&INF_POS).unwrap() == 0);
1060        assert!(INF_NEG.clamp(&TWO, &INF_POS).cmp(&TWO).unwrap() == 0);
1061
1062        assert!(TWO.clamp(&INF_NEG, &ONE).cmp(&ONE).unwrap() == 0);
1063        assert!(INF_POS.clamp(&INF_NEG, &ONE).cmp(&ONE).unwrap() == 0);
1064        assert!(INF_NEG.clamp(&INF_NEG, &ONE).cmp(&INF_NEG).unwrap() == 0);
1065
1066        assert!(ZERO.clamp(&INF_POS, &INF_NEG).is_nan());
1067        assert!(ZERO.clamp(&TWO, &ONE).is_nan());
1068        assert!(ZERO.clamp(&INF_POS, &ONE).is_nan());
1069        assert!(ZERO.clamp(&TWO, &INF_NEG).is_nan());
1070
1071        for i in 1..0b111 {
1072            let n1 = if i & 1 == 0 { NAN } else { ONE };
1073            let n2 = if i & 0b10 == 0 { NAN } else { ONE };
1074            let n3 = if i & 0b100 == 0 { NAN } else { ONE };
1075            assert!(n1.clamp(&n2, &n3).is_nan());
1076        }
1077
1078        // signum
1079        assert!(TWO.signum().cmp(&ONE).unwrap() == 0);
1080        assert!(ZERO.signum().cmp(&ONE).unwrap() == 0);
1081        assert!(TWO.inv_sign().signum().cmp(&ONE.inv_sign()).unwrap() == 0);
1082        assert!(NAN.signum().is_nan());
1083        assert!(INF_POS.signum().cmp(&ONE).unwrap() == 0);
1084        assert!(INF_NEG.signum().cmp(&ONE.inv_sign()).unwrap() == 0);
1085    }
1086
1087    #[test]
1088    fn test_lib_stable() {
1089        let mut d1: BigFloat;
1090        let mut d2: BigFloat;
1091        let niter = 10000;
1092
1093        for _ in 0..niter {
1094            d1 = random_normal_float(4, 30);
1095            d2 = random_normal_float(4, 34);
1096            check_stable1(
1097                BigFloat::add,
1098                BigFloat::sub,
1099                &d1,
1100                &d2,
1101                1,
1102                "check stable for add/sub",
1103            );
1104        }
1105
1106        for _ in 0..niter {
1107            d1 = random_normal_float(40, 40);
1108            d2 = random_normal_float(40, 40);
1109            check_stable1(
1110                BigFloat::mul,
1111                BigFloat::div,
1112                &d1,
1113                &d2,
1114                1,
1115                "check stable for mul/div",
1116            );
1117        }
1118
1119        for _ in 0..niter {
1120            d1 = random_normal_float(255, 127).abs();
1121            check_stable2(
1122                BigFloat::sqrt,
1123                |x| x.mul(x),
1124                &d1,
1125                5,
1126                "check stable for sqrt",
1127            );
1128        }
1129
1130        for _ in 0..niter {
1131            d1 = random_normal_float(255, 127).abs();
1132            check_stable2(
1133                BigFloat::cbrt,
1134                |x| x.mul(x).mul(x),
1135                &d1,
1136                5,
1137                "check stable for cbrt",
1138            );
1139        }
1140
1141        for _ in 0..niter {
1142            d1 = random_normal_float(255, 127).abs();
1143            check_stable2(
1144                BigFloat::ln,
1145                BigFloat::exp,
1146                &d1,
1147                5,
1148                "check stable for ln/exp",
1149            );
1150        }
1151
1152        for _ in 0..niter {
1153            d1 = random_normal_float(255, 127).abs();
1154            let d2 = random_normal_float(255, 127).abs();
1155            if !d1.is_zero() && !d1.is_zero() {
1156                check_stable1(
1157                    BigFloat::log,
1158                    |p1, p2| BigFloat::pow(p2, p1),
1159                    &d1,
1160                    &d2,
1161                    1,
1162                    "check stable for log/pow",
1163                );
1164            }
1165        }
1166
1167        for _ in 0..niter {
1168            d1 = random_normal_float(90, 127).abs();
1169            check_stable2(
1170                BigFloat::sin,
1171                BigFloat::asin,
1172                &d1,
1173                10,
1174                "check stable for sin/asin",
1175            );
1176        }
1177
1178        for _ in 0..niter {
1179            d1 = random_normal_float(90, 127).abs();
1180            check_stable2(
1181                BigFloat::cos,
1182                BigFloat::acos,
1183                &d1,
1184                10,
1185                "check stable for cos/acos",
1186            );
1187        }
1188
1189        for _ in 0..niter {
1190            d1 = random_normal_float(3, 40).abs();
1191            check_stable2(
1192                BigFloat::tan,
1193                BigFloat::atan,
1194                &d1,
1195                10,
1196                "check stable for tan/atan",
1197            );
1198        }
1199
1200        for _ in 0..niter {
1201            d1 = random_normal_float(91, 127);
1202            check_stable2(
1203                BigFloat::sinh,
1204                BigFloat::asinh,
1205                &d1,
1206                5,
1207                "check stable for sinh/asinh",
1208            );
1209        }
1210
1211        for _ in 0..niter {
1212            d1 = random_normal_float(91, 127);
1213            check_stable2(
1214                BigFloat::cosh,
1215                BigFloat::acosh,
1216                &d1,
1217                5,
1218                "check stable for cosh/acosh",
1219            );
1220        }
1221
1222        for _ in 0..niter {
1223            d1 = random_normal_float(88, 127);
1224            check_stable2(
1225                BigFloat::tanh,
1226                BigFloat::atanh,
1227                &d1,
1228                5,
1229                "check stable for tanh/atanh",
1230            );
1231        }
1232    }
1233
1234    fn random_f64_exp(exp_range: i32, exp_shift: i32) -> f64 {
1235        let mut f: f64 = random::<f64>() + 1.0;
1236        f *= 10.0f64.powi(random::<i32>().abs() % exp_range - exp_shift);
1237        if random::<i8>() & 1 == 0 {
1238            f = -f;
1239        }
1240        f
1241    }
1242
1243    fn random_normal_float(exp_range: i32, exp_shift: i32) -> BigFloat {
1244        let mut mantissa = [0i16; DECIMAL_PARTS];
1245        for i in 0..DECIMAL_PARTS {
1246            mantissa[i] = (random::<u16>() % DECIMAL_BASE as u16) as i16;
1247        }
1248        if mantissa[DECIMAL_PARTS - 1] == 0 {
1249            mantissa[DECIMAL_PARTS - 1] = (DECIMAL_BASE - 1) as i16;
1250        }
1251        while mantissa[DECIMAL_PARTS - 1] / 1000 == 0 {
1252            mantissa[DECIMAL_PARTS - 1] *= 10;
1253        }
1254        let sign = if random::<i8>() & 1 == 0 { DECIMAL_SIGN_POS } else { DECIMAL_SIGN_NEG };
1255        let exp = (if exp_range != 0 { random::<i32>().abs() % exp_range } else { 0 }) - exp_shift;
1256        BigFloat::from_raw_parts(mantissa, DECIMAL_POSITIONS as i16, sign, exp as i8)
1257    }
1258
1259    // test function extremum
1260    // sides: 3 - both, 2 - left, 1 - right, 0 - none (just exact value).
1261    // x - input value, y - expected value.
1262    // err - number of digits in mantissa allowed for error at point near extremum.
1263    // exact_err - number of digits in mantissa allowed for error at extremum.
1264    // exp_err - compare exponent of difference to exp_err
1265    fn test_extremum(
1266        f: fn(&BigFloat) -> BigFloat,
1267        x: &BigFloat,
1268        y: &BigFloat,
1269        sides: u8,
1270        exact_err: usize,
1271        err: usize,
1272        eps: &BigFloat,
1273        exp_err: i8,
1274    ) {
1275        let n = f(x);
1276        assert_close(&n, y, exact_err, exp_err);
1277
1278        if sides & 1 != 0 {
1279            let x1 = x.add(eps);
1280            let n = f(&x1);
1281            assert_close(&n, y, err, exp_err);
1282        }
1283
1284        if sides & 2 != 0 {
1285            let x1 = x.sub(eps);
1286            let n = f(&x1);
1287            assert_close(&n, y, err, exp_err);
1288        }
1289    }
1290
1291    // assert n is close to y
1292    fn assert_close(n: &BigFloat, y: &BigFloat, err: usize, exp_err: i8) {
1293        assert!(
1294            n.cmp(y).unwrap() == 0
1295                || n.sub(y).get_mantissa_len() < err
1296                || n.sub(y).get_exponent() < exp_err
1297        );
1298    }
1299
1300    // assert v is as large as ndigits digits of scale_num's mantissa, i.e. make sure v << scale_num
1301    fn assert_small(v: &BigFloat, ndigits: i32, scale_num: &BigFloat) {
1302        let eps = get_assert_small_eps(ndigits, scale_num);
1303        assert!(v.abs().cmp(&eps).unwrap() < 0);
1304    }
1305
1306    // assert_small if not other condition true
1307    fn assert_small_if_not(cond: bool, v: &BigFloat, ndigits: i32, scale_num: &BigFloat) {
1308        if !cond {
1309            assert_small(v, ndigits, scale_num);
1310        }
1311    }
1312
1313    // assert_small either for v1 or for v2
1314    fn assert_small_any(v1: &BigFloat, v2: &BigFloat, ndigits: i32, scale_num: &BigFloat) {
1315        let eps = get_assert_small_eps(ndigits, scale_num);
1316        assert!(v1.abs().cmp(&eps).unwrap() < 0 || v2.abs().cmp(&eps).unwrap() < 0);
1317    }
1318
1319    // prepare epsilon value
1320    fn get_assert_small_eps(ndigits: i32, scale_num: &BigFloat) -> BigFloat {
1321        let mut eps = ONE;
1322        let e = ndigits - eps.get_mantissa_len() as i32 - (DECIMAL_POSITIONS as i32)
1323            + scale_num.get_exponent() as i32
1324            + scale_num.get_mantissa_len() as i32;
1325        if e < DECIMAL_MIN_EXPONENT as i32 {
1326            eps = MIN_POSITIVE;
1327            if e + DECIMAL_POSITIONS as i32 > DECIMAL_MIN_EXPONENT as i32 {
1328                eps.set_exponent((e + DECIMAL_POSITIONS as i32) as i8);
1329            }
1330        } else if e > DECIMAL_MAX_EXPONENT as i32 {
1331            eps = MAX;
1332        } else {
1333            eps.set_exponent(e as i8);
1334        }
1335        eps
1336    }
1337
1338    // make sure consecutive application of f and inverse f do not diverge
1339    fn check_stable1(
1340        f: fn(&BigFloat, &BigFloat) -> BigFloat,
1341        finv: fn(&BigFloat, &BigFloat) -> BigFloat,
1342        d1: &BigFloat,
1343        d2: &BigFloat,
1344        niter: usize,
1345        #[allow(unused_variables)] msg: &str,
1346    ) {
1347        let mut n1 = *d1;
1348        let mut n2;
1349        for _ in 0..niter {
1350            let p = f(&n1, d2);
1351            n2 = finv(&p, d2);
1352            n1 = n2;
1353        }
1354        n2 = f(&n1, d2);
1355        n2 = finv(&n2, d2);
1356        let assertion = (n1.is_nan() && n2.is_nan()) || n1.cmp(&n2).unwrap() == 0;
1357        #[cfg(feature = "std")]
1358        if !assertion {
1359            println!("{}", msg);
1360        }
1361        assert!(assertion);
1362    }
1363
1364    // make sure consecutive application of f and inverse f do not diverge
1365    fn check_stable2(
1366        f: fn(&BigFloat) -> BigFloat,
1367        finv: fn(&BigFloat) -> BigFloat,
1368        d1: &BigFloat,
1369        niter: usize,
1370        #[allow(unused_variables)] msg: &str,
1371    ) {
1372        let mut n1 = *d1;
1373        let mut n2;
1374        for _ in 0..niter {
1375            let p = f(&n1);
1376            n2 = finv(&p);
1377            n1 = n2;
1378        }
1379        n2 = f(&n1);
1380        n2 = finv(&n2);
1381        let assertion = (n1.is_nan() && n2.is_nan()) || n1.cmp(&n2).unwrap() == 0;
1382        #[cfg(feature = "std")]
1383        if !assertion {
1384            println!("{}", msg);
1385        }
1386        assert!(assertion);
1387    }
1388}
num_bigfloat/lib.rs

num_bigfloat/
lib.rs
