/*

 *  R : A Computer Language for Statistical Data Analysis

 *  Copyright (C) 1997--2021  The R Core Team

 *  Copyright (C) 2002--2009  The R Foundation

 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka

 *

 *  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 3 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/

 */

#ifdef HAVE_CONFIG_H

# include <config.h>

#endif

#include <Defn.h>   // Rexp10  (et al)

#include <float.h>  /* for DBL_MAX */

#include <Graphics.h>

#include <Print.h>

#include <Rmath.h> // for imax2

/* used in graphics and grid */

SEXP CreateAtVector(double axp[], const double usr[], int nint, Rboolean logflag)

{

/*  Create an  'at = ...' vector for  axis(.)

 *  i.e., the vector of tick mark locations,

 *  when none has been specified (= default).

 *

 *  axp[0:2] = (x1, x2, nInt), where x1..x2 are the extreme tick marks

 *       {unless in log case, where nInt \in {1,2,3 ; -1,-2,....}

 *        and the `nint' argument is used *instead*.}

 *

 * only if(logflag && axp[2] >= 0)

 *      usr[0:1] is used, additionally

 *

 *  The resulting REAL vector must have length >= 1, ideally >= 2

 */

    SEXP at = R_NilValue;/* -Wall*/

    double dn, rng, small;

    int i, n;

    // "arbitrary" threshold: |delta_tick| / SMALL  is "barely visible" in plot

#define SMALL_F 100.

    if (!logflag || axp[2] < 0) { /* --- linear axis --- Only use axp[] arg. */

  n = (int)(fabs(axp[2]) + 0.25);/* >= 0 */

  dn = imax2(1, n);

  rng = axp[1] - axp[0];

  at = allocVector(REALSXP, n + 1);

  double a_i;

  if(!R_FINITE(rng)) { // need to carefully work around overflow

      double at_ = axp[0]/dn; // 2021-07: "/dn" avoids overflow

      rng = axp[1]/dn - at_;

      small = fabs(rng)/SMALL_F;

#ifdef DEBUG_axis

      REprintf("CreateAtVector(axp=(%g,%g, %g), log=F, diff(*)=Inf: at_=%g, rng=%g, small=%g\n",

         axp[0],axp[1], axp[2], at_, rng, small);

#endif

      int n2 = n/2; // integer division

      for (i = 0; i <= n2; i++) { // from the left

    a_i = axp[0] + i * rng;

    // REprintf(" at[i=%2d]=%g\n", i+1, a_i);

    REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i;

      }

      for (int i2 = 0; i2 < n-n2; i2++) { // from the right

    i = n-i2; // for(i in n:k) where k = n-(n-n2-1) = n2+1

    a_i = axp[1] - i2 * rng;

    // REprintf(" at[i=%2d]=%g\n", i+1, a_i);

    REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i;

      }

  }

  else { // rng is finite (normal case):

      small = fabs(rng)/SMALL_F/dn;

      for (i = 0; i <= n; i++) {

    a_i = axp[0] + (i / dn) * rng;

    REAL(at)[i] = (fabs(a_i) < small) ? 0. : a_i;

      }

  }

    }

    else { /* ------ log axis ----- */

  bool reversed = false;

  double

      umin = usr[0],

      umax = usr[1];

  n = (int)(axp[2] + 0.5);

  /* {xy}axp[2] for 'log': GLpretty() [./graphics.c] sets

     n < 0: very small scale ==> linear axis, above, or

     n = 1,2,3.  see switch() below */

#ifdef DEBUG_axis

  REprintf("CreateAtVector(axp=(%g,%g,%g), usr=(%g,%g), _log_):",

     axp[0],axp[1],axp[2],  usr[0],usr[1]);

#endif

  if (umin > umax) {

      reversed = (axp[0] > axp[1]);

      if (reversed) {

    /* have *reversed* log axis -- whereas

     * the switch(n) { .. } below assumes *increasing* values

     * --> reverse axis direction here, and reverse back at end */

    umin = usr[1];

    umax = usr[0];

    dn = axp[0]; axp[0] = axp[1]; axp[1] = dn;

      }

      else {

    /* can the following still happen... ? */

    warning("CreateAtVector \"log\"(from axis()): "

      "usr[0] = %g > %g = usr[1] !", umin, umax);

      }

  }

  /* allow a fuzz (iff we don't under-/over-flow) since we will do things like 0.2*dn >= umin */

  dn = 1 - 1e-12; if(fabs(umin*dn) >     0.  ) umin *= dn;

  dn = 1 + 1e-12; if(fabs(umax*dn) <= DBL_MAX) umax *= dn;

  dn = axp[0];

  if (dn < DBL_MIN) {/* was 1e-300; now seems too cautious */

      if (dn <= 0) /* real trouble (once for Solaris) later on */

    error("CreateAtVector [log-axis()]: axp[0] = %g < 0!", dn);

      else

    warning("CreateAtVector [log-axis()]: small axp[0] = %g", dn);

  }

  /* You get the 3 cases below by

   *  for (y in 1e-5*c(1,2,8))  plot(y, log = "y")

   */

  switch(n) {

  case 1: /* large range: 1  * 10^k */

  {

      i = (int)(floor(log10(axp[1])) - ceil(log10(axp[0])) + 0.25);

      // want nint intervals, i.e. typically nint+1 breaks :

      int ne = i / nint;

      /* for nint breaks, i.e. typically nint-1 intervals, would be

       * ne = i / imax2(1, nint - 1);  *PLUS* replace s/nint/nint-1/ below !! */

#ifdef DEBUG_axis

      REprintf(" .. case 1: umin,umax= %g,%g;\n  (nint=%d, ne=%d); ",

         umin, umax, nint, ne);

      if (ne < 1) {

    REprintf("ne = %d <= 0 !!\n\t axp[0:1]=(%g,%g) ==> i = %d, nint = %d; ",

       ne, axp[0],axp[1], i, nint);

      }

#endif

      double l10_max = log10(umax),

    d0 = l10_max - log10(dn);

#ifdef DEBUG_axis

      REprintf("exponent diff d0=%g\n", d0);

#endif

      if(ne < 1) ne = 1;

      else // if ne is too large, i.e, the "final tick" is beyond umax, reduce it :

    while(ne > 1 && nint*ne > d0) {

        ne--;

#ifdef DEBUG_axis

        REprintf(" last > umax ==> ne--: ne=%d\n", ne);

#endif

    }

      int k = 1 + ne / 308; // >= 1, typically == 1.

      if(k > 1) {// i.e. ne > 308: 10^ne overflows; must split the multiplication

    ne = k*(ne/k); // <= ne_{previous}

#ifdef DEBUG_axis

    REprintf(" original ne > 308: split in k=%d parts; new ne=%d\n", k,ne);

#endif

      }

      /* Now, still in exponent-10 range: nint*ne <= d0 = l10_max - log10(dn)

       * If difference (=: d1) is "large", say > 3, increase the first at[] =: d0

       */

      double d1 = d0 - nint*ne; // >= 0

#ifdef DEBUG_axis

      REprintf("expo.diff d0 - nint*ne =: d1=%g\n", d1);

#endif

      d0 = dn;

#define Large_D1 5

//              === was '3' all up into R 4.1.0

      if(d1 > Large_D1) {

    d0 = dn * Rexp10(floor(d1/2));

#ifdef DEBUG_axis

    REprintf("large d1 => d0 := dn * 10 ^ fl(d1/2) = dn * 10^%d = %g\n",

       (int)floor(d1/2), d0);

#endif

      }

      rng = Rexp10((double)ne/k); // = 10^(ne/k) >= 10

      n=0;

      dn=d0;

      while(dn < umax) {

    for(int j=0; j < k; j++)

        dn *= rng;

    n++;

      }

#ifdef DEBUG_axis

      REprintf(" rng:=10^(ne/(k=%d)) = %g => n=%d, final dn=%g\n", k, rng, n, dn);

#endif

      if (!n)

    error("log - axis(), 'at' creation, _LARGE_ range: "

          "invalid {xy}axp or par; nint=%d\n"

          "  axp[0:1]=(%g,%g), usr[0:1]=(%g,%g); i=%d, ni=%d",

          nint, axp[0],axp[1], umin,umax, i,ne);

      at = allocVector(REALSXP, n);

      dn=d0;

      for(int i=0; i < n; i++) {

    REAL(at)[i] = dn;

    for(int j=0; j < k; j++)

        dn *= rng;

      }

      break;

  }

  case 2: /* medium range:  1, 5    * 10^k */

      n = 0;

      if (0.5 * dn >= umin) n++;

#ifdef DEBUG_axis

      REprintf(" .. case 2: (dn, umin,umax, n) = (%g, %g,%g, %d)\n",

         dn, umin, umax, n);

#endif

      for (;;) {

    if (dn > umax) break;

    n++;

    if (5 * dn > umax) break;

    n++;

    dn *= 10;

      }

      if (!n)

    error("log - axis(), 'at' creation, _MEDIUM_ range: "

          "invalid {xy}axp or par;\n"

          "  axp[0]= %g, usr[0:1]=(%g,%g)",

          axp[0], umin,umax);

      at = allocVector(REALSXP, n);

      dn = axp[0];

      n = 0;

      if (0.5 * dn >= umin) REAL(at)[n++] = 0.5 * dn;

      for (;;) {

    if (dn > umax) break;

    REAL(at)[n++] = dn;

    if (5 * dn > umax) break;

    REAL(at)[n++] = 5 * dn;

    dn *= 10;

      }

      break;

  case 3: /* small range:  1,2,5,10 * 10^k */

      n = 0;

      if (0.2 * dn >= umin) n++;

      if (0.5 * dn >= umin) n++;

      for (;;) {

    if (dn > umax) break;

    n++;

    if (2 * dn > umax) break;

    n++;

    if (5 * dn > umax) break;

    n++;

    dn *= 10;

      }

#ifdef DEBUG_axis

      REprintf(" .. case 3: (umin,umax)-usr[*] = (%g, %g); n=%d, dn=%g\n",

         umin-usr[reversed? 1: 0],

         umax-usr[reversed? 0: 1], n, dn);

#endif

      if (!n)

    error("log - axis(), 'at' creation, _SMALL_ range: "

          "invalid {xy}axp or par;\n"

          "  axp[0]= %g, usr[0:1]=(%g,%g)",

          axp[0], umin,umax);

      at = allocVector(REALSXP, n);

      dn = axp[0];

      n = 0;

      if (0.2 * dn >= umin) REAL(at)[n++] = 0.2 * dn;

      if (0.5 * dn >= umin) REAL(at)[n++] = 0.5 * dn;

      for (;;) {

    if (dn > umax) break;

    REAL(at)[n++] = dn;

    if (2 * dn > umax) break;

    REAL(at)[n++] = 2 * dn;

    if (5 * dn > umax) break;

    REAL(at)[n++] = 5 * dn;

    dn *= 10;

      }

      break;

  default:

      error("log - axis(), 'at' creation: INVALID {xy}axp[3] = %g",

      axp[2]);

  }

  if (reversed) {/* reverse back again - last assignment was at[n++]= . */

      for (i = 0; i < n/2; i++) { /* swap( at[i], at[n-i-1] ) : */

    dn = REAL(at)[i];

         REAL(at)[i] = REAL(at)[n-i-1];

                       REAL(at)[n-i-1] = dn;

      }

  }

    } /* linear / log */

    return at;

}