r2rs-nmath 0.1.1

Statistics programming for Rust based on R's nmath package
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

// Mathlib : A C Library of Special Functions
// Copyright (C) 1998 Ross Ihaka
// Copyright (C) 2000-2018 The R Core Team
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, a copy is available at
// https://www.R-project.org/Licenses/

use std::ops::Rem;

use num_traits::Float;
use strafe_consts::LOG10_2;

use crate::traits::DPQ;

/// Rounds "x" to "digits" decimal digits.
pub fn fround(mut x: f64, digits: isize) -> f64 {
    let max_digits = f64::MAX_10_EXP as usize + f64::DIGITS as usize; /* typically = 308+15 = 323
                                                                       * was DBL_MAX_10_EXP (= 308, IEEE) till R 3.6.x; before,
                                                                       * was (DBL_DIG - 1)  till R 0.99  */
    let max10e = f64::MAX_10_EXP; // == 308 ("IEEE")

    /* Note that large digits make sense for very small numbers */
    if x.is_nan() {
        return f64::nan();
    }
    if !x.is_finite() {
        return x;
    }

    if digits as usize > max_digits || x == 0.0 {
        return x;
    } else if digits < -max10e as isize {
        // includes -Inf {aka ML_NEGINF}
        return 0.;
    } else if digits == 0 {
        // common
        return x.round();
    }

    let mut sgn = 1.0;
    if x < 0.0 {
        sgn = -1.0;
        x = -x;
    } // now  x > 0
    let l10x = LOG10_2 * (0.5 + x.log2()); // ~= log10(x), but cheaper (presumably)
    if l10x as u32 + digits as u32 > f64::DIGITS {
        // rounding to so many digits that no rounding is needed
        return sgn * x;
    } else {
        let pow10;
        let x10;
        let i10;
        let xd;
        let xu; // x, rounded _d_own or _u_p
        if digits <= max10e as isize {
            // both pow10 := 10^d and x10 := x * pow10 do *not* overflow
            pow10 = 10.0.pow_di(digits);
            x10 = x * pow10;
            i10 = x10.floor();
            xd = i10 / pow10;
            xu = x10.ceil() / pow10;
        } else {
            // DBL_MAX_10_EXP =: max10e < dig <= DBL_DIG - l10x: case of |x| << 1; ~ 10^-305
            let e10 = digits - max10e as isize; // > 0
            let p10 = 10.0.pow_di(e10);
            pow10 = 10.0.pow_di(max10e as isize);
            x10 = (x * pow10) * p10;
            i10 = x10.floor();
            xd = i10 / pow10 / p10;
            xu = x10.ceil() / pow10 / p10;
        }
        let du = xu - x;
        let dd = x - xd;
        //  D =  du - dd
        //  return sgn * ((D < 0 || (is_odd_i10 && D == 0)) ? xu : xd);
        return sgn
            * if du < dd || (i10.rem(2.0) == 1.0 && du == dd) {
                xu
            } else {
                xd
            };
    }
}