/*
  ---------------------------------

   rboxlib.c
     Generate input points

   notes:
     For documentation, see prompt[] of rbox.c
     50 points generated for 'rbox D4'

   WARNING:
     incorrect range if qh_RANDOMmax is defined wrong (user.h)
*/

#include "libqhull.h"  /* First for user.h */
#include "random.h"

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#ifdef _MSC_VER  /* Microsoft Visual C++ */
#pragma warning( disable : 4706)  /* assignment within conditional expression. */
#pragma warning( disable : 4996)  /* this function (strncat,sprintf,strcpy) or variable may be unsafe. */
#endif

#define MAXdim 200
#define PI 3.1415926535897932384

/* ------------------------------ prototypes ----------------*/
int qh_roundi(double a);
void qh_out1(double a);
void qh_out2n(double a, double b);
void qh_out3n(double a, double b, double c);
void qh_outcoord(int iscdd, double *coord, int dim);
void qh_outcoincident(int coincidentpoints, double radius, int iscdd, double *coord, int dim);
void qh_rboxpoints2(char* rbox_command, double **simplex);

void    qh_fprintf_rbox(FILE *fp, int msgcode, const char *fmt, ... );
void    qh_free(void *mem);
void   *qh_malloc(size_t size);
int     qh_rand(void);
void    qh_srand(int seed);


/* ------------------------------ globals -------------------*/

/* No state is carried between rbox requests */
typedef struct rboxT rboxT;
struct rboxT {
  FILE *fout;
  FILE *ferr;
  int isinteger;
  double out_offset;
  jmp_buf errexit;        /* exit label for rboxpoints, defined by setjmp(), called by qh_errexit_rbox() */
  char  jmpXtra[40];      /* extra bytes in case jmp_buf is defined wrong by compiler */
};


int rbox_inuse= 0;
rboxT rbox;

/*---------------------------------

  qh_rboxpoints( fout, ferr, rbox_command )
    Generate points to fout according to rbox options
    Report errors on ferr

  returns:
    0 (qh_ERRnone) on success
    1 (qh_ERRinput) on input error
    4 (qh_ERRmem) on memory error
    5 (qh_ERRqhull) on internal error

  notes:
    To avoid using stdio, redefine qh_malloc, qh_free, and qh_fprintf_rbox (user.c)
    Split out qh_rboxpoints2() to avoid -Wclobbered

  design:
    Straight line code (consider defining a struct and functions):

    Parse arguments into variables
    Determine the number of points
    Generate the points
*/
int qh_rboxpoints(FILE* fout, FILE* ferr, char* rbox_command) {
  int exitcode;
  double *simplex;

  if (rbox_inuse) {
    qh_fprintf_stderr(6188, "rbox error: rbox in use by another process.  Please lock calls to rbox or use libqhull_r/rboxlib_r.c\n");
    return qh_ERRqhull;
  }
  rbox_inuse= True;
  rbox.ferr= ferr;
  rbox.fout= fout;
  
  simplex= NULL;
  exitcode= setjmp(rbox.errexit);
  if (exitcode) {
    /* same code for error exit and normal return.  qh.NOerrexit is set */
    if (simplex)
      qh_free(simplex);
    rbox_inuse= False;
    return exitcode;
  }
  qh_rboxpoints2(rbox_command, &simplex);
  /* same code for error exit and normal return */
  if (simplex)
    qh_free(simplex);
  rbox_inuse= False;
  return qh_ERRnone;
} /* rboxpoints */

void qh_rboxpoints2(char* rbox_command, double **simplex) {
  int i,j,k;
  int gendim;
  int coincidentcount=0, coincidenttotal=0, coincidentpoints=0;
  int cubesize, diamondsize, seed=0, count, apex;
  int dim=3, numpoints=0, totpoints, addpoints=0;
  int issphere=0, isaxis=0,  iscdd=0, islens=0, isregular=0, iswidth=0, addcube=0;
  int isgap=0, isspiral=0, NOcommand=0, adddiamond=0;
  int israndom=0, istime=0;
  int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
  double width=0.0, gap=0.0, radius=0.0, coincidentradius=0.0;
  double coord[MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
  double *coordp, *simplexp;
  int nthroot, mult[MAXdim];
  double norm, factor, randr, rangap, tempr, lensangle=0, lensbase=1;
  double anglediff, angle, x, y, cube=0.0, diamond=0.0;
  double box= qh_DEFAULTbox; /* scale all numbers before output */
  double randmax= qh_RANDOMmax;
  char command[250], seedbuf[50];
  char *s=command, *t, *first_point=NULL;
  time_t timedata;

  *command= '\0';
  strncat(command, rbox_command, sizeof(command)-sizeof(seedbuf)-strlen(command)-1);

  while (*s && !isspace(*s))  /* skip program name */
    s++;
  while (*s) {
    while (*s && isspace(*s))
      s++;
    if (*s == '-')
      s++;
    if (!*s)
      break;
    if (isdigit(*s)) {
      numpoints= qh_strtol(s, &s);
      continue;
    }
    /* ============= read flags =============== */
    switch (*s++) {
    case 'c':
      addcube= 1;
      t= s;
      while (isspace(*t))
        t++;
      if (*t == 'G')
        cube= qh_strtod(++t, &s);
      break;
    case 'd':
      adddiamond= 1;
      t= s;
      while (isspace(*t))
        t++;
      if (*t == 'G')
        diamond= qh_strtod(++t, &s);
      break;
    case 'h':
      iscdd= 1;
      break;
    case 'l':
      isspiral= 1;
      break;
    case 'n':
      NOcommand= 1;
      break;
    case 'r':
      isregular= 1;
      break;
    case 's':
      issphere= 1;
      break;
    case 't':
      istime= 1;
      if (isdigit(*s)) {
        seed= qh_strtol(s, &s);
        israndom= 0;
      }else
        israndom= 1;
      break;
    case 'x':
      issimplex= 1;
      break;
    case 'y':
      issimplex2= 1;
      break;
    case 'z':
      rbox.isinteger= 1;
      break;
    case 'B':
      box= qh_strtod(s, &s);
      isbox= 1;
      break;
    case 'C':
      if (*s)
        coincidentpoints=  qh_strtol(s, &s);
      if (*s == ',') {
        ++s;
        coincidentradius=  qh_strtod(s, &s);
      }
      if (*s == ',') {
        ++s;
        coincidenttotal=  qh_strtol(s, &s);
      }
      if (*s && !isspace(*s)) {
        qh_fprintf_rbox(rbox.ferr, 7080, "rbox error: arguments for 'Cn,r,m' are not 'int', 'float', and 'int'.  Remaining string is '%s'\n", s);
        qh_errexit_rbox(qh_ERRinput);
      }
      if (coincidentpoints==0){
        qh_fprintf_rbox(rbox.ferr, 6268, "rbox error: missing arguments for 'Cn,r,m' where n is the number of coincident points, r is the radius (default 0.0), and m is the number of points\n");
        qh_errexit_rbox(qh_ERRinput);
      }
      if (coincidentpoints<0 || coincidenttotal<0 || coincidentradius<0.0){
        qh_fprintf_rbox(rbox.ferr, 6269, "rbox error: negative arguments for 'Cn,m,r' where n (%d) is the number of coincident points, m (%d) is the number of points, and r (%.2g) is the radius (default 0.0)\n", coincidentpoints, coincidenttotal, coincidentradius);
        qh_errexit_rbox(qh_ERRinput);
      }
      break;
    case 'D':
      dim= qh_strtol(s, &s);
      if (dim < 1
      || dim > MAXdim) {
        qh_fprintf_rbox(rbox.ferr, 6189, "rbox error: dimension, D%d, out of bounds (>=%d or <=0)\n", dim, MAXdim);
        qh_errexit_rbox(qh_ERRinput);
      }
      break;
    case 'G':
      if (isdigit(*s))
        gap= qh_strtod(s, &s);
      else
        gap= 0.5;
      isgap= 1;
      break;
    case 'L':
      if (isdigit(*s))
        radius= qh_strtod(s, &s);
      else
        radius= 10;
      islens= 1;
      break;
    case 'M':
      ismesh= 1;
      if (*s)
        meshn= qh_strtod(s, &s);
      if (*s == ',') {
        ++s;
        meshm= qh_strtod(s, &s);
      }else
        meshm= 0.0;
      if (*s == ',') {
        ++s;
        meshr= qh_strtod(s, &s);
      }else
        meshr= sqrt(meshn*meshn + meshm*meshm);
      if (*s && !isspace(*s)) {
        qh_fprintf_rbox(rbox.ferr, 7069, "rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
        meshn= 3.0, meshm=4.0, meshr=5.0;
      }
      break;
    case 'O':
      rbox.out_offset= qh_strtod(s, &s);
      break;
    case 'P':
      if (!first_point)
        first_point= s - 1;
      addpoints++;
      while (*s && !isspace(*s))   /* read points later */
        s++;
      break;
    case 'W':
      width= qh_strtod(s, &s);
      iswidth= 1;
      break;
    case 'Z':
      if (isdigit(*s))
        radius= qh_strtod(s, &s);
      else
        radius= 1.0;
      isaxis= 1;
      break;
    default:
      qh_fprintf_rbox(rbox.ferr, 6352, "rbox error: unknown flag at '%s'.\nExecute 'rbox' without arguments for documentation.\n", s - 1);
      qh_errexit_rbox(qh_ERRinput);
    }
    if (*s && !isspace(*s)) {
      qh_fprintf_rbox(rbox.ferr, 6353, "rbox error: missing space between flags at %s.\n", s);
      qh_errexit_rbox(qh_ERRinput);
    }
  }

  /* ============= defaults, constants, and sizes =============== */
  if (rbox.isinteger && !isbox)
    box= qh_DEFAULTzbox;
  if (addcube) {
    tempr= floor(ldexp(1.0,dim)+0.5);
    cubesize= (int)tempr;
    if (cube == 0.0)
      cube= box;
  }else
    cubesize= 0;
  if (adddiamond) {
    diamondsize= 2*dim;
    if (diamond == 0.0)
      diamond= box;
  }else
    diamondsize= 0;
  if (islens) {
    if (isaxis) {
        qh_fprintf_rbox(rbox.ferr, 6190, "rbox error: can not combine 'Ln' with 'Zn'\n");
        qh_errexit_rbox(qh_ERRinput);
    }
    if (radius <= 1.0) {
        qh_fprintf_rbox(rbox.ferr, 6191, "rbox error: lens radius %.2g should be greater than 1.0\n",
               radius);
        qh_errexit_rbox(qh_ERRinput);
    }
    lensangle= asin(1.0/radius);
    lensbase= radius * cos(lensangle);
  }

  if (!numpoints) {
    if (issimplex2)
        ; /* ok */
    else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
        qh_fprintf_rbox(rbox.ferr, 6192, "rbox error: missing count\n");
        qh_errexit_rbox(qh_ERRinput);
    }else if (adddiamond + addcube + addpoints)
        ; /* ok */
    else {
        numpoints= 50;  /* ./rbox D4 is the test case */
        issphere= 1;
    }
  }
  if ((issimplex + islens + isspiral + ismesh > 1)
  || (issimplex + issphere + isspiral + ismesh > 1)) {
    qh_fprintf_rbox(rbox.ferr, 6193, "rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
    qh_errexit_rbox(qh_ERRinput);
  }
  if (coincidentpoints>0 && (numpoints == 0 || coincidenttotal > numpoints)) {
    qh_fprintf_rbox(rbox.ferr, 6270, "rbox error: 'Cn,r,m' requested n coincident points for each of m points.  Either there is no points or m (%d) is greater than the number of points (%d).\n", coincidenttotal, numpoints);
    qh_errexit_rbox(qh_ERRinput);
  }
  if (coincidentpoints > 0 && isregular) {
    qh_fprintf_rbox(rbox.ferr, 6423, "rbox error: 'Cn,r,m' is not implemented for regular points ('r')\n");
    qh_errexit_rbox(qh_ERRinput);
  }

  if (coincidenttotal == 0)
    coincidenttotal= numpoints;

  /* ============= print header with total points =============== */
  if (issimplex || ismesh)
    totpoints= numpoints;
  else if (issimplex2)
    totpoints= numpoints+dim+1;
  else if (isregular) {
    totpoints= numpoints;
    if (dim == 2) {
        if (islens)
          totpoints += numpoints - 2;
    }else if (dim == 3) {
        if (islens)
          totpoints += 2 * numpoints;
      else if (isgap)
        totpoints += 1 + numpoints;
      else
        totpoints += 2;
    }
  }else
    totpoints= numpoints + isaxis;
  totpoints += cubesize + diamondsize + addpoints;
  totpoints += coincidentpoints*coincidenttotal;

  /* ============= seed randoms =============== */
  if (istime == 0) {
    for (s=command; *s; s++) {
      if (issimplex2 && *s == 'y') /* make 'y' same seed as 'x' */
        i= 'x';
      else
        i= *s;
      seed= 11*seed + i;
    }
  }else if (israndom) {
    seed= (int)time(&timedata);
    sprintf(seedbuf, " t%d", seed);  /* appends an extra t, not worth removing */
    strncat(command, seedbuf, sizeof(command) - strlen(command) - 1);
    t= strstr(command, " t ");
    if (t)
      strcpy(t+1, t+3); /* remove " t " */
  } /* else, seed explicitly set to n */
  qh_RANDOMseed_(seed);

  /* ============= print header =============== */

  if (iscdd)
      qh_fprintf_rbox(rbox.fout, 9391, "%s\nbegin\n        %d %d %s\n",
      NOcommand ? "" : command,
      totpoints, dim+1,
      rbox.isinteger ? "integer" : "real");
  else if (NOcommand)
      qh_fprintf_rbox(rbox.fout, 9392, "%d\n%d\n", dim, totpoints);
  else
      /* qh_fprintf_rbox special cases 9393 to append 'command' to the RboxPoints.comment() */
      qh_fprintf_rbox(rbox.fout, 9393, "%d %s\n%d\n", dim, command, totpoints);

  /* ============= explicit points =============== */
  if ((s= first_point)) {
    while (s && *s) { /* 'P' */
      count= 0;
      if (iscdd)
        qh_out1(1.0);
      while (*++s) {
        qh_out1(qh_strtod(s, &s));
        count++;
        if (isspace(*s) || !*s)
          break;
        if (*s != ',') {
          qh_fprintf_rbox(rbox.ferr, 6194, "rbox error: missing comma after coordinate in %s\n\n", s);
          qh_errexit_rbox(qh_ERRinput);
        }
      }
      if (count < dim) {
        for (k=dim-count; k--; )
          qh_out1(0.0);
      }else if (count > dim) {
        qh_fprintf_rbox(rbox.ferr, 6195, "rbox error: %d coordinates instead of %d coordinates in %s\n\n",
                  count, dim, s);
        qh_errexit_rbox(qh_ERRinput);
      }
      qh_fprintf_rbox(rbox.fout, 9394, "\n");
      while ((s= strchr(s, 'P'))) {
        if (isspace(s[-1]))
          break;
      }
    }
  }

  /* ============= simplex distribution =============== */
  if (issimplex+issimplex2) {
    if (!(*simplex= (double *)qh_malloc( (size_t)(dim * (dim+1)) * sizeof(double)))) {
      qh_fprintf_rbox(rbox.ferr, 6196, "rbox error: insufficient memory for simplex\n");
      qh_errexit_rbox(qh_ERRmem); /* qh_ERRmem */
    }
    simplexp= *simplex;
    if (isregular) {
      for (i=0; i randmax/2)
          coord[dim-1]= -coord[dim-1];
      /* ============= project 'Wn' point toward boundary =============== */
      }else if (iswidth && !issphere) {
        j= qh_RANDOMint % gendim;
        if (coord[j] < 0)
          coord[j]= -1.0 - coord[j] * width;
        else
          coord[j]= 1.0 - coord[j] * width;
      }
      /* ============= scale point to box =============== */
      for (k=0; k=0; k--) {
        if (j & ( 1 << k))
          qh_out1(cube);
        else
          qh_out1(-cube);
      }
      qh_fprintf_rbox(rbox.fout, 9400, "\n");
    }
  }

  /* ============= write diamond vertices =============== */
  if (adddiamond) {
    for (j=0; j=0; k--) {
        if (j/2 != k)
          qh_out1(0.0);
        else if (j & 0x1)
          qh_out1(diamond);
        else
          qh_out1(-diamond);
      }
      qh_fprintf_rbox(rbox.fout, 9401, "\n");
    }
  }

  if (iscdd)
    qh_fprintf_rbox(rbox.fout, 9402, "end\nhull\n");
} /* rboxpoints2 */

/*------------------------------------------------
outxxx - output functions for qh_rboxpoints
*/
int qh_roundi(double a) {
  if (a < 0.0) {
    if (a - 0.5 < INT_MIN) {
      qh_fprintf_rbox(rbox.ferr, 6200, "rbox input error: negative coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
      qh_errexit_rbox(qh_ERRinput);
    }
    return (int)(a - 0.5);
  }else {
    if (a + 0.5 > INT_MAX) {
      qh_fprintf_rbox(rbox.ferr, 6201, "rbox input error: coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
      qh_errexit_rbox(qh_ERRinput);
    }
    return (int)(a + 0.5);
  }
} /* qh_roundi */

void qh_out1(double a) {

  if (rbox.isinteger)
    qh_fprintf_rbox(rbox.fout, 9403, "%d ", qh_roundi(a+rbox.out_offset));
  else
    qh_fprintf_rbox(rbox.fout, 9404, qh_REAL_1, a+rbox.out_offset);
} /* qh_out1 */

void qh_out2n(double a, double b) {

  if (rbox.isinteger)
    qh_fprintf_rbox(rbox.fout, 9405, "%d %d\n", qh_roundi(a+rbox.out_offset), qh_roundi(b+rbox.out_offset));
  else
    qh_fprintf_rbox(rbox.fout, 9406, qh_REAL_2n, a+rbox.out_offset, b+rbox.out_offset);
} /* qh_out2n */

void qh_out3n(double a, double b, double c) {

  if (rbox.isinteger)
    qh_fprintf_rbox(rbox.fout, 9407, "%d %d %d\n", qh_roundi(a+rbox.out_offset), qh_roundi(b+rbox.out_offset), qh_roundi(c+rbox.out_offset));
  else
    qh_fprintf_rbox(rbox.fout, 9408, qh_REAL_3n, a+rbox.out_offset, b+rbox.out_offset, c+rbox.out_offset);
} /* qh_out3n */

void qh_outcoord(int iscdd, double *coord, int dim) {
    double *p= coord;
    int k;

    if (iscdd)
      qh_out1(1.0);
    for (k=0; k < dim; k++)
      qh_out1(*(p++));
    qh_fprintf_rbox(rbox.fout, 9396, "\n");
} /* qh_outcoord */

void qh_outcoincident(int coincidentpoints, double radius, int iscdd, double *coord, int dim) {
  double *p;
  double randr, delta;
  int i,k;
  double randmax= qh_RANDOMmax;

  for (i=0; i