The R Project SVN R

/*
 *  Mathlib : A C Library of Special Functions
 *  Copyright (C) 2000-2025 The R Core Team
 *  Copyright (C) 1998 Ross Ihaka
 *
 *  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/
 *
 *  SYNOPSIS
 *
 *    #include <Rmath.h>
 *    double fprec(double x, double digits);
 *
 *  DESCRIPTION
 *
 *    Returns the value of x rounded to "digits" significant
 *    decimal digits.
 *
 *  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.
 */

#include <config.h>
#include "nmath.h"

/* Improvements by Martin Maechler, May 1997;
   further ones, Feb.2000:
   Replace  pow(x, (double)i) by  R_pow_di(x, i) {and use  int dig} */

#define MAX_DIGITS 22
/* was till R 0.99: DBL_DIG := digits of precision of a double, usually 15 */
/* FIXME: Hmm, have quite a host of these:

       1) ./fround.c   uses much more (sensibly!) ``instead''
       2) ../main/coerce.c   & ../main/deparse.c have  DBL_DIG  directly
       3) ../main/options.c has   #define MAX_DIGITS 22  for options(digits)

       Really should decide on a (config.h dependent?) global MAX_DIGITS.
       --MM--
     */

// R's  signif(x, digits)   via   Math2(args, fprec) in  ../main/arithmetic.c :
double fprec(double x, double digits)
{
    if (ISNAN(x) || ISNAN(digits))
    return x + digits;
    if (!R_FINITE(x)) return x;
    if (!R_FINITE(digits)) {
    if(digits > 0.) return x;
    else digits = 1.;
    }
    if(x == 0) return x;

    int dig = (int)round(digits);
    if (dig > MAX_DIGITS) {
    return x;
    } else if (dig < 1)
    dig = 1;

    double sgn = 1.;
    if(x < 0.0) {
    sgn = -sgn;
    x = -x;
    }
    double l10 = log10(x);
    int e10 = dig-1 - (int)floor(l10);
    // Max.expon. of 10 (w/o denormalizing or overflow; = R's trunc( log10(.Machine$double.xmax) ):
    const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE")
    if(fabs(l10) < max10e - 2) {
    double p10 = 1.;
    if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */
        p10 =  R_pow_di(10., e10-max10e);
        e10 = max10e;
    }
    double pow10;
    if(e10 > 0) { /* Try always to have pow >= 1
             and so exactly representable */
        pow10 = R_pow_di(10., e10);
        return(sgn*(nearbyint((x*pow10)*p10)/pow10)/p10);
    } else {
        pow10 = R_pow_di(10., -e10);
        return(sgn*(nearbyint((x/pow10))*pow10));
    }
    } else { /* -- LARGE or small -- */
    bool do_round = log10(DBL_MAX) - l10 >= R_pow_di(10., -dig); /* e.g. signif(1.09e308, 2) */
    int e2 = dig + ((e10 > 0)? 1 : -1) * MAX_DIGITS;
    double
        p10 = R_pow_di(10., e2),
        P10 = R_pow_di(10., e10-e2);
    x *= p10;
    x *= P10;
    /*-- p10 * P10 = 10 ^ e10 */
    if(do_round) x += 0.5;
    x = floor(x) / p10;
    return(sgn*x/P10);
    }
}