Replace rand/random/lrand48/etc with lib/gis random number functions

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@61415 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/defs/gmath.h

@@ -22,14 +22,16 @@ extern int fft(int, double *[2], int, int, int);
 extern int fft2(int, double (*)[2], int, int, int);
 
 /* gauss.c */
-extern double G_math_rand_gauss(int, double);
+extern double G_math_rand_gauss(double);
 
 /* max_pow2.c */
 extern long G_math_max_pow2 (long n);
 extern long G_math_min_pow2 (long n);
 
 /* rand1.c */
-extern float G_math_rand(int);
+extern void G_math_srand(int);
+extern int G_math_srand_auto(void);
+extern float G_math_rand(void);
 
 /* del2g.c */
 extern int del2g(double *[2], int, double);

lib/gmath/gauss.c

@@ -14,14 +14,14 @@
  * \return double
  */
 
-double G_math_rand_gauss(const int seed, const double sigma)
+double G_math_rand_gauss(double sigma)
 {
     double x, y, r2;
 
     do {
 	/* choose x,y in uniform square (-1,-1) to (+1,+1) */
-	x = -1 + 2 * G_math_rand(seed);
-	y = -1 + 2 * G_math_rand(seed);
+	x = -1 + 2 * G_math_rand();
+	y = -1 + 2 * G_math_rand();
 
 	/* see if it is in the unit circle */
 	r2 = x * x + y * y;

lib/gmath/rand1.c

@@ -1,5 +1,5 @@
-#include <grass/config.h>
 #include <stdlib.h>
+#include <grass/gis.h>
 #include <grass/gmath.h>
 
 
@@ -14,17 +14,29 @@
  * \return float
  */
 
-float G_math_rand(int seed)
+float G_math_rand(void)
 {
-#if defined(HAVE_DRAND48)
-    if (seed < 0)
-	srand48(-seed);
+    return G_drand48();
+}
+
+/*!
+ * \brief Seed the pseudo-random number generator
+ *
+ * \param seedval 32-bit integer used to seed the PRNG
+ */
+
+void G_math_srand(int seed)
+{
+    G_srand48(seed);
+}
 
-    return (float)drand48();
-#else
-    if (seed < 0)
-	srand(-seed);
+/*!
+ * \brief Seed the pseudo-random number generator from the time and PID
+ *
+ * \return generated seed value passed to G_srand48()
+ */
 
-    return 1.0f * rand() / RAND_MAX;
-#endif
+int G_math_srand_auto(void)
+{
+    return (int) G_srand48_auto();
 }

lib/raster/color_rand.c

@@ -33,16 +33,17 @@ void Rast_make_random_colors(struct Colors *colors, CELL min, CELL max)
 	G_fatal_error(_("Rast_make_random_colors: min (%d) > max (%d)"),
 		      min, max);
 
-    srand(time(NULL));
+    /* FIXME - allow seed to be specified for repeatability */
+    G_srand48_auto();
 
-    count = MAX_COLORS - DEVIATION + rand() % DEVIATION;
+    count = MAX_COLORS - DEVIATION + G_lrand48() % DEVIATION;
     if (count > max - min + 1)
 	count = max - min + 1;
 
     for (n = 1; n <= count; n++) {
-	red = rand() & 0377;
-	grn = rand() & 0377;
-	blu = rand() & 0377;
+	red = G_lrand48() & 0xff;
+	grn = G_lrand48() & 0xff;
+	blu = G_lrand48() & 0xff;
 	Rast_add_modular_c_color_rule(&n, red, grn, blu,
 				      &n, red, grn, blu, colors);
     }

raster/r.flow/main.c

@@ -26,8 +26,6 @@
  **  Update MN: commented line 387
  */
 
-#include <stdlib.h>		/* for the random number generation */
-#include <time.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
 #include <grass/glocale.h>
@@ -320,7 +318,8 @@ static void calculate(void)
 
     ystep = region.ns_res * (double)loopstep;
 
-    srand(time(0));
+    /* FIXME - allow seed to be specified for repeatability */
+    G_srand48_auto();
 
     for (row = 0, y = (double)region.north - (region.ns_res * .5);
 	 row < region.rows; row += loopstep, y -= ystep) {
@@ -337,9 +336,9 @@ static void calculate(void)
 #ifdef OFFSET
 		/* disabled by helena June 2005 */
 		roffset = parm.offset * (double)region.ew_res
-		    * ((2. * (double)rand() / (double)RAND_MAX) - 1.);
+		    * ((2. * G_drand48()) - 1.);
 		coffset = parm.offset * (double)region.ns_res
-		    * ((2. * (double)rand() / (double)RAND_MAX) - 1.);
+		    * ((2. * G_drand48()) - 1.);
 #endif
 		pts.x = x;
 		pts.y = y;

raster/r.li/r.li.daemon/daemon.c

@@ -15,7 +15,6 @@
  * \include
  * 
  */
-#include <grass/config.h>
 #include <stdlib.h>
 #include <stddef.h>
 #include <fcntl.h>
@@ -36,11 +35,6 @@
 #include <grass/glocale.h>
 #include "daemon.h"
 
-#ifdef __MINGW32__
-#define srandom srand
-#define random rand
-#endif
-
 int calculateIndex(char *file, rli_func *f,
 		   char **parameters, char *raster, char *output)
 {
@@ -524,11 +518,11 @@ int disposeAreas(struct list *l, struct g_area *g, char *def)
 	    G_fatal_error(_("Too many units to place"));
 	assigned = G_malloc(units * sizeof(int));
 	i = 0;
-	srandom(0);
+	G_srand48(0);
 	while (i < units) {
 	    int j, position, found = FALSE;
 
-	    position = random() % max_units;
+	    position = G_lrand48() % max_units;
 	    for (j = 0; j < i; j++) {
 		if (assigned[j] == position)
 		    found = TRUE;
@@ -592,7 +586,7 @@ int disposeAreas(struct list *l, struct g_area *g, char *def)
 	if (r_strat_len < g->rl || c_strat_len < g->cl)
 	    G_fatal_error(_("Too many stratified random sample for raster map"));
 	loop = r_strat * c_strat;
-	srandom(0);
+	G_srand48(0);
 	for (i = 0; i < loop; i++) {
 	    msg m;
 
@@ -600,9 +594,9 @@ int disposeAreas(struct list *l, struct g_area *g, char *def)
 		m.type = AREA;
 		m.f.f_a.aid = i;
 		m.f.f_a.x = (int)g->sf_x + ((i % c_strat) * c_strat_len) +
-		    (random() % (c_strat_len - g->cl));
+		    (G_lrand48() % (c_strat_len - g->cl));
 		m.f.f_a.y = (int)g->sf_y + (rint(i / c_strat) * r_strat_len) +
-		    (random() % (r_strat_len - g->rl));
+		    (G_lrand48() % (r_strat_len - g->rl));
 		m.f.f_a.rl = g->rl;
 		m.f.f_a.cl = g->cl;
 		insertNode(l, m);
@@ -611,10 +605,10 @@ int disposeAreas(struct list *l, struct g_area *g, char *def)
 		m.type = MASKEDAREA;
 		m.f.f_ma.aid = i;
 		m.f.f_ma.x = (int)g->sf_x + ((i % c_strat) * c_strat_len) +
-		    (random() % (c_strat_len - g->cl));
+		    (G_lrand48() % (c_strat_len - g->cl));
 		m.f.f_ma.y =
 		    (int)g->sf_y + (rint(i / c_strat) * r_strat_len) +
-		    (random() % (r_strat_len - g->rl));
+		    (G_lrand48() % (r_strat_len - g->rl));
 		m.f.f_ma.rl = g->rl;
 		m.f.f_ma.cl = g->cl;
 		strcpy(m.f.f_ma.mask, g->maskname);

raster/r.random/creat_rand.c

@@ -1,43 +1,14 @@
-#include <time.h>
-#include <stdlib.h>
-#include <unistd.h>
-#include <sys/types.h>
-#include <grass/config.h>
+#include <grass/gis.h>
 
 
-#ifndef USE_RAND
-
-#ifndef HAVE_DRAND48
-#define lrand48() (((long) rand() ^ ((long) rand() << 16)) & 0x7FFFFFFF)
-#define srand48(sv) (srand((unsigned)(sv)))
-#endif
-
-
-long make_rand(void)
-{
-    return lrand48();
-}
-
-void init_rand(void)
-{
-    srand48((long)time((time_t *) 0));
-}
-
-#else
-
-static long labs(int n)
-{
-    return n < 0 ? (-n) : n;
-}
-
 long make_rand(void)
 {
-    return (labs(rand() + (rand() << 16)));
+    return G_lrand48();
 }
 
 void init_rand(void)
 {
-    srand(getpid());
+    /* FIXME - allow seed to be specified for repeatability */
+    G_srand48_auto();
 }
 
-#endif

raster/r.spread/main.c

@@ -336,7 +336,8 @@ int main(int argc, char *argv[])
 	exit(EXIT_FAILURE);
 
 
-    srand(getpid());
+    /* FIXME - allow seed to be specified for repeatability */
+    G_srand48_auto();
 
     display = flag.display->answer;
 #if 1

raster/r.spread/pick_dist.c

@@ -18,13 +18,14 @@
 #include <stdlib.h>
 #include <limits.h>
 #include <math.h>
+#include <grass/gis.h>
 #include "local_proto.h"
 
 int pick_dist(int u)
 {
     int v;
 
-    v = (int)((u + 0.99999999999) * rand() / INT_MAX);
-    u = (int)((v + 0.99999999999) * rand() / INT_MAX);
-    return ((int)((u + 0.99999999999) * rand() / INT_MAX));	/*4th for a test */
+    v = (int)((u + 0.99999999999) * G_drand48());
+    u = (int)((v + 0.99999999999) * G_drand48());
+    return ((int)((u + 0.99999999999) * G_drand48()));	/*4th for a test */
 }

raster/r.spread/pick_ignite.c

@@ -13,6 +13,7 @@
  ******************************************************************************/
 #include <stdlib.h>
 #include <limits.h>
+#include <grass/gis.h>
 #include "local_proto.h"
 
 int p[18] = { 100, 100, 100, 90, 80, 70, 60, 50, 40,	/*adapt the "average" */
@@ -21,5 +22,5 @@ int p[18] = { 100, 100, 100, 90, 80, 70, 60, 50, 40,	/*adapt the "average" */
 
 int pick_ignite(int i)
 {
-    return ((100.0 * rand() / INT_MAX) <= p[i]);
+    return ((100.0 * G_drand48()) <= p[i]);
 }

raster/r.surf.fractal/spec_syn.c

@@ -38,14 +38,13 @@ int specsyn(double *data[2],	/* Array holding complex data to transform. */
 
     int row, col,		/* Counts through half the data array.  */
       row0, col0,		/* 'other half' of the array.           */
-      coeff,			/* No. of Fourier coeffficents to calc. */
-      seed;			/* Random number seed.                  */
+      coeff;			/* No. of Fourier coeffficents to calc. */
 
     double phase, rad,		/* polar coordinates of Fourier coeff.  */
      *temp[2];
 
-    seed = -1 * getpid();
-    G_math_rand(seed);		/* Reset random number generator.       */
+    /* FIXME - allow seed to be specified for repeatability */
+    G_math_srand_auto();	/* Reset random number generator.       */
 
     temp[0] = (double *)G_malloc(nn * nn * sizeof(double));
     temp[1] = (double *)G_malloc(nn * nn * sizeof(double));
@@ -65,12 +64,12 @@ int specsyn(double *data[2],	/* Array holding complex data to transform. */
 	for (col = 0; col <= nn / 2; col++) {
 	    /* Generate random Fourier coefficients. */
 
-	    phase = TWOPI * G_math_rand(0);
+	    phase = TWOPI * G_math_rand();
 
 	    if ((row != 0) || (col != 0))
 		rad =
 		    pow(row * row + col * col,
-			-(H + 1) / 2.0) * G_math_rand_gauss(0, 1.);
+			-(H + 1) / 2.0) * G_math_rand_gauss(1.);
 	    else
 		rad = 0.0;
 
@@ -101,10 +100,10 @@ int specsyn(double *data[2],	/* Array holding complex data to transform. */
 
     for (row = 1; row < nn / 2; row++)
 	for (col = 1; col < nn / 2; col++) {
-	    phase = TWOPI * G_math_rand(0);
+	    phase = TWOPI * G_math_rand();
 	    rad =
 		pow(row * row + col * col,
-		    -(H + 1) / 2.0) * G_math_rand_gauss(0, 1.);
+		    -(H + 1) / 2.0) * G_math_rand_gauss(1.);
 
 	    *(data[0] + row * nn + nn - col) = rad * cos(phase);
 	    *(data[1] + row * nn + nn - col) = rad * sin(phase);

raster/r.surf.gauss/gaussurf.c

@@ -32,7 +32,8 @@ int gaussurf(char *out,		/* Name of raster maps to be opened.    */
 
 	/****** INITIALISE RANDOM NUMBER GENERATOR ******/
 
-    G_math_rand(-1 * getpid());
+    /* FIXME - allow seed to be specified for repeatability */
+    G_math_srand_auto();
 
 	/****** OPEN CELL FILES AND GET CELL DETAILS ******/
 
@@ -50,7 +51,7 @@ int gaussurf(char *out,		/* Name of raster maps to be opened.    */
 	G_percent(row_count, nrows, 5);
 	for (col_count = 0; col_count < ncols; col_count++)
 	    *(row_out + col_count) =
-		(DCELL) (G_math_rand_gauss(2742, sigma) + mean);
+		(DCELL) (G_math_rand_gauss(sigma) + mean);
 
 	/* Write contents row by row */
 	Rast_put_d_row(fd_out, (DCELL *) row_out);

raster/r.surf.random/randsurf.c

@@ -1,4 +1,3 @@
-#include <unistd.h>
 #include <math.h>
 #include <grass/gis.h>
 #include <grass/raster.h>
@@ -21,7 +20,8 @@ int randsurf(char *out,		/* Name of raster maps to be opened.    */
     int row_count, col_count;
 
 	/****** INITIALISE RANDOM NUMBER GENERATOR ******/
-    G_math_rand(-1 * getpid());
+    /* FIXME - allow seed to be specified for repeatability */
+    G_math_srand_auto();
 
 	/****** OPEN CELL FILES AND GET CELL DETAILS ******/
     fd_out = Rast_open_new(out, int_map ? CELL_TYPE : DCELL_TYPE);
@@ -40,12 +40,12 @@ int randsurf(char *out,		/* Name of raster maps to be opened.    */
 	for (col_count = 0; col_count < ncols; col_count++) {
 	    if (int_map)
 		*(row_out_C + col_count) =
-		    (CELL) (G_math_rand(2742) * (max + 1 - min) + min);
+		    (CELL) (G_math_rand() * (max + 1 - min) + min);
 	    /* under represents first and last bin */
 	    /*                  *(row_out_C + col_count) = (CELL) floor(rand1(2742)*(max-min)+min +0.5); */
 	    else
 		*(row_out_D + col_count) =
-		    (DCELL) (G_math_rand(2742) * (max - min) + min);
+		    (DCELL) (G_math_rand() * (max - min) + min);
 	}
 
 	/* Write contents row by row */

vector/v.extract/main.c

@@ -317,19 +317,15 @@ int main(int argc, char **argv)
 	    G_fatal_error(_("Random category count is larger or equal to "
 			    "uniq <%s> feature category count %d"),
 			  opt.type->answer, ucat_count);
-
-	if (ucat_count >= RAND_MAX)
-	    G_fatal_error(_("There are more categories than random number "
-			    "generator can reach. Report this as a GRASS bug."));
 	
-	seed = getpid();
 	/* Initialise random number generator */
-	G_math_rand(-1 * seed);
+	/* FIXME - allow seed to be specified for repeatability */
+	G_math_srand_auto();
 
 	/* Fill cat_array with list of valid random numbers */
 	while (cat_count < nrandom) {
 	    /* Random number in range from 0 to largest CAT value */
-	    prnd = (int) floor(G_math_rand(seed) *
+	    prnd = (int) floor(G_math_rand() *
 			       (ucat_array[ucat_count - 1] + 1));
 	    qsort(cat_array, cat_count, sizeof(int), cmp);
 	    /* Check if generated number isn't already in 

vector/v.kcv/main.c

@@ -19,25 +19,19 @@
 
 
 #include <stdlib.h>
-#include <unistd.h>
 #include <math.h>
 #include <grass/gis.h>
 #include <grass/dbmi.h>
 #include <grass/vector.h>
 #include <grass/glocale.h>
 
-static double myrand(void)
-{
-    return rand() / (1.0 + RAND_MAX);
-}
-
 struct Cell_head window;
 
 int main(int argc, char *argv[]) 
 {
     int i, line, nlines, nlinks;
     int idx, *line_idx, lines_left;
-    double (*rng) ();
+    double (*rng)(void) = G_drand48;
     int nsites, p, np, min_count, spill;
     struct partition {
 	int id, count, max;
@@ -48,7 +42,6 @@ int main(int argc, char *argv[])
     struct line_cats *Cats;
     struct GModule *module;
     struct Option *map_opt, *col_opt, *npart_opt, *field_opt;
-    struct Flag *drand48_flag;
  
     /* Attributes */
     struct field_info *Fi;
@@ -81,15 +74,6 @@ int main(int argc, char *argv[])
     col_opt->description =
 	_("Name for new column to which partition number is written");
 
-    drand48_flag = G_define_flag();
-    drand48_flag->key = 'd';
-#ifdef HAVE_DRAND48
-    drand48_flag->description = _("Use drand48()");
-#else
-    drand48_flag->description = _("Use drand48() (ignored)");
-#endif
-
-
     G_gisinit(argv[0]);
 
     if (G_parser(argc, argv))
@@ -99,17 +83,8 @@ int main(int argc, char *argv[])
     if (np < 2)
 	G_fatal_error(_("'%s' must be > 1"), npart_opt->key);
 
-#ifdef HAVE_DRAND48
-    if (drand48_flag->answer) {
-	rng = drand48;
-	srand48((long)getpid());
-    }
-    else
-#endif
-    {
-	rng = myrand;
-	srand(getpid());
-    }
+    /* FIXME - allow seed to be specified for repeatability */
+    G_srand48_auto();
 
     Points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();

vector/v.qcount/findquads.c

@@ -6,26 +6,13 @@
  * Read the file GPL.TXT coming with GRASS for details.
  */
 
-#include <sys/types.h>		/* for getpid() */
-#include <unistd.h>
 #include <stdlib.h>
 #include <math.h>
 #include <grass/gis.h>
 #include <grass/vector.h>
 #include "quaddefs.h"
 
-#ifndef RAND_MAX
-#define RAND_MAX (pow(2.0,31.0)-1)
-#endif
-
-#define R_INIT srand
-#define RANDOM(lo,hi) ((double) rand() / (double)RAND_MAX * (hi-lo) + lo)
-
-/*
-   #define R_MAX 1.0
-   #define R_INIT srand48
-   #define RANDOM(lo,hi) (drand48()/R_MAX*(hi-lo)+lo)
- */
+#define RANDOM(lo,hi) (G_drand48() * (hi-lo) + lo)
 
 /*
  * returns Z struct filled with centers of n non-overlapping circles of
@@ -42,8 +29,8 @@ COOR *find_quadrats(int n, double r, struct Cell_head window)
     if (quads == NULL)
 	G_fatal_error("cannot allocate memory for quadrats");
 
-    srand((unsigned int)getpid());
-    /* R_INIT (1); */
+    /* FIXME - allow seed to be specified for repeatability */
+    G_srand48_auto();
 
     e_max = window.east - r;
     e_min = window.west + r;

vector/v.random/main.c

@@ -32,8 +32,6 @@
 
 #include <stdlib.h>
 #include <math.h>
-#include <sys/types.h>
-#include <unistd.h>
 #include <string.h>
 
 #include <grass/gis.h>
@@ -41,20 +39,6 @@
 #include <grass/dbmi.h>
 #include <grass/glocale.h>
 
-#ifndef RAND_MAX
-#define RAND_MAX (pow(2.0,31.0)-1)
-#endif
-double myrand(void);
-
-#if defined(__CYGWIN__) || defined(__APPLE__) || defined(__MINGW32__)
-double drand48()
-{
-    return (rand() / (1.0 * RAND_MAX));
-}
-
-#define srand48(sv) (srand((unsigned)(sv)))
-#endif
-
 /* for qsort */
 
 typedef struct {
@@ -77,10 +61,10 @@ static int sort_by_size(const void *a, const void *b)
 int main(int argc, char *argv[])
 {
     char *output, buf[DB_SQL_MAX];
-    double (*rng) ();
-    double max, zmin, zmax;
+    double (*rng)(void) = G_drand48;
+    double zmin, zmax;
     int seed;
-    int i, j, k, n, b, type, usefloat;
+    int i, j, k, n, type, usefloat;
     int area, nareas, field;
     struct boxlist *List = NULL;
     BOX_SIZE *size_list = NULL;
@@ -99,7 +83,7 @@ int main(int argc, char *argv[])
     } parm;
     struct
     {
-	struct Flag *rand, *drand48, *z, *notopo, *a;
+	struct Flag *z, *notopo, *a;
     } flag;
     struct field_info *Fi;
     dbDriver *driver;
@@ -187,10 +171,6 @@ int main(int argc, char *argv[])
     flag.a->key = 'a';
     flag.a->description = _("Generate n points for each individual area");
 
-    flag.drand48 = G_define_flag();
-    flag.drand48->key = 'd';
-    flag.drand48->description = _("Use drand48() function instead of rand()");
-
     flag.notopo = G_define_standard_flag(G_FLG_V_TOPO);
 
     if (G_parser(argc, argv))
@@ -198,7 +178,6 @@ int main(int argc, char *argv[])
 
     output = parm.output->answer;
     n = atoi(parm.nsites->answer);
-    b = flag.drand48->answer ? TRUE : FALSE;
     
     if(parm.seed->answer)
         seed = atoi(parm.seed->answer);
@@ -304,24 +283,11 @@ int main(int argc, char *argv[])
 
     Vect_hist_command(&Out);
 
-    if (b) {
-	rng = drand48;
-	max = 1.0;
-	/* Init the random seed */
-	if(parm.seed->answer)
-	    srand48((long)seed);
-	else
-	    srand48((long)getpid());
-    }
-    else {  	/* default is rand() */
-	rng = myrand;
-	max = RAND_MAX;
-	/* Init the random seed */
-	if(parm.seed->answer)
-	    srand(seed);
-	else
-	    srand(getpid());
-    }
+    /* Init the random seed */
+    if(parm.seed->answer)
+	G_srand48(seed);
+    else
+	G_srand48_auto();
 
     G_get_window(&window);
 
@@ -462,9 +428,9 @@ int main(int argc, char *argv[])
 		Vect_reset_cats(Cats);
 
 		while (outside) {
-		    x = rng() / max * (bbox.W - bbox.E) + bbox.E;
-		    y = rng() / max * (bbox.N - bbox.S) + bbox.S;
-		    z = rng() / max * (zmax - zmin) + zmin;
+		    x = rng() * (bbox.W - bbox.E) + bbox.E;
+		    y = rng() * (bbox.N - bbox.S) + bbox.S;
+		    z = rng() * (zmax - zmin) + zmin;
 
 		    ret = Vect_point_in_area(x, y, &In, area, &abox);
 
@@ -512,9 +478,9 @@ int main(int argc, char *argv[])
 	    Vect_reset_line(Points);
 	    Vect_reset_cats(Cats);
 
-	    x = rng() / max * (window.west - window.east) + window.east;
-	    y = rng() / max * (window.north - window.south) + window.south;
-	    z = rng() / max * (zmax - zmin) + zmin;
+	    x = rng() * (window.west - window.east) + window.east;
+	    y = rng() * (window.north - window.south) + window.south;
+	    z = rng() * (zmax - zmin) + zmin;
 	    
 	    if (nareas) {
 		int outside = 1;
@@ -588,9 +554,9 @@ int main(int argc, char *argv[])
 			}
 		    }
 		    if (outside) {
-			x = rng() / max * (window.west - window.east) + window.east;
-			y = rng() / max * (window.north - window.south) + window.south;
-			z = rng() / max * (zmax - zmin) + zmin;
+			x = rng() * (window.west - window.east) + window.east;
+			y = rng() * (window.north - window.south) + window.south;
+			z = rng() * (zmax - zmin) + zmin;
 		    }
 		} while (outside);
 	    }
@@ -635,8 +601,3 @@ int main(int argc, char *argv[])
 
     exit(EXIT_SUCCESS);
 }
-
-double myrand()
-{
-    return (double)rand();
-}