/* * R : A Computer Language for Statistical Data Analysis * Copyright (C) 1995-2017 The R Core Team * * 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/ */ /* Pretty Intervals * ---------------- * Constructs m "pretty" values which cover the given interval *lo <= *up * m ~= *ndiv + 1 (i.e., ndiv := approximate number of INTERVALS) * * It is not quite clear what should happen for *lo = *up; * S itself behaves quite funilly, then. * * In my opinion, a proper 'pretty' should always ensure * *lo < *up, and hence *ndiv >=1 in the result. * However, in S and here, we allow *lo == *up, and *ndiv = 0. * Note however, that we are NOT COMPATIBLE to S. [Martin M.] * * NEW (0.63.2): ns, nu are double (==> no danger of integer overflow) * * We determine * if the interval (up - lo) is ``small'' [<==> i_small == TRUE, below]. * For the ``i_small'' situation, there is a parameter shrink_sml, * the factor by which the "scale" is shrunk. ~~~~~~~~~~ * It is advisable to set it to some (smaller) integer power of 2, * since this enables exact floating point division. */ #ifdef HAVE_CONFIG_H #include #endif #ifdef ENABLE_NLS #include #define _(String) gettext (String) #else #define _(String) (String) #endif #include #include #include #include #include #ifdef DEBUGpr # include #endif #ifdef HAVE_VISIBILITY_ATTRIBUTE # define attribute_hidden __attribute__ ((visibility ("hidden"))) #else # define attribute_hidden #endif attribute_hidden double R_pretty(double *lo, double *up, int *ndiv, int min_n, double shrink_sml, double high_u_fact[], int eps_correction, int return_bounds) { /* From version 0.65 on, we had rounding_eps := 1e-5, before, r..eps = 0 * then, 1e-7 was consistent with seq.default() and seq.int() till 2010-02-03, * where it was changed to 1e-10 for seq*(), and in 2017-08-14 for pretty(): */ #define rounding_eps 1e-10 #define h high_u_fact[0] #define h5 high_u_fact[1] double dx, cell, unit, base, U; double ns, nu; int k; Rboolean i_small; dx = *up - *lo; /* cell := "scale" here */ if(dx == 0 && *up == 0) { /* up == lo == 0 */ cell = 1; i_small = TRUE; } else { cell = fmax2(fabs(*lo),fabs(*up)); /* U = upper bound on cell/unit */ U = 1 + ((h5 >= 1.5*h+.5) ? 1/(1+h) : 1.5/(1+h5)); U *= imax2(1,*ndiv) * DBL_EPSILON; // avoid overflow for large ndiv /* added times 3, as several calculations here */ i_small = dx < cell * U * 3; } /*OLD: cell = FLT_EPSILON+ dx / *ndiv; FLT_EPSILON = 1.192e-07 */ if(i_small) { if(cell > 10) cell = 9 + cell/10; cell *= shrink_sml; if(min_n > 1) cell /= min_n; } else { cell = dx; if(*ndiv > 1) cell /= *ndiv; } if(cell < 20*DBL_MIN) { warning(_("Internal(pretty()): very small range.. corrected")); cell = 20*DBL_MIN; } else if(cell * 10 > DBL_MAX) { warning(_("Internal(pretty()): very large range.. corrected")); cell = .1*DBL_MAX; } /* NB: the power can be negative and this relies on exact calculation, which glibc's exp10 does not achieve */ base = pow(10.0, floor(log10(cell))); /* base <= cell < 10*base */ /* unit : from { 1,2,5,10 } * base * such that |u - cell| is small, * favoring larger (if h > 1, else smaller) u values; * favor '5' more than '2' if h5 > h (default h5 = .5 + 1.5 h) */ unit = base; if((U = 2*base)-cell < h*(cell-unit)) { unit = U; if((U = 5*base)-cell < h5*(cell-unit)) { unit = U; if((U =10*base)-cell < h*(cell-unit)) unit = U; }} /* Result: c := cell, u := unit, b := base * c in [ 1, (2+ h) /(1+h) ] b ==> u= b * c in ( (2+ h)/(1+h), (5+2h5)/(1+h5)] b ==> u= 2b * c in ( (5+2h)/(1+h), (10+5h) /(1+h) ] b ==> u= 5b * c in ((10+5h)/(1+h), 10 ) b ==> u=10b * * ===> 2/5 *(2+h)/(1+h) <= c/u <= (2+h)/(1+h) */ ns = floor(*lo/unit+rounding_eps); nu = ceil (*up/unit-rounding_eps); #ifdef DEBUGpr REprintf("pretty(lo=%g,up=%g,ndiv=%d,min_n=%d,shrink=%g,high_u=(%g,%g)," "eps=%d)\n\t dx=%g; is.small:%d. ==> cell=%g; unit=%g\n", *lo, *up, *ndiv, min_n, shrink_sml, h, h5, eps_correction, dx, (int)i_small, cell, unit); #endif if(eps_correction && (eps_correction > 1 || !i_small)) { if(*lo != 0.) *lo *= (1- DBL_EPSILON); else *lo = -DBL_MIN; if(*up != 0.) *up *= (1+ DBL_EPSILON); else *up = +DBL_MIN; } #ifdef DEBUGpr if(ns*unit > *lo) REprintf("\t ns= %.0f -- while(ns*unit > *lo) ns--;\n", ns); #endif while(ns*unit > *lo + rounding_eps*unit) ns--; #ifdef DEBUGpr if(nu*unit < *up) REprintf("\t nu= %.0f -- while(nu*unit < *up) nu++;\n", nu); #endif while(nu*unit < *up - rounding_eps*unit) nu++; k = (int)(0.5 + nu - ns); if(k < min_n) { /* ensure that nu - ns == min_n */ #ifdef DEBUGpr REprintf("\tnu-ns=%.0f-%.0f=%d SMALL -> ensure nu-ns= min_n=%d\n", nu,ns, k, min_n); #endif k = min_n - k; if(ns >= 0.) { nu += k/2; ns -= k/2 + k%2;/* ==> nu-ns = old(nu-ns) + min_n -k = min_n */ } else { ns -= k/2; nu += k/2 + k%2; } *ndiv = min_n; } else { *ndiv = k; } if(return_bounds) { /* if()'s to ensure that result covers original range */ if(ns * unit < *lo) *lo = ns * unit; if(nu * unit > *up) *up = nu * unit; } else { *lo = ns; *up = nu; } #ifdef DEBUGpr REprintf("\t ns=%.0f ==> lo=%g\n", ns, *lo); REprintf("\t nu=%.0f ==> up=%g ==> ndiv = %d\n", nu, *up, *ndiv); #endif return unit; #undef h #undef h5 }