The R Project SVN R

Rev

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

Rev 71083 Rev 74317
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) 1998 Ross Ihaka
 3 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000-11 The R Core Team
 4 
 *  Copyright (C) 2000-2018 The R Core Team
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
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 29... Line 29...
29 
 */
 29 
 */
30 
 
 30 
 
31 
#include <config.h> /* needed for HAVE_* */
 31 
#include <config.h> /* needed for HAVE_* */
32 
#include "nmath.h"
 32 
#include "nmath.h"
33 
 
 33 
 
34 
 
 - 
 
35 
/*  nearbyint is C99, so all platforms should have it (and AFAIK, all do) */
 - 
 
36 
#ifdef HAVE_NEARBYINT
 - 
 
37 
# define R_rint nearbyint
 - 
 
38 
#elif defined(HAVE_RINT)
 - 
 
39 
# define R_rint rint
 - 
 
40 
#else
 - 
 
41 
# define R_rint private_rint
 - 
 
42 
# include "nmath2.h" // for private_rint
 - 
 
43 
 
 - 
 
44 
/* also used potentially in fprec.c and main/format.c */
 - 
 
45 
double attribute_hidden private_rint(double x)
 - 
 
46 
{
 - 
 
47 
    double tmp, sgn = 1.0;
 - 
 
48 
    long ltmp;
 - 
 
49 
 
 - 
 
50 
    if (x != x) return x;			/* NaN */
 - 
 
51 
 
 - 
 
52 
    if (x < 0.0) {
 - 
 
53 
	x = -x;
 - 
 
54 
	sgn = -1.0;
 - 
 
55 
    }
 - 
 
56 
 
 - 
 
57 
    if(x < (double) LONG_MAX) { /* in <limits.h> is architecture dependent */
 - 
 
58 
	ltmp = x + 0.5;
 - 
 
59 
	/* implement round to even */
 - 
 
60 
	if(fabs(x + 0.5 - ltmp) < 10*DBL_EPSILON
 - 
 
61 
	   && (ltmp % 2 == 1)) ltmp--;
 - 
 
62 
	tmp = ltmp;
 - 
 
63 
    } else {
 - 
 
64 
	/* ignore round to even: too small a point to bother */
 - 
 
65 
	tmp = floor(x + 0.5);
 - 
 
66 
    }
 - 
 
67 
    return sgn * tmp;
 - 
 
68 
}
 - 
 
69 
#endif
 - 
 
70 
 
 - 
 
71 
double fround(double x, double digits) {
 34 
double fround(double x, double digits) {
72 
#define MAX_DIGITS DBL_MAX_10_EXP
 35 
#define MAX_DIGITS DBL_MAX_10_EXP
73 
    /* = 308 (IEEE); was till R 0.99: (DBL_DIG - 1) */
 36 
    /* = 308 (IEEE); was till R 0.99: (DBL_DIG - 1) */
74 
    /* Note that large digits make sense for very small numbers */
 37 
    /* Note that large digits make sense for very small numbers */
75 
    LDOUBLE pow10, sgn, intx;
 38 
    LDOUBLE pow10, sgn, intx;
Line 88... Line 51...
88 
	sgn = -1.;
 51 
	sgn = -1.;
89 
	x = -x;
 52 
	x = -x;
90 
    } else
 53 
    } else
91 
	sgn = 1.;
 54 
	sgn = 1.;
92 
    if (dig == 0) {
 55 
    if (dig == 0) {
93 
	return (double)(sgn * R_rint(x));
 56 
	return (double)(sgn * nearbyint(x));
94 
    } else if (dig > 0) {
 57 
    } else if (dig > 0) {
95 
        pow10 = R_pow_di(10., dig);
 58 
        pow10 = R_pow_di(10., dig);
96 
	intx = floor(x);
 59 
	intx = floor(x);
97 
	return (double)(sgn * (intx + R_rint((double)((x-intx) * pow10)) / pow10));
 60 
	return (double)(sgn * (intx + nearbyint((double)((x-intx) * pow10)) / pow10));
98 
    } else {
 61 
    } else {
99 
        pow10 = R_pow_di(10., -dig);
 62 
        pow10 = R_pow_di(10., -dig);
100 
        return (double)(sgn * R_rint((double)(x/pow10)) * pow10);
 63 
        return (double)(sgn * nearbyint((double)(x/pow10)) * pow10);
101 
    }
 64 
    }
102 
}
 65 
}