The R Project SVN R

Rev

Rev 42302 | Rev 54763 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 1 
/*
- 
 
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 1998 Ross Ihaka
- 
 
 4 
 *  Copyright (C) 2000, 2001, 2005-2006 The R Development Core Team
- 
 
 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 
 *  http://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 
#ifndef HAVE_RINT
- 
 
 42 
#define USE_BUILTIN_RINT
- 
 
 43 
#endif
- 
 
 44 
 
- 
 
 45 
#ifdef USE_BUILTIN_RINT
- 
 
 46 
#define R_rint private_rint
- 
 
 47 
extern double private_rint(double x);
- 
 
 48 
#else
- 
 
 49 
#define R_rint rint
- 
 
 50 
#endif
- 
 
 51 
/* Improvements by Martin Maechler, May 1997;
- 
 
 52 
   further ones, Feb.2000:
- 
 
 53 
   Replace  pow(x, (double)i) by  R_pow_di(x, i) {and use  int dig} */
- 
 
 54 
 
- 
 
 55 
#define MAX_DIGITS 22
- 
 
 56 
/* was till R 0.99: DBL_DIG := digits of precision of a double, usually 15 */
- 
 
 57 
/* FIXME: Hmm, have quite a host of these:
- 
 
 58 
 
- 
 
 59 
       1) ./fround.c   uses much more (sensibly!) ``instead''
- 
 
 60 
       2) ../main/coerce.c   & ../main/deparse.c have  DBL_DIG	directly
- 
 
 61 
       3) ../main/options.c has	  #define MAX_DIGITS 22	 for options(digits)
- 
 
 62 
 
- 
 
 63 
       Really should decide on a (config.h dependent?) global MAX_DIGITS.
- 
 
 64 
       --MM--
- 
 
 65 
     */
- 
 
 66 
 
- 
 
 67 
 
- 
 
 68 
double fprec(double x, double digits)
- 
 
 69 
{
- 
 
 70 
    double l10, pow10, sgn, p10, P10;
- 
 
 71 
    int e10, e2, do_round, dig;
- 
 
 72 
    /* Max.expon. of 10 (=308.2547) */
- 
 
 73 
    const static int max10e = DBL_MAX_EXP * M_LOG10_2;
- 
 
 74 
 
- 
 
 75 
#ifdef IEEE_754
- 
 
 76 
    if (ISNAN(x) || ISNAN(digits))
- 
 
 77 
	return x + digits;
- 
 
 78 
    if (!R_FINITE(x)) return x;
- 
 
 79 
    if (!R_FINITE(digits)) {
- 
 
 80 
	if(digits > 0) return x;
- 
 
 81 
	else return 0;
- 
 
 82 
    }
- 
 
 83 
#endif
- 
 
 84 
    if(x == 0) return x;
- 
 
 85 
    dig = (int)floor(digits+0.5);
- 
 
 86 
    if (dig > MAX_DIGITS) {
- 
 
 87 
	return x;
- 
 
 88 
    } else if (dig < 1)
- 
 
 89 
	dig = 1;
- 
 
 90 
 
- 
 
 91 
    sgn = 1.0;
- 
 
 92 
    if(x < 0.0) {
- 
 
 93 
	sgn = -sgn;
- 
 
 94 
	x = -x;
- 
 
 95 
    }
- 
 
 96 
    l10 = log10(x);
- 
 
 97 
    e10 = (int)(dig-1-floor(l10));
- 
 
 98 
    if(fabs(l10) < max10e - 2) {
- 
 
 99 
	p10 = 1.0;
- 
 
 100 
	if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */
- 
 
 101 
	    p10 =  R_pow_di(10., e10-max10e);
- 
 
 102 
	    e10 = max10e;
- 
 
 103 
	}
- 
 
 104 
	if(e10 > 0) { /* Try always to have pow >= 1
- 
 
 105 
			 and so exactly representable */
- 
 
 106 
	    pow10 = R_pow_di(10., e10);
- 
 
 107 
	    return(sgn*(R_rint((x*pow10)*p10)/pow10)/p10);
- 
 
 108 
	} else {
- 
 
 109 
	    pow10 = R_pow_di(10., -e10);
- 
 
 110 
	    return(sgn*(R_rint((x/pow10))*pow10));
- 
 
 111 
	}
- 
 
 112 
    } else { /* -- LARGE or small -- */
- 
 
 113 
	do_round = max10e - l10	 >= R_pow_di(10., -dig);
- 
 
 114 
	e2 = dig + ((e10>0)? 1 : -1) * MAX_DIGITS;
- 
 
 115 
	p10 = R_pow_di(10., e2);	x *= p10;
- 
 
 116 
	P10 = R_pow_di(10., e10-e2);	x *= P10;
- 
 
 117 
	/*-- p10 * P10 = 10 ^ e10 */
- 
 
 118 
	if(do_round) x += 0.5;
- 
 
 119 
	x = floor(x) / p10;
- 
 
 120 
	return(sgn*x/P10);
- 
 
 121 
    }
- 
 
 122 
}