From: Ben Pfaff Date: Tue, 16 Nov 2004 08:08:00 +0000 (+0000) Subject: Adopt GSL random number generators, paving the way for providing the X-Git-Tag: v0.4.0~227 X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d5e5cf282cc56899913654a62c425e584acecfb2;p=pspp-builds.git Adopt GSL random number generators, paving the way for providing the complete suite of random number generators on expressions. --- diff --git a/src/ChangeLog b/src/ChangeLog index 5d353632..88dd8b6e 100644 --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,35 @@ +Mon Nov 15 23:47:40 2004 Ben Pfaff + + Adopt GSL random number generators, paving the way for providing + the complete suite of random number generators on expressions. + + * Makefile.am: Remove random.c, random.h. + + * random.c: Removed. + + * random.h: Removed. + + * algorithm.c: (algo_default_random) Use GSL functions. + + * casefile.c: (test_casefile) Use GSL RNG functions. + + * expr-evl.c: (expr_evaluate) Use GSL RNG functions for OP_NORMAL, + OP_UNIFORM. + + * sample.c: (cmd_sample) Use GSL RNG functions. + (sample_trns_proc) Ditto. + + * set.q: (static var set_seed) Removed. + (static var seed_flag) Removed. + (static var rng) New variable. + (aux_stc_custom_seed) No seed value anymore, don't print anything. + (stc_custom_seed) Use new seed functions. + (seed_is_set) Removed. + (get_rng) New function that composes the entire external + interface. + (set_rng) New function. + (random_seed) New function. + Mon Nov 15 22:08:25 2004 Ben Pfaff * expr-evl.c: (expr_evaluate) Fix XDATE.JDAY formula. Thanks to diff --git a/src/Makefile.am b/src/Makefile.am index 9176897d..af860861 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -67,7 +67,7 @@ matrix-data.c mis-val.c misc.c misc.h modify-vars.c \ moments.c moments.h numeric.c output.c output.h permissions.c \ pfm-read.c pfm-read.h \ pfm-write.c pfm-write.h \ -pool.c pool.h postscript.c print.c random.c random.h recode.c \ +pool.c pool.h postscript.c print.c recode.c \ rename-vars.c repeat.c repeat.h sample.c sel-if.c settings.h \ sfm-read.c sfm-read.h sfm-write.c sfm-write.h sfmP.h som.c som.h \ sort.c sort.h \ diff --git a/src/algorithm.c b/src/algorithm.c index 230b1437..e6b9132b 100644 --- a/src/algorithm.c +++ b/src/algorithm.c @@ -91,11 +91,12 @@ #include #include "algorithm.h" +#include #include #include #include #include "alloc.h" -#include "random.h" +#include "settings.h" /* Some of the assertions in this file are very expensive. We don't use them by default. */ @@ -310,7 +311,8 @@ is_partitioned (const void *array, size_t count, size_t size, unsigned algo_default_random (unsigned max, void *aux UNUSED) { - return rng_get_unsigned (pspp_rng ()) % max; + unsigned long r_min = gsl_rng_min (get_rng ()); + return (gsl_rng_get (get_rng ()) - r_min) % max; } /* Randomly reorders ARRAY, which contains COUNT elements of SIZE diff --git a/src/casefile.c b/src/casefile.c index 03828b71..6d74e8e2 100644 --- a/src/casefile.c +++ b/src/casefile.c @@ -739,9 +739,9 @@ exit_handler (void) casefile_destroy (casefiles); } +#include #include #include "command.h" -#include "random.h" #include "lexer.h" static void test_casefile (int pattern, size_t value_cnt, size_t case_cnt); @@ -800,11 +800,10 @@ test_casefile (int pattern, size_t value_cnt, size_t case_cnt) struct casefile *cf; struct casereader *r1, *r2; struct ccase c; - struct rng *rng; + gsl_rng *rng; size_t i, j; - rng = rng_create (); - rng_seed (rng, &zero, sizeof zero); + rng = gsl_rng_alloc (gsl_rng_mt19937); cf = casefile_create (value_cnt); for (i = 0; i < case_cnt; i++) write_random_case (cf, i); @@ -831,7 +830,7 @@ test_casefile (int pattern, size_t value_cnt, size_t case_cnt) for (i = j = 0; i < case_cnt; i++) { read_and_verify_random_case (cf, r1, i); - if (rng_get_int (rng) % pattern == 0) + if (gsl_rng_get (rng) % pattern == 0) read_and_verify_random_case (cf, r2, j++); if (i == case_cnt / 2) casefile_to_disk (cf); @@ -871,7 +870,7 @@ test_casefile (int pattern, size_t value_cnt, size_t case_cnt) casereader_destroy (r1); } casefile_destroy (cf); - rng_destroy (rng); + gsl_rng_free (rng); } static void diff --git a/src/expr-evl.c b/src/expr-evl.c index 8628cfd3..d3e02ea1 100644 --- a/src/expr-evl.c +++ b/src/expr-evl.c @@ -34,6 +34,7 @@ #include "expr.h" #include "exprP.h" #include "error.h" +#include #include #include #include @@ -46,7 +47,7 @@ #include "misc.h" #include "moments.h" #include "pool.h" -#include "random.h" +#include "settings.h" #include "str.h" #include "var.h" #include "vfm.h" @@ -1122,11 +1123,11 @@ expr_evaluate (const struct expression *e, const struct ccase *c, int case_idx, break; case OP_NORMAL: if (sp->f != SYSMIS) - sp->f *= rng_get_double_normal (pspp_rng ()); + sp->f = gsl_ran_gaussian (get_rng (), sp->f); break; case OP_UNIFORM: if (sp->f != SYSMIS) - sp->f *= rng_get_double (pspp_rng ()); + sp->f *= gsl_rng_uniform (get_rng ()); break; case OP_SYSMIS: sp->f = sp->f == SYSMIS || !finite (sp->f); diff --git a/src/random.c b/src/random.c deleted file mode 100644 index e6db94c1..00000000 --- a/src/random.c +++ /dev/null @@ -1,247 +0,0 @@ -/* PSPP - computes sample statistics. - Copyright (C) 1997-9, 2000 Free Software Foundation, Inc. - Written by Ben Pfaff . - - 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, write to the Free Software - Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA - 02111-1307, USA. */ - -#include -#include "random.h" -#include "error.h" -#include -#include -#include -#include -#include -#include "alloc.h" -#include "magic.h" -#include "settings.h" - -/* Random number generator. */ -struct rng - { - /* RC4-based random bytes. */ - uint8_t s[256]; - uint8_t i, j; - - /* Normal distribution. */ - double next_normal; - }; - - -/* Return a `random' seed by using the real time clock */ -unsigned long -random_seed(void) -{ - time_t t; - - time(&t); - - return (unsigned long) t; -} - -/* Creates a new random number generator, seeds it based on - the current time, and returns it. */ -struct rng * -rng_create (void) -{ - struct rng *rng; - static unsigned long seed=0; - unsigned long s; - - rng = xmalloc (sizeof *rng); - - - if ( seed_is_set(&s) ) - { - seed = s; - } - else if ( seed == 0 ) - { - seed = random_seed(); - } - assert(seed); - /* - if (t == 0 || set_seed_used) - { - if (set_seed == NOT_LONG) - time (&t); - else - t = set_seed; - set_seed_used=0; - } - else - t++; - */ - rng_seed (rng, &seed, sizeof seed); - rng->next_normal = NOT_DOUBLE; - return rng; -} - -/* Destroys RNG. */ -void -rng_destroy (struct rng *rng) -{ - free (rng); -} - -/* Swap bytes. */ -static void -swap_byte (uint8_t *a, uint8_t *b) -{ - uint8_t t = *a; - *a = *b; - *b = t; -} - -/* Seeds RNG based on the SIZE bytes in BUF. - At most the first 256 bytes of BUF are used. */ -void -rng_seed (struct rng *rng, const void *key_, size_t size) -{ - const uint8_t *key = key_; - size_t key_idx; - uint8_t *s; - int i, j; - - assert (rng != NULL); - - s = rng->s; - rng->i = rng->j = 0; - for (i = 0; i < 256; i++) - s[i] = i; - for (key_idx = 0, i = j = 0; i < 256; i++) - { - j = (j + s[i] + key[key_idx]) & 255; - swap_byte (s + i, s + j); - if (++key_idx >= size) - key_idx = 0; - } -} - -/* Reads SIZE random bytes from RNG into BUF. */ -void -rng_get_bytes (struct rng *rng, void *buf_, size_t size) -{ - uint8_t *buf = buf_; - uint8_t *s; - uint8_t i, j; - - assert (rng != 0); - - s = rng->s; - i = rng->i; - j = rng->j; - while (size-- > 0) - { - i += 1; - j += s[i]; - swap_byte (s + i, s + j); - *buf++ = s[(s[i] + s[j]) & 255]; - } - rng->i = i; - rng->j = j; -} - -/* Returns a random int in the range [0, INT_MAX]. */ -int -rng_get_int (struct rng *rng) -{ - int value; - - do - { - rng_get_bytes (rng, &value, sizeof value); - value = abs (value); - } - while (value < 0); - - return value; -} - -/* Returns a random unsigned in the range [0, UINT_MAX]. */ -unsigned -rng_get_unsigned (struct rng *rng) -{ - unsigned value; - - rng_get_bytes (rng, &value, sizeof value); - return value; -} - -/* Returns a random number from the uniform distribution with - range [0,1). */ -double -rng_get_double (struct rng *rng) -{ - for (;;) - { - unsigned long ulng; - double dbl; - - rng_get_bytes (rng, &ulng, sizeof ulng); - dbl = ulng / (ULONG_MAX + 1.0); - if (dbl >= 0 && dbl < 1) - return dbl; - } -} - -/* Returns a random number from the distribution with mean 0 and - standard deviation 1. (Multiply the result by the desired - standard deviation, then add the desired mean.) */ -double -rng_get_double_normal (struct rng *rng) -{ - /* Knuth, _The Art of Computer Programming_, Vol. 2, 3.4.1C, - Algorithm P. */ - double this_normal; - - if (rng->next_normal != NOT_DOUBLE) - { - this_normal = rng->next_normal; - rng->next_normal = NOT_DOUBLE; - } - else - { - double v1, v2, s; - - do - { - double u1 = rng_get_double (rng); - double u2 = rng_get_double (rng); - v1 = 2.0 * u1 - 1.0; - v2 = 2.0 * u2 - 1.0; - s = v1 * v1 + v2 * v2; - } - while (s >= 1); - - this_normal = v1 * sqrt (-2. * log (s) / s); - rng->next_normal = v2 * sqrt (-2. * log (s) / s); - } - - return this_normal; -} - -/* Gets an initialized RNG for use in PSPP transformations and - procedures. */ -struct rng * -pspp_rng (void) -{ - static struct rng *rng; - - if (rng == NULL) - rng = rng_create (); - return rng; -} diff --git a/src/random.h b/src/random.h deleted file mode 100644 index 7ab72377..00000000 --- a/src/random.h +++ /dev/null @@ -1,40 +0,0 @@ -/* PSPP - computes sample statistics. - Copyright (C) 1997-9, 2000 Free Software Foundation, Inc. - Written by Ben Pfaff . - - 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, write to the Free Software - Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA - 02111-1307, USA. */ - -#if !random_h -#define random_h 1 - -#include - -struct rng *rng_create (void); -void rng_destroy (struct rng *); -void rng_seed (struct rng *, const void *, size_t); -void rng_get_bytes (struct rng *, void *, size_t); -int rng_get_int (struct rng *); -unsigned rng_get_unsigned (struct rng *); -double rng_get_double (struct rng *); -double rng_get_double_normal (struct rng *); - -struct rng *pspp_rng (void); - -/* Return a `random' seed by using the real time clock */ -unsigned long random_seed(void); - - -#endif /* random.h */ diff --git a/src/sample.c b/src/sample.c index 8db7c837..96c5f4d9 100644 --- a/src/sample.c +++ b/src/sample.c @@ -18,6 +18,7 @@ 02111-1307, USA. */ #include +#include #include #include #include @@ -25,7 +26,7 @@ #include "command.h" #include "error.h" #include "lexer.h" -#include "random.h" +#include "settings.h" #include "str.h" #include "var.h" @@ -63,6 +64,9 @@ cmd_sample (void) return CMD_FAILURE; if (!lex_integer_p ()) { + unsigned long min = gsl_rng_min (get_rng ()); + unsigned long max = gsl_rng_max (get_rng ()); + type = TYPE_FRACTION; if (tokval <= 0 || tokval >= 1) { @@ -71,7 +75,7 @@ cmd_sample (void) return CMD_FAILURE; } - frac = tokval * UINT_MAX; + frac = tokval * (max - min) + min; a = b = 0; } else @@ -119,7 +123,7 @@ sample_trns_proc (struct trns_header * trns, struct ccase *c UNUSED, if (t->type == TYPE_FRACTION) { - if (rng_get_unsigned (pspp_rng ()) <= t->frac) + if (gsl_rng_get (get_rng ()) <= t->frac) return -1; else return -2; @@ -128,7 +132,7 @@ sample_trns_proc (struct trns_header * trns, struct ccase *c UNUSED, if (t->m >= t->n) return -2; - U = rng_get_double (pspp_rng ()); + U = gsl_rng_uniform (get_rng ()); if ((t->N - t->t) * U >= t->n - t->m) { t->t++; diff --git a/src/set.q b/src/set.q index 4c2d66a3..c291c1a8 100644 --- a/src/set.q +++ b/src/set.q @@ -64,6 +64,7 @@ #include #include #include +#include #include "alloc.h" #include "command.h" #include "lexer.h" @@ -74,7 +75,6 @@ #include "var.h" #include "format.h" #include "copyleft.h" -#include "random.h" #include "signal.h" @@ -106,8 +106,7 @@ static int set_listing=1; static char *set_pager=0; #endif /* !USE_INTERNAL_PAGER */ -static unsigned long set_seed; -static int seed_flag=0; +static gsl_rng *rng; static int long_view=0; int set_testing_mode=0; @@ -122,6 +121,9 @@ static void set_routing (int q, int *setting); static int set_ccx (const char *cc_string, struct set_cust_currency * cc, int cc_name); +static void set_rng (unsigned long); +static unsigned long random_seed (void); + /* (specification) "SET" (stc_): automenu=automenu:on/off; @@ -301,7 +303,6 @@ aux_stc_custom_results(struct cmd_set *cmd UNUSED) static int aux_stc_custom_seed(struct cmd_set *cmd UNUSED) { - msg(MM, "%ld",set_seed); return 0; } @@ -704,15 +705,14 @@ stc_custom_seed (struct cmd_set *cmd UNUSED) { lex_match ('='); if (lex_match_id ("RANDOM")) - set_seed = random_seed(); + set_rng (random_seed ()); else { if (!lex_force_num ()) return 0; - set_seed = tokval; + set_rng (tokval); lex_get (); } - seed_flag = 1; return 1; } @@ -1313,26 +1313,28 @@ get_pager(void) return set_pager; } -/* Return 1 if the seed has been set since the last time this function - was called. - Fill the value pointed to by seed with the seed . -*/ -int -seed_is_set(unsigned long *seed) +gsl_rng * +get_rng (void) { - int result = 0; - - *seed = set_seed ; - - if ( seed_flag ) - result = 1; - - seed_flag = 0; + if (rng == NULL) + set_rng (random_seed ()); + return rng; +} - return result; - +static void +set_rng (unsigned long seed) +{ + rng = gsl_rng_alloc (gsl_rng_mt19937); + if (rng == NULL) + out_of_memory (); + gsl_rng_set (rng, seed); } +static unsigned long +random_seed (void) +{ + return time (0); +} static int global_algorithm = ENHANCED; static int cmd_algorithm = ENHANCED; diff --git a/src/settings.h b/src/settings.h index 14f902ec..4c196b7e 100644 --- a/src/settings.h +++ b/src/settings.h @@ -227,11 +227,8 @@ const char *get_pager(void); #endif /* !USE_INTERNAL_PAGER */ -/* Return 1 if the seed has been set since the last time this function - was called. - Fill the value pointed to by seed with the seed . -*/ -int seed_is_set(unsigned long *seed); +#include +gsl_rng *get_rng (void); enum {ENHANCED,COMPATIBLE}; diff --git a/src/vfm.c b/src/vfm.c index 89624165..797167e6 100644 --- a/src/vfm.c +++ b/src/vfm.c @@ -35,7 +35,6 @@ #include "error.h" #include "expr.h" #include "misc.h" -#include "random.h" #include "settings.h" #include "som.h" #include "str.h" diff --git a/tests/ChangeLog b/tests/ChangeLog index ded692b8..9bc159a3 100644 --- a/tests/ChangeLog +++ b/tests/ChangeLog @@ -1,3 +1,8 @@ +Mon Nov 15 23:52:55 2004 Ben Pfaff + + * bugs/random.sh: Update expected random values to reflect the GSL + random number generator. + Sat Nov 6 14:49:27 WST 2004 John Darrington * command/oneway-with-splits.sh Added. diff --git a/tests/bugs/random.sh b/tests/bugs/random.sh index 1e0e55e2..0f4acd77 100755 --- a/tests/bugs/random.sh +++ b/tests/bugs/random.sh @@ -72,26 +72,27 @@ activity="compare output" diff -b -B -w $TEMPDIR/pspp.list - << EOF R1 -------- - 2.36 - 3.13 - 1.76 - .15 - 5.88 - 8.74 - 2.19 - 6.53 - 5.69 - 6.77 - 7.20 - 4.01 - .03 - 4.67 - 5.10 - .44 - 8.27 - 6.81 - 9.55 - 8.74 + 7.71 + 2.99 + .21 + 4.95 + 6.34 + 4.43 + 7.49 + 8.32 + 4.99 + 5.83 + 2.25 + .25 + 1.98 + 7.09 + 7.61 + 2.66 + 1.69 + 2.64 + .88 + 1.50 + EOF if [ $? -ne 0 ] ; then fail ; fi