complete suite of random number generators on expressions.
+Mon Nov 15 23:47:40 2004 Ben Pfaff <blp@gnu.org>
+
+ 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 <blp@gnu.org>
* expr-evl.c: (expr_evaluate) Fix XDATE.JDAY formula. Thanks to
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 \
#include <config.h>
#include "algorithm.h"
+#include <gsl/gsl_rng.h>
#include <limits.h>
#include <stdlib.h>
#include <string.h>
#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. */
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
casefile_destroy (casefiles);
}
\f
+#include <gsl/gsl_rng.h>
#include <stdarg.h>
#include "command.h"
-#include "random.h"
#include "lexer.h"
static void 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);
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);
casereader_destroy (r1);
}
casefile_destroy (cf);
- rng_destroy (rng);
+ gsl_rng_free (rng);
}
static void
#include "expr.h"
#include "exprP.h"
#include "error.h"
+#include <gsl/gsl_randist.h>
#include <math.h>
#include <errno.h>
#include <stdio.h>
#include "misc.h"
#include "moments.h"
#include "pool.h"
-#include "random.h"
+#include "settings.h"
#include "str.h"
#include "var.h"
#include "vfm.h"
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);
+++ /dev/null
-/* PSPP - computes sample statistics.
- Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
- Written by Ben Pfaff <blp@gnu.org>.
-
- 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 <config.h>
-#include "random.h"
-#include "error.h"
-#include <inttypes.h>
-#include <limits.h>
-#include <math.h>
-#include <stdlib.h>
-#include <time.h>
-#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;
-}
+++ /dev/null
-/* PSPP - computes sample statistics.
- Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
- Written by Ben Pfaff <blp@gnu.org>.
-
- 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 <stddef.h>
-
-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 */
02111-1307, USA. */
#include <config.h>
+#include <gsl/gsl_rng.h>
#include <limits.h>
#include <stdio.h>
#include <math.h>
#include "command.h"
#include "error.h"
#include "lexer.h"
-#include "random.h"
+#include "settings.h"
#include "str.h"
#include "var.h"
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)
{
return CMD_FAILURE;
}
- frac = tokval * UINT_MAX;
+ frac = tokval * (max - min) + min;
a = b = 0;
}
else
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;
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++;
#include <stdio.h>
#include <errno.h>
#include <stdlib.h>
+#include <time.h>
#include "alloc.h"
#include "command.h"
#include "lexer.h"
#include "var.h"
#include "format.h"
#include "copyleft.h"
-#include "random.h"
#include "signal.h"
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;
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;
static int
aux_stc_custom_seed(struct cmd_set *cmd UNUSED)
{
- msg(MM, "%ld",set_seed);
return 0;
}
{
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;
}
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;
#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/gsl_rng.h>
+gsl_rng *get_rng (void);
enum {ENHANCED,COMPATIBLE};
#include "error.h"
#include "expr.h"
#include "misc.h"
-#include "random.h"
#include "settings.h"
#include "som.h"
#include "str.h"
+Mon Nov 15 23:52:55 2004 Ben Pfaff <blp@gnu.org>
+
+ * bugs/random.sh: Update expected random values to reflect the GSL
+ random number generator.
+
Sat Nov 6 14:49:27 WST 2004 John Darrington <john@darrington.wattle.id.au>
* command/oneway-with-splits.sh Added.
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
projects / pspp-builds.git / commitdiff