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

// 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/
//
// NOTES
//
// This routine is a translation into C of a Fortran subroutine
// by W. Fullerton of Los Alamos Scientific Laboratory.
// Some modifications have been made so that the routines
// conform to the IEEE 754 standard.

/// Returns the value of x rounded to "digits" significant
/// decimal digits.
/// R's  signif(x, digits)   via   Math2(args, fprec) in  ../main/arithmetic.c :
pub fn fprec(mut x: f64, digits: usize) -> f64 {
    // Max.expon. of 10 (w/o denormalizing or overflow; = R's  trunc( log10(.Machine$double.xmax) )
    let max10e = f64::MAX_10_EXP; // == 308 ("IEEE")

    if x.is_nan() {
        return x;
    }
    if !x.is_finite() {
        return x;
    }
    if x == 0.0 {
        return x;
    }
    let mut dig = digits as i32;
    if dig > 22 {
        return x;
    } else if dig < 1 {
        dig = 1
    }
    let mut sgn = 1.0;
    if x < 0.0 {
        sgn = -sgn;
        x = -x
    }
    let l10 = x.log10();
    let mut e10 = (dig - 1) - l10.floor() as i32;
    if l10.abs() < max10e as f64 - 2.0 {
        let mut p10 = 1.0;
        if e10 > max10e {
            /* numbers less than 10^(dig-1) * 1e-308 */
            p10 = 10.0_f64.powi(e10 - max10e);
            e10 = max10e
        }
        if e10 > 0 {
            /* Try always to have pow >= 1
            and so exactly representable */
            let pow10 = 10.0_f64.powi(e10);
            sgn * ((x * pow10 * p10).round() / pow10) / p10
        } else {
            let pow10 = 10.0_f64.powi(-e10);
            sgn * ((x / pow10).round() * pow10)
        }
    } else {
        let do_round = max10e as f64 - l10 >= 10.0_f64.powi(-dig);
        let e2 = dig + (if e10 > 0 { 1 } else { -(1) }) * 22;
        let p10 = 10.0_f64.powi(e2);
        x *= p10;
        let p_10 = 10.0_f64.powi(e10 - e2);
        x *= p_10;
        /*-- p10 * p_10 = 10 ^ e10 */
        if do_round {
            x += 0.5
        }
        x = x.floor() / p10;
        sgn * x / p_10
    }
}