-/* Copyright (C) 1999,2000,2001 Free Software Foundation, Inc.
- * 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, or (at your option)
- * any later version.
+/* Copyright (C) 1999,2000,2001, 2003 Free Software Foundation, Inc.
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
*
- * This program is distributed in the hope that it will be useful,
+ * This library 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.
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * You should have received a copy of the GNU General Public License
- * along with this software; see the file COPYING. If not, write to
- * the Free Software Foundation, Inc., 59 Temple Place, Suite 330,
- * Boston, MA 02111-1307 USA
- *
- * As a special exception, the Free Software Foundation gives permission
- * for additional uses of the text contained in its release of GUILE.
- *
- * The exception is that, if you link the GUILE library with other files
- * to produce an executable, this does not by itself cause the
- * resulting executable to be covered by the GNU General Public License.
- * Your use of that executable is in no way restricted on account of
- * linking the GUILE library code into it.
- *
- * This exception does not however invalidate any other reasons why
- * the executable file might be covered by the GNU General Public License.
- *
- * This exception applies only to the code released by the
- * Free Software Foundation under the name GUILE. If you copy
- * code from other Free Software Foundation releases into a copy of
- * GUILE, as the General Public License permits, the exception does
- * not apply to the code that you add in this way. To avoid misleading
- * anyone as to the status of such modified files, you must delete
- * this exception notice from them.
- *
- * If you write modifications of your own for GUILE, it is your choice
- * whether to permit this exception to apply to your modifications.
- * If you do not wish that, delete this exception notice. */
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
/* Author: Mikael Djurfeldt <djurfeldt@nada.kth.se> */
+#if HAVE_CONFIG_H
+# include <config.h>
+#endif
+
#include "libguile/_scm.h"
+#include <gmp.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#define M_PI 3.14159265359
#endif
-#if SIZEOF_LONG > 4
-#if SIZEOF_INT > 4
-#define LONG32 unsigned short
-#else
-#define LONG32 unsigned int
-#endif
-#define LONG64 unsigned long
-#else
-#define LONG32 unsigned long
-#ifdef __MINGW32__
-#define LONG64 unsigned __int64
-#else
-#define LONG64 unsigned long long
-#endif
-#endif
-
-#if SIZEOF_LONG > 4 || defined (HAVE_LONG_LONGS)
+#ifdef SCM_HAVE_T_INT64
unsigned long
scm_i_uniform32 (scm_t_i_rstate *state)
{
- LONG64 x = (LONG64) A * state->w + state->c;
- LONG32 w = x & 0xffffffffUL;
+ scm_t_int64 x = (scm_t_int64) A * state->w + state->c;
+ scm_t_int32 w = x & 0xffffffffUL;
state->w = w;
state->c = x >> 32L;
return w;
unsigned long
scm_i_uniform32 (scm_t_i_rstate *state)
{
- LONG32 x1 = L (A) * L (state->w);
- LONG32 x2 = L (A) * H (state->w);
- LONG32 x3 = H (A) * L (state->w);
- LONG32 w = L (x1) + L (state->c);
- LONG32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w);
- LONG32 x4 = H (A) * H (state->w);
+ scm_t_int32 x1 = L (A) * L (state->w);
+ scm_t_int32 x2 = L (A) * H (state->w);
+ scm_t_int32 x3 = H (A) * L (state->w);
+ scm_t_int32 w = L (x1) + L (state->c);
+ scm_t_int32 m = H (x1) + L (x2) + L (x3) + H (state->c) + H (w);
+ scm_t_int32 x4 = H (A) * H (state->w);
state->w = w = (L (m) << 16) + L (w);
state->c = H (x2) + H (x3) + x4 + H (m);
return w;
void
scm_i_init_rstate (scm_t_i_rstate *state, char *seed, int n)
{
- LONG32 w = 0L;
- LONG32 c = 0L;
+ scm_t_int32 w = 0L;
+ scm_t_int32 c = 0L;
int i, m;
for (i = 0; i < n; ++i)
{
scm_c_default_rstate ()
#define FUNC_NAME "scm_c_default_rstate"
{
- SCM state = SCM_CDR (scm_var_random_state);
+ SCM state = SCM_VARIABLE_REF (scm_var_random_state);
if (!SCM_RSTATEP (state))
SCM_MISC_ERROR ("*random-state* contains bogus random state", SCM_EOL);
return SCM_RSTATE (state);
return r;
}
+/*
+ SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
+
+ Takes a random state (source of random bits) and a bignum m.
+ Returns a bignum b, 0 <= b < m.
+
+ It does this by allocating a bignum b with as many base 65536 digits
+ as m, filling b with random bits (in 32 bit chunks) up to the most
+ significant 1 in m, and, finally checking if the resultant b is too
+ large (>= m). If too large, we simply repeat the process again. (It
+ is important to throw away all generated random bits if b >= m,
+ otherwise we'll end up with a distorted distribution.)
+
+*/
+
SCM
scm_c_random_bignum (scm_t_rstate *state, SCM m)
{
- SCM b;
- int i, nd;
- LONG32 *bits, mask, w;
- nd = SCM_NUMDIGS (m);
- /* calculate mask for most significant digit */
-#if SIZEOF_INT == 4
- /* 16 bit digits */
- if (nd & 1)
- {
- /* fix most significant 16 bits */
- unsigned short s = SCM_BDIGITS (m)[nd - 1];
- mask = s < 0x100 ? scm_masktab[s] : scm_masktab[s >> 8] << 8 | 0xff;
- }
- else
-#endif
- {
- /* fix most significant 32 bits */
-#if SIZEOF_INT == 4
- w = SCM_BDIGITS (m)[nd - 1] << 16 | SCM_BDIGITS (m)[nd - 2];
-#else
- w = SCM_BDIGITS (m)[nd - 1];
-#endif
- mask = (w < 0x10000
- ? (w < 0x100
- ? scm_masktab[w]
- : scm_masktab[w >> 8] << 8 | 0xff)
- : (w < 0x1000000
- ? scm_masktab[w >> 16] << 16 | 0xffff
- : scm_masktab[w >> 24] << 24 | 0xffffff));
- }
- b = scm_i_mkbig (nd, 0);
- bits = (LONG32 *) SCM_BDIGITS (b);
+ SCM result = scm_i_mkbig ();
+ const size_t m_bits = mpz_sizeinbase (SCM_I_BIG_MPZ (m), 2);
+ /* how many bits would only partially fill the last unsigned long? */
+ const size_t end_bits = m_bits % (sizeof (unsigned long) * SCM_CHAR_BIT);
+ unsigned long *random_chunks = NULL;
+ const unsigned long num_full_chunks =
+ m_bits / (sizeof (unsigned long) * SCM_CHAR_BIT);
+ const unsigned long num_chunks = num_full_chunks + ((end_bits) ? 1 : 0);
+
+ /* we know the result will be this big */
+ mpz_realloc2 (SCM_I_BIG_MPZ (result), m_bits);
+
+ random_chunks =
+ (unsigned long *) scm_gc_calloc (num_chunks * sizeof (unsigned long),
+ "random bignum chunks");
+
do
{
- i = nd;
- /* treat most significant digit specially */
-#if SIZEOF_INT == 4
- /* 16 bit digits */
- if (i & 1)
- {
- ((SCM_BIGDIG*) bits)[i - 1] = scm_the_rng.random_bits (state) & mask;
- i /= 2;
- }
- else
-#endif
- {
- /* fix most significant 32 bits */
-#if SIZEOF_INT == 4
- w = scm_the_rng.random_bits (state) & mask;
- ((SCM_BIGDIG*) bits)[i - 2] = w & 0xffff;
- ((SCM_BIGDIG*) bits)[i - 1] = w >> 16;
- i = i / 2 - 1;
-#else
- i /= 2;
- bits[--i] = scm_the_rng.random_bits (state) & mask;
-#endif
- }
- /* now fill up the rest of the bignum */
- while (i)
- bits[--i] = scm_the_rng.random_bits (state);
- b = scm_i_normbig (b);
- if (SCM_INUMP (b))
- return b;
- } while (scm_bigcomp (b, m) <= 0);
- return b;
+ unsigned long *current_chunk = random_chunks + (num_chunks - 1);
+ unsigned long chunks_left = num_chunks;
+
+ mpz_set_ui (SCM_I_BIG_MPZ (result), 0);
+
+ if (end_bits)
+ {
+ /* generate a mask with ones in the end_bits position, i.e. if
+ end_bits is 3, then we'd have a mask of ...0000000111 */
+ const unsigned long rndbits = scm_the_rng.random_bits (state);
+ int rshift = (sizeof (unsigned long) * SCM_CHAR_BIT) - end_bits;
+ unsigned long mask = ((unsigned long) ULONG_MAX) >> rshift;
+ unsigned long highest_bits = rndbits & mask;
+ *current_chunk-- = highest_bits;
+ chunks_left--;
+ }
+
+ while (chunks_left)
+ {
+ /* now fill in the remaining unsigned long sized chunks */
+ *current_chunk-- = scm_the_rng.random_bits (state);
+ chunks_left--;
+ }
+ mpz_import (SCM_I_BIG_MPZ (result),
+ num_chunks,
+ -1,
+ sizeof (unsigned long),
+ 0,
+ 0,
+ random_chunks);
+ /* if result >= m, regenerate it (it is important to regenerate
+ all bits in order not to get a distorted distribution) */
+ } while (mpz_cmp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (m)) >= 0);
+ scm_gc_free (random_chunks,
+ num_chunks * sizeof (unsigned long),
+ "random bignum chunks");
+ return scm_i_normbig (result);
}
/*
if (SCM_REALP (n))
return scm_make_real (SCM_REAL_VALUE (n)
* scm_c_uniform01 (SCM_RSTATE (state)));
- SCM_VALIDATE_SMOB (1, n, big);
+
+ SCM_VALIDATE_BIGINT (1, n);
return scm_c_random_bignum (SCM_RSTATE (state), n);
}
#undef FUNC_NAME
}
#undef FUNC_NAME
-#ifdef HAVE_ARRAYS
+#if SCM_HAVE_ARRAYS
static void
vector_scale (SCM v, double c)
}
#undef FUNC_NAME
-#endif /* HAVE_ARRAYS */
+#endif /* SCM_HAVE_ARRAYS */
SCM_DEFINE (scm_random_exp, "random:exp", 0, 1, 0,
(SCM state),
#include "libguile/random.x"
- /* Check that the assumptions about bits per bignum digit are correct. */
-#if SIZEOF_INT == 4
- m = 16;
-#else
- m = 32;
-#endif
- if (m != SCM_BITSPERDIG)
- {
- fprintf (stderr, "Internal inconsistency: Confused about bignum digit size in random.c\n");
- exit (1);
- }
-
scm_add_feature ("random");
}
HCoop / bpt / guile.git / blobdiff