1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 2000-2019 The R Core Team
4 * Copyright (C) 1998 Ross Ihaka
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * https://www.R-project.org/Licenses/
19 *
20 * SYNOPSIS
21 *
22 * #include <Rmath.h>
23 * double fprec(double x, double digits);
24 *
25 * DESCRIPTION
26 *
27 * Returns the value of x rounded to "digits" significant
28 * decimal digits.
29 *
30 * NOTES
31 *
32 * This routine is a translation into C of a Fortran subroutine
33 * by W. Fullerton of Los Alamos Scientific Laboratory.
34 * Some modifications have been made so that the routines
35 * conform to the IEEE 754 standard.
36 */
37
38 #include <config.h>
39 #include "nmath.h"
40
41 /* Improvements by Martin Maechler, May 1997;
42 further ones, Feb.2000:
43 Replace pow(x, (double)i) by R_pow_di(x, i) {and use int dig} */
44
45 #define MAX_DIGITS 22
46 /* was till R 0.99: DBL_DIG := digits of precision of a double, usually 15 */
47 /* FIXME: Hmm, have quite a host of these:
48
49 1) ./fround.c uses much more (sensibly!) ``instead''
50 2) ../main/coerce.c & ../main/deparse.c have DBL_DIG directly
51 3) ../main/options.c has #define MAX_DIGITS 22 for options(digits)
52
53 Really should decide on a (config.h dependent?) global MAX_DIGITS.
54 --MM--
55 */
56
57 // R's signif(x, digits) via Math2(args, fprec) in ../main/arithmetic.c :
fprec(double x,double digits)58 double fprec(double x, double digits)
59 {
60 double l10, pow10, sgn, p10, P10;
61 int e10, e2, do_round, dig;
62 // Max.expon. of 10 (w/o denormalizing or overflow; = R's trunc( log10(.Machine$double.xmax) )
63 const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE")
64
65 if (ISNAN(x) || ISNAN(digits))
66 return x + digits;
67 if (!R_FINITE(x)) return x;
68 if (!R_FINITE(digits)) {
69 if(digits > 0.0) return x;
70 else digits = 1.0;
71 }
72 if(x == 0) return x;
73 dig = (int)round(digits);
74 if (dig > MAX_DIGITS) {
75 return x;
76 } else if (dig < 1)
77 dig = 1;
78
79 sgn = 1.0;
80 if(x < 0.0) {
81 sgn = -sgn;
82 x = -x;
83 }
84 l10 = log10(x);
85 e10 = (int)(dig-1-floor(l10));
86 if(fabs(l10) < max10e - 2) {
87 p10 = 1.0;
88 if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */
89 p10 = R_pow_di(10., e10-max10e);
90 e10 = max10e;
91 }
92 if(e10 > 0) { /* Try always to have pow >= 1
93 and so exactly representable */
94 pow10 = R_pow_di(10., e10);
95 return(sgn*(nearbyint((x*pow10)*p10)/pow10)/p10);
96 } else {
97 pow10 = R_pow_di(10., -e10);
98 return(sgn*(nearbyint((x/pow10))*pow10));
99 }
100 } else { /* -- LARGE or small -- */
101 do_round = max10e - l10 >= R_pow_di(10., -dig);
102 e2 = dig + ((e10>0)? 1 : -1) * MAX_DIGITS;
103 p10 = R_pow_di(10., e2); x *= p10;
104 P10 = R_pow_di(10., e10-e2); x *= P10;
105 /*-- p10 * P10 = 10 ^ e10 */
106 if(do_round) x += 0.5;
107 x = floor(x) / p10;
108 return(sgn*x/P10);
109 }
110 }
111