| /* |
| * 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 <config.h> |
| #endif |
| |
| #ifdef ENABLE_NLS |
| #include <libintl.h> |
| #define _(String) gettext (String) |
| #else |
| #define _(String) (String) |
| #endif |
| |
| #include <math.h> |
| #include <float.h> /* DBL_MAX */ |
| #include <R_ext/Arith.h> /* NA handling */ |
| #include <Rmath.h> |
| #include <R_ext/Error.h> |
| |
| #include <R_ext/Applic.h> |
| |
| #ifdef DEBUGpr |
| # include <R_ext/Print.h> |
| #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 |
| } |