/* * R : A Computer Language for Statistical Data Analysis * Copyright (C) 1995-2021 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 /* DBL_MAX */ #include /* NA handling */ #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, const double high_u_fact[], // = (h, h5, f_min) below 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 // (h, h5, f_min) = c(high.u.bias, u5.bias, f.min) in base::pretty.default(): #define h high_u_fact[0] #define h5 high_u_fact[1] #define f_min high_u_fact[2] double // save input boundaries lo_ = *lo, up_ = *up, dx = up_ - lo_, cell, U; Rboolean i_small; /* 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; } #ifdef DEBUGpr REprintf("R_pretty(lo=%g,up=%g,ndiv=%d,min_n=%d,shrink=%g,high_u=(%g,%g,%g),eps=%d,bnds=%d)" "\n\t => dx=%g; i_small:%s. ==> first cell=%g\n", lo_, up_, *ndiv, min_n, shrink_sml, h, h5, min_f, eps_correction, return_bounds, dx, i_small ? "TRUE" : "F", cell); #endif /*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(R_FINITE(dx)) { if(*ndiv > 1) cell /= *ndiv; } else { // up - lo = +Inf (overflow; both are finite) if(*ndiv < 2) { warning(_("R_pretty(): infinite range; *ndiv=%d, should have ndiv >= 2"), *ndiv); } else { cell = up_/(*ndiv) - lo_/(*ndiv); } } } // f_min: arg, default = 2^-20, was 20. till R 4.1.0 (2021-05) #define MAX_F 1.25 // was 10. " " " double subsmall = f_min*DBL_MIN; if(subsmall == 0.) // subnormals underflowing to zero (not yet seen!) subsmall = DBL_MIN; if(cell < subsmall) { // possibly subnormal warning(_("R_pretty(): very small range 'cell'=%g, corrected to %g"), cell, subsmall); cell = subsmall; } else if(cell > DBL_MAX/MAX_F) { warning(_("R_pretty(): very large range 'cell'=%g, corrected to %g"), cell, DBL_MAX/MAX_F); cell = DBL_MAX/MAX_F; } #undef MAX_F /* NB: the power can be negative and this relies on exact calculation, which glibc's exp10 does not achieve */ double 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) */ double 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, b := base, u := unit): * 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+2h5)/(1+h5), (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) */ double ns = floor(lo_/unit+rounding_eps), nu = ceil (up_/unit-rounding_eps); #ifdef DEBUGpr REprintf("\t => final cell=%g; base=%g unit=%g; (ns,nu) = (%g,%g)\n", cell, base, unit, ns, nu); #endif if(eps_correction && (eps_correction > 1 || !i_small)) { // FIXME?: assumes 0 <= lo <= up (what if lo <= up < 0 ?) 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 + rounding_eps*unit) REprintf("\t ns= %.0f -- while(ns*unit > lo + r_eps * u) ns--;\n", ns); #endif while(ns*unit > *lo + rounding_eps*unit) ns--; #ifdef DEBUGpr if(!R_FINITE(ns*unit)) REprintf("\t infinite (ns=%.0f)*(unit=%g) ==> ns++\n", ns, unit); #endif while(!R_FINITE(ns*unit)) ns++; #ifdef DEBUGpr if(nu*unit < *up - rounding_eps*unit) REprintf("\t nu= %.0f -- while(nu*unit < up - r_eps * u) nu++;\n", nu); #endif while(nu*unit < *up - rounding_eps*unit) nu++; #ifdef DEBUGpr if(!R_FINITE(nu*unit)) REprintf("\t infinite (nu=%.0f)*(unit=%g) ==> nu--\n", nu, unit); #endif while(!R_FINITE(nu*unit)) nu--; int k = (int)(0.5 + nu - ns); #ifdef DEBUGpr REprintf(" possibly adjusted (ns,nu) = (%g,%g) ==> k=%d\n", ns, nu, k); #endif 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(lo_ == 0. && ns == 0. && up_ != 0.) { nu += k; } else if(up_ == 0. && nu == 0. && lo_ != 0.) { ns -= k; } else 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) {// used in pretty.default(), ensure result covers original range #ifdef DEBUGpr if(ns * unit < *lo) {REprintf("lo=%g too large, set to ns*unit\n", *lo); *lo = ns * unit;} if(nu * unit > *up) {REprintf("up=%g too small, set to nu*unit\n", *up); *up = nu * unit;} #else if(ns * unit < *lo) *lo = ns * unit; if(nu * unit > *up) *up = nu * unit; #endif } else { // used in graphics GEpretty(), hence grid::grid.pretty() *lo = ns; *up = nu; } #ifdef DEBUGpr REprintf("\t ns=%5.0g ==> lo=%g\n", ns, *lo); REprintf("\t nu=%5.0g ==> up=%g ==> ndiv = %d\n", nu, *up, *ndiv); #endif return unit; #undef h #undef h5 }