The R Project SVN R

Rev

Rev 79203 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 79203 Rev 88192
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 2000-2019 The R Core Team
 3 
 *  Copyright (C) 2000-2025 The R Core Team
4 
 *  Copyright (C) 1998 Ross Ihaka
 4 
 *  Copyright (C) 1998 Ross Ihaka
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 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
 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
 8 
 *  the Free Software Foundation; either version 2 of the License, or
Line 55... Line 55...
55 
     */
 55 
     */
56 
 
 56 
 
57 
// R's  signif(x, digits)   via   Math2(args, fprec) in  ../main/arithmetic.c :
 57 
// R's  signif(x, digits)   via   Math2(args, fprec) in  ../main/arithmetic.c :
58 
double fprec(double x, double digits)
 58 
double fprec(double x, double digits)
59 
{
 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))
 60 
    if (ISNAN(x) || ISNAN(digits))
66 
	return x + digits;
 61 
	return x + digits;
67 
    if (!R_FINITE(x)) return x;
 62 
    if (!R_FINITE(x)) return x;
68 
    if (!R_FINITE(digits)) {
 63 
    if (!R_FINITE(digits)) {
69 
	if(digits > 0.0) return x;
 64 
	if(digits > 0.) return x;
70 
	else digits = 1.0;
 65 
	else digits = 1.;
71 
    }
 66 
    }
72 
    if(x == 0) return x;
 67 
    if(x == 0) return x;
- 
 
 68 
 
73 
    dig = (int)round(digits);
 69 
    int dig = (int)round(digits);
74 
    if (dig > MAX_DIGITS) {
 70 
    if (dig > MAX_DIGITS) {
75 
	return x;
 71 
	return x;
76 
    } else if (dig < 1)
 72 
    } else if (dig < 1)
77 
	dig = 1;
 73 
	dig = 1;
78 
 
 74 
 
79 
    sgn = 1.0;
 75 
    double sgn = 1.;
80 
    if(x < 0.0) {
 76 
    if(x < 0.0) {
81 
	sgn = -sgn;
 77 
	sgn = -sgn;
82 
	x = -x;
 78 
	x = -x;
83 
    }
 79 
    }
84 
    l10 = log10(x);
 80 
    double l10 = log10(x);
85 
    e10 = (int)(dig-1-floor(l10));
 81 
    int e10 = dig-1 - (int)floor(l10);
- 
 
 82 
    // Max.expon. of 10 (w/o denormalizing or overflow; = R's trunc( log10(.Machine$double.xmax) ):
- 
 
 83 
    const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE")
86 
    if(fabs(l10) < max10e - 2) {
 84 
    if(fabs(l10) < max10e - 2) {
87 
	p10 = 1.0;
 85 
	double p10 = 1.;
88 
	if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */
 86 
	if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */
89 
	    p10 =  R_pow_di(10., e10-max10e);
 87 
	    p10 =  R_pow_di(10., e10-max10e);
90 
	    e10 = max10e;
 88 
	    e10 = max10e;
91 
	}
 89 
	}
- 
 
 90 
	double pow10;
92 
	if(e10 > 0) { /* Try always to have pow >= 1
 91 
	if(e10 > 0) { /* Try always to have pow >= 1
93 
			 and so exactly representable */
 92 
			 and so exactly representable */
94 
	    pow10 = R_pow_di(10., e10);
 93 
	    pow10 = R_pow_di(10., e10);
95 
	    return(sgn*(nearbyint((x*pow10)*p10)/pow10)/p10);
 94 
	    return(sgn*(nearbyint((x*pow10)*p10)/pow10)/p10);
96 
	} else {
 95 
	} else {
97 
	    pow10 = R_pow_di(10., -e10);
 96 
	    pow10 = R_pow_di(10., -e10);
98 
	    return(sgn*(nearbyint((x/pow10))*pow10));
 97 
	    return(sgn*(nearbyint((x/pow10))*pow10));
99 
	}
 98 
	}
100 
    } else { /* -- LARGE or small -- */
 99 
    } else { /* -- LARGE or small -- */
101 
	do_round = max10e - l10	 >= R_pow_di(10., -dig);
 100 
	bool do_round = log10(DBL_MAX) - l10 >= R_pow_di(10., -dig); /* e.g. signif(1.09e308, 2) */
102 
	e2 = dig + ((e10>0)? 1 : -1) * MAX_DIGITS;
 101 
	int e2 = dig + ((e10 > 0)? 1 : -1) * MAX_DIGITS;
- 
 
 102 
	double
103 
	p10 = R_pow_di(10., e2);	x *= p10;
 103 
	    p10 = R_pow_di(10., e2),
104 
	P10 = R_pow_di(10., e10-e2);	x *= P10;
 104 
	    P10 = R_pow_di(10., e10-e2);
- 
 
 105 
	x *= p10;
- 
 
 106 
	x *= P10;
105 
	/*-- p10 * P10 = 10 ^ e10 */
 107 
	/*-- p10 * P10 = 10 ^ e10 */
106 
	if(do_round) x += 0.5;
 108 
	if(do_round) x += 0.5;
107 
	x = floor(x) / p10;
 109 
	x = floor(x) / p10;
108 
	return(sgn*x/P10);
 110 
	return(sgn*x/P10);
109 
    }
 111 
    }