[gnumeric] compilation: split random number generation into its own file.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] compilation: split random number generation into its own file.
- Date: Mon, 2 Aug 2010 18:37:47 +0000 (UTC)
commit 56f02689d378fb4d3b70049a87d4bf3e58486e0e
Author: Morten Welinder <terra gnome org>
Date: Mon Aug 2 14:08:04 2010 -0400
compilation: split random number generation into its own file.
ChangeLog | 4 +
plugins/fn-random/functions.c | 1 +
src/Makefile.am | 2 +
src/dialogs/dialog-about.c | 2 +-
src/gnm-random.c | 1377 +++++++++++++++++++++++++++++++++++++++++
src/gnm-random.h | 44 ++
src/mathfunc.c | 1370 +----------------------------------------
src/mathfunc.h | 36 +-
src/tools/data-shuffling.c | 2 +-
src/tools/goal-seek.c | 2 +-
src/tools/random-generator.c | 2 +-
src/wbc-gtk-actions.c | 2 +-
12 files changed, 1436 insertions(+), 1408 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index dcc12fb..6076b49 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2010-08-02 Morten Welinder <terra gnome org>
+
+ * src/gnm-random.c: Split from src/mathfunc.c.
+
2010-08-01 Jean Brefort <jean brefort normalesup org>
* src/item-bar.c (ib_draw_cell): use theme foreground color for row and
diff --git a/plugins/fn-random/functions.c b/plugins/fn-random/functions.c
index 214effd..920a43d 100644
--- a/plugins/fn-random/functions.c
+++ b/plugins/fn-random/functions.c
@@ -27,6 +27,7 @@
#include <gnumeric.h>
#include <func.h>
#include <mathfunc.h>
+#include <gnm-random.h>
#include <rangefunc.h>
#include <value.h>
#include <expr.h>
diff --git a/src/Makefile.am b/src/Makefile.am
index 4957b84..47e6740 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -117,6 +117,7 @@ libspreadsheet_la_SOURCES = \
gnm-graph-window.c \
gnm-pane.c \
gnm-pane-impl.h \
+ gnm-random.c \
gnumeric-simple-canvas.c \
graph.c \
gutils.c \
@@ -228,6 +229,7 @@ libspreadsheet_include_HEADERS = \
gnm-format.h \
gnm-graph-window.h \
gnm-pane.h \
+ gnm-random.h \
gnm-sheet-slicer.h \
gnm-style-impl.h \
gnumeric.h \
diff --git a/src/dialogs/dialog-about.c b/src/dialogs/dialog-about.c
index e8e150e..9f1b95a 100644
--- a/src/dialogs/dialog-about.c
+++ b/src/dialogs/dialog-about.c
@@ -27,7 +27,7 @@
#include <string.h>
#include <gui-util.h>
-#include <mathfunc.h>
+#include <gnm-random.h>
#include <wbc-gtk.h>
#include <gnm-format.h>
#include <goffice/goffice.h>
diff --git a/src/gnm-random.c b/src/gnm-random.c
new file mode 100644
index 0000000..ee73287
--- /dev/null
+++ b/src/gnm-random.c
@@ -0,0 +1,1377 @@
+/* for random() */
+#define _SVID_SOURCE 1
+#define _BSD_SOURCE 1
+
+#include <gnumeric-config.h>
+#include "gnumeric.h"
+#include "gnm-random.h"
+#include "mathfunc.h"
+#include <glib/gstdio.h>
+
+#include <fcntl.h>
+#include <unistd.h>
+#include <string.h>
+
+/*
+ * mathfunc.c: Mathematical functions.
+ *
+ * Authors:
+ * Morten Welinder <terra gnome org>
+ * Makoto Matsumoto and Takuji Nishimura (Mersenne Twister, see note in code)
+ * James Theiler. (See note 1.)
+ * Brian Gough. (See note 1.)
+ *
+ *
+ * NOTE 1: most of the random distribution code comes from the GNU Scientific
+ * Library (GSL), notably version 1.1.1. GSL is distributed under GPL licence,
+ * see COPYING. The relevant parts are copyright (C) 1996, 1997, 1998, 1999,
+ * 2000 James Theiler and Brian Gough.
+ *
+ * Thank you!
+ */
+
+/* ------------------------------------------------------------------------- */
+/* http://www.math.keio.ac.jp/matumoto/CODES/MT2002/mt19937ar.c */
+/* Imported by hand -- MW. */
+
+/*
+ A C-program for MT19937, with initialization improved 2002/1/26.
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+
+ Before using, initialize the state by using init_genrand(seed)
+ or init_by_array(init_key, key_length).
+
+ Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ 3. The names of its contributors may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+ Any feedback is very welcome.
+ http://www.math.keio.ac.jp/matumoto/emt.html
+ email: matumoto math keio ac jp
+*/
+
+#if 0
+#include <stdio.h>
+#endif
+
+/* Period parameters */
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL /* constant vector a */
+#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+
+static unsigned long mt[N]; /* the array for the state vector */
+static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
+
+/* initializes mt[N] with a seed */
+static void init_genrand(unsigned long s)
+{
+ mt[0]= s & 0xffffffffUL;
+ for (mti=1; mti<N; mti++) {
+ mt[mti] =
+ (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
+ /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+ /* In the previous versions, MSBs of the seed affect */
+ /* only MSBs of the array mt[]. */
+ /* 2002/01/09 modified by Makoto Matsumoto */
+ mt[mti] &= 0xffffffffUL;
+ /* for >32 bit machines */
+ }
+}
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+static void mt_init_by_array(unsigned long init_key[], int key_length)
+{
+ int i, j, k;
+ init_genrand(19650218UL);
+ i=1; j=0;
+ k = (N>key_length ? N : key_length);
+ for (; k; k--) {
+ mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ + init_key[j] + j; /* non linear */
+ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++; j++;
+ if (i>=N) { mt[0] = mt[N-1]; i=1; }
+ if (j>=key_length) j=0;
+ }
+ for (k=N-1; k; k--) {
+ mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
+ - i; /* non linear */
+ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++;
+ if (i>=N) { mt[0] = mt[N-1]; i=1; }
+ }
+
+ mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
+}
+
+/* generates a random number on [0,0xffffffff]-interval */
+static unsigned long genrand_int32(void)
+{
+ unsigned long y;
+ static unsigned long mag01[2]={0x0UL, MATRIX_A};
+ /* mag01[x] = x * MATRIX_A for x=0,1 */
+
+ if (mti >= N) { /* generate N words at one time */
+ int kk;
+
+ if (mti == N+1) /* if init_genrand() has not been called, */
+ init_genrand(5489UL); /* a default initial seed is used */
+
+ for (kk=0;kk<N-M;kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+ }
+ for (;kk<N-1;kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+ }
+ y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+ mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+ mti = 0;
+ }
+
+ y = mt[mti++];
+
+ /* Tempering */
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9d2c5680UL;
+ y ^= (y << 15) & 0xefc60000UL;
+ y ^= (y >> 18);
+
+ return y;
+}
+
+#if 0
+/* generates a random number on [0,0x7fffffff]-interval */
+long genrand_int31(void)
+{
+ return (long)(genrand_int32()>>1);
+}
+
+/* generates a random number on [0,1]-real-interval */
+double genrand_real1(void)
+{
+ return genrand_int32()*(1.0/4294967295.0);
+ /* divided by 2^32-1 */
+}
+
+/* generates a random number on [0,1)-real-interval */
+double genrand_real2(void)
+{
+ return genrand_int32()*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/* generates a random number on (0,1)-real-interval */
+double genrand_real3(void)
+{
+ return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+#endif
+
+/* generates a random number on [0,1) with 53-bit resolution*/
+static double genrand_res53(void)
+{
+ unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
+ return(a*67108864.0+b)*(1.0/9007199254740992.0);
+}
+/* These real versions are due to Isaku Wada, 2002/01/09 added */
+
+#if 0
+int main(void)
+{
+ int i;
+ unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
+ init_by_array(init, length);
+ printf("1000 outputs of genrand_int32()\n");
+ for (i=0; i<1000; i++) {
+ printf("%10lu ", genrand_int32());
+ if (i%5==4) printf("\n");
+ }
+ printf("\n1000 outputs of genrand_real2()\n");
+ for (i=0; i<1000; i++) {
+ printf("%10.8f ", genrand_real2());
+ if (i%5==4) printf("\n");
+ }
+ return 0;
+}
+#endif
+
+
+#undef N
+#undef M
+#undef MATRIX_A
+#undef UPPER_MASK
+#undef LOWER_MASK
+
+/* ------------------------------------------------------------------------ */
+
+/* FIXME: we need something that catches partials and EAGAIN. */
+#define fullread read
+
+#define RANDOM_DEVICE "/dev/urandom"
+
+/*
+ * Conservative random number generator. The result is (supposedly) uniform
+ * and between 0 and 1. (0 possible, 1 not.) The result should have about
+ * 64 bits randomness.
+ */
+gnm_float
+random_01 (void)
+{
+ static int device_fd = -2;
+ static int seeded = -2;
+
+ while (seeded) {
+ if (seeded == -2) {
+ char const *seed = g_getenv ("GNUMERIC_PRNG_SEED");
+ if (seed) {
+ int len = strlen (seed);
+ int i;
+ unsigned long *longs = g_new (unsigned long, len + 1);
+
+ /* We drop only one character into each long. */
+ for (i = 0; i < len; i++)
+ longs[i] = (unsigned char)seed[i];
+ mt_init_by_array (longs, len);
+ g_free (longs);
+ seeded = 1;
+
+ g_warning ("Using pseudo-random numbers.");
+ } else {
+ seeded = 0;
+ break;
+ }
+ }
+
+ /*
+ * Only 52-bit precision. But hey, if you are using pseudo
+ * random numbers that ought to be good enough to you.
+ */
+ return genrand_res53 ();
+ }
+
+ if (device_fd == -2) {
+ device_fd = g_open (RANDOM_DEVICE, O_RDONLY, 0);
+ /*
+ * We could check that we really have a device, but it hard
+ * to come up with a non-paranoid reason to.
+ */
+ }
+
+ if (device_fd >= 0) {
+ static ssize_t bytes_left = 0;
+ static unsigned char data[32 * sizeof (gnm_float)];
+ gnm_float res = 0;
+ size_t i;
+
+ if (bytes_left < (ssize_t)sizeof (gnm_float)) {
+ ssize_t got = fullread (device_fd, &data, sizeof (data));
+ if (got < (ssize_t)sizeof (gnm_float))
+ goto failure;
+ bytes_left += got;
+ }
+
+ bytes_left -= sizeof (gnm_float);
+ for (i = 0; i < sizeof (gnm_float); i++)
+ res = (res + data[bytes_left + i]) / 256;
+ return res;
+
+ failure:
+ /* It failed when it shouldn't. Disable. */
+ g_warning ("Reading from %s failed; reverting to pseudo-random.",
+ RANDOM_DEVICE);
+ close (device_fd);
+ device_fd = -1;
+ }
+
+ return genrand_res53 ();
+}
+
+/*
+ * Generate a N(0,1) distributed number.
+ */
+gnm_float
+random_normal (void)
+{
+ static gboolean has_saved = FALSE;
+ static gnm_float saved;
+
+ if (has_saved) {
+ has_saved = FALSE;
+ return saved;
+ } else {
+ gnm_float u, v, r2, rsq;
+ do {
+ u = 2 * random_01 () - 1;
+ v = 2 * random_01 () - 1;
+ r2 = u * u + v * v;
+ } while (r2 > 1 || r2 == 0);
+
+ rsq = gnm_sqrt (-2 * gnm_log (r2) / r2);
+
+ has_saved = TRUE;
+ saved = v * rsq;
+
+ return u * rsq;
+ }
+}
+
+gnm_float
+random_lognormal (gnm_float zeta, gnm_float sigma)
+{
+ return gnm_exp (sigma * random_normal () + zeta);
+}
+
+static gnm_float
+random_gaussian (gnm_float sigma)
+{
+ return sigma * random_normal ();
+}
+
+/*
+ * Generate a poisson distributed number.
+ */
+gnm_float
+random_poisson (gnm_float lambda)
+{
+ /*
+ * This may not be optimal code, but it sure is easy to
+ * understand compared to R's code.
+ */
+ return qpois (random_01 (), lambda, TRUE, FALSE);
+}
+
+/*
+ * Generate a binomial distributed number.
+ */
+gnm_float
+random_binomial (gnm_float p, gnm_float trials)
+{
+ return qbinom (random_01 (), trials, p, TRUE, FALSE);
+}
+
+/*
+ * Generate a negative binomial distributed number.
+ */
+gnm_float
+random_negbinom (gnm_float p, gnm_float f)
+{
+ return qnbinom (random_01 (), f, p, TRUE, FALSE);
+}
+
+/*
+ * Generate an exponential distributed number.
+ */
+gnm_float
+random_exponential (gnm_float b)
+{
+ return -b * gnm_log (random_01 ());
+}
+
+/*
+ * Generate a bernoulli distributed number.
+ */
+gnm_float
+random_bernoulli (gnm_float p)
+{
+ gnm_float r = random_01 ();
+
+ return (r <= p) ? 1.0 : 0.0;
+}
+
+/*
+ * Generate a cauchy distributed number. From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_cauchy (gnm_float a)
+{
+ gnm_float u;
+
+ do {
+ u = random_01 ();
+ } while (u == 0.5);
+
+ return a * gnm_tan (M_PIgnum * u);
+}
+
+/*
+ * Generate a Weibull distributed number. From the GNU Scientific
+ * library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_weibull (gnm_float a, gnm_float b)
+{
+ gnm_float x, z;
+
+ do {
+ x = random_01 ();
+ } while (x == 0.0);
+
+ z = gnm_pow (-gnm_log (x), 1 / b);
+
+ return a * z;
+}
+
+/*
+ * Generate a Laplace (two-sided exponential probability) distributed number.
+ * From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_laplace (gnm_float a)
+{
+ gnm_float u;
+
+ do {
+ u = 2 * random_01 () - 1.0;
+ } while (u == 0.0);
+
+ if (u < 0)
+ return a * gnm_log (-u);
+ else
+ return -a * gnm_log (u);
+}
+
+gnm_float
+random_laplace_pdf (gnm_float x, gnm_float a)
+{
+ return (1 / (2 * a)) * gnm_exp (-gnm_abs (x) / a);
+}
+
+/*
+ * Generate a Rayleigh distributed number. From the GNU Scientific library
+ * 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_rayleigh (gnm_float sigma)
+{
+ gnm_float u;
+
+ do {
+ u = random_01 ();
+ } while (u == 0.0);
+
+ return sigma * gnm_sqrt (-2.0 * gnm_log (u));
+}
+
+/*
+ * Generate a Rayleigh tail distributed number. From the GNU Scientific library
+ * 1.1.1. The Rayleigh tail distribution has the form
+ * p(x) dx = (x / sigma^2) exp((a^2 - x^2)/(2 sigma^2)) dx
+ *
+ * for x = a ... +infty
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_rayleigh_tail (gnm_float a, gnm_float sigma)
+{
+ gnm_float u;
+
+ do {
+ u = random_01 ();
+ } while (u == 0.0);
+
+ return gnm_sqrt (a * a - 2.0 * sigma * sigma * gnm_log (u));
+}
+
+/* The Gamma distribution of order a>0 is defined by:
+ *
+ * p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx
+ *
+ * for x>0. If X and Y are independent gamma-distributed random
+ * variables of order a1 and a2 with the same scale parameter b, then
+ * X+Y has gamma distribution of order a1+a2.
+ *
+ * The algorithms below are from Knuth, vol 2, 2nd ed, p. 129.
+ */
+
+static gnm_float
+gamma_frac (gnm_float a)
+{
+ /* This is exercise 16 from Knuth; see page 135, and the solution is
+ * on page 551. */
+
+ gnm_float x, q;
+ gnm_float p = M_Egnum / (a + M_Egnum);
+ do {
+ gnm_float v;
+ gnm_float u = random_01 ();
+ do {
+ v = random_01 ();
+ } while (v == 0.0);
+
+ if (u < p) {
+ x = gnm_pow (v, 1 / a);
+ q = gnm_exp (-x);
+ } else {
+ x = 1 - gnm_log (v);
+ q = gnm_pow (x, a - 1);
+ }
+ } while (random_01 () >= q);
+
+ return x;
+}
+
+static gnm_float
+gamma_large (gnm_float a)
+{
+ /*
+ * Works only if a > 1, and is most efficient if a is large
+ *
+ * This algorithm, reported in Knuth, is attributed to Ahrens. A
+ * faster one, we are told, can be found in: J. H. Ahrens and
+ * U. Dieter, Computing 12 (1974) 223-246.
+ */
+
+ gnm_float sqa, x, y, v;
+ sqa = gnm_sqrt (2 * a - 1);
+ do {
+ do {
+ y = gnm_tan (M_PIgnum * random_01 ());
+ x = sqa * y + a - 1;
+ } while (x <= 0);
+ v = random_01 ();
+ } while (v > (1 + y * y) * gnm_exp ((a - 1) * gnm_log (x / (a - 1)) -
+ sqa * y));
+
+ return x;
+}
+
+static gnm_float
+ran_gamma_int (gnm_float a)
+{
+ if (a < 12) {
+ gnm_float prod;
+
+ do {
+ unsigned int i, ua;
+ prod = 1;
+ ua = (unsigned int)a;
+
+ for (i = 0; i < ua; i++)
+ prod *= random_01 ();
+
+ /*
+ * This handles the 0-probability event of getting
+ * an actual zero as well as the possibility of
+ * underflow.
+ */
+ } while (prod == 0);
+
+ return -gnm_log (prod);
+ } else
+ return gamma_large (a);
+}
+
+/*
+ * Generate a Gamma distributed number. From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_gamma (gnm_float a, gnm_float b)
+{
+ gnm_float na;
+
+ if (gnm_isnan (a) || gnm_isnan (b) || a <= 0 || b <= 0)
+ return gnm_nan;
+
+ na = gnm_floor (a);
+
+ if (a == na)
+ return b * ran_gamma_int (na);
+ else if (na == 0)
+ return b * gamma_frac (a);
+ else
+ return b * (ran_gamma_int (na) + gamma_frac (a - na));
+}
+
+/*
+ * Generate a Pareto distributed number. From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_pareto (gnm_float a, gnm_float b)
+{
+ gnm_float x;
+
+ do {
+ x = random_01 ();
+ } while (x == 0.0);
+
+ return b * gnm_pow (x, -1 / a);
+}
+
+/*
+ * Generate a F-distributed number. From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_fdist (gnm_float nu1, gnm_float nu2)
+{
+ gnm_float Y1 = random_gamma (nu1 / 2, 2.0);
+ gnm_float Y2 = random_gamma (nu2 / 2, 2.0);
+
+ return (Y1 * nu2) / (Y2 * nu1);
+}
+
+/*
+ * Generate a Beta-distributed number. From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_beta (gnm_float a, gnm_float b)
+{
+ gnm_float x1 = random_gamma (a, 1.0);
+ gnm_float x2 = random_gamma (b, 1.0);
+
+ return x1 / (x1 + x2);
+}
+
+/*
+ * Generate a Chi-Square-distributed number. From the GNU Scientific library
+ * 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_chisq (gnm_float nu)
+{
+ return 2 * random_gamma (nu / 2, 1.0);
+}
+
+/*
+ * Generate a logistic-distributed number. From the GNU Scientific library
+ * 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_logistic (gnm_float a)
+{
+ gnm_float x;
+
+ do {
+ x = random_01 ();
+ } while (x == 0 || x == 1);
+
+ return a * gnm_log (x / (1 - x));
+}
+
+/*
+ * Generate a geometric-distributed number. From the GNU Scientific library
+ * 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_geometric (gnm_float p)
+{
+ gnm_float u;
+
+ if (p == 1)
+ return 1;
+ do {
+ u = random_01 ();
+ } while (u == 0);
+
+ return gnm_floor (gnm_log (u) / gnm_log1p (-p) + 1);
+}
+
+gnm_float
+random_hypergeometric (gnm_float n1, gnm_float n2, gnm_float t)
+{
+ return qhyper (random_01 (), n1, n2, t, TRUE, FALSE);
+}
+
+
+/*
+ * Generate a logarithmic-distributed number. From the GNU Scientific library
+ * 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_logarithmic (gnm_float p)
+{
+ gnm_float c, v;
+
+ c = gnm_log1p (-p);
+ do {
+ v = random_01 ();
+ } while (v == 0);
+
+ if (v >= p)
+ return 1;
+ else {
+ gnm_float u, q;
+
+ do {
+ u = random_01 ();
+ } while (u == 0);
+ q = gnm_expm1 (c * u);
+
+ if (v <= q * q)
+ return gnm_floor (1 + gnm_log (v) / gnm_log (q));
+ else if (v <= q)
+ return 2;
+ else
+ return 1;
+ }
+}
+
+/*
+ * Generate a T-distributed number. From the GNU Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_tdist (gnm_float nu)
+{
+ if (nu <= 2) {
+ gnm_float Y1 = random_normal ();
+ gnm_float Y2 = random_chisq (nu);
+
+ gnm_float t = Y1 / gnm_sqrt (Y2 / nu);
+
+ return t;
+ } else {
+ gnm_float Y1, Y2, Z, t;
+ do {
+ Y1 = random_normal ();
+ Y2 = random_exponential (1 / (nu / 2 - 1));
+
+ Z = Y1 * Y1 / (nu - 2);
+ } while (1 - Z < 0 || gnm_exp (-Y2 - Z) > (1 - Z));
+
+ /* Note that there is a typo in Knuth's formula, the line below
+ * is taken from the original paper of Marsaglia, Mathematics
+ * of Computation, 34 (1980), p 234-256. */
+
+ t = Y1 / gnm_sqrt ((1 - 2 / nu) * (1 - Z));
+ return t;
+ }
+}
+
+/*
+ * Generate a Type I Gumbel-distributed random number. From the GNU
+ * Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_gumbel1 (gnm_float a, gnm_float b)
+{
+ gnm_float x;
+
+ do {
+ x = random_01 ();
+ } while (x == 0.0);
+
+ return (gnm_log (b) - gnm_log (-gnm_log (x))) / a;
+}
+
+/*
+ * Generate a Type II Gumbel-distributed random number. From the GNU
+ * Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ */
+gnm_float
+random_gumbel2 (gnm_float a, gnm_float b)
+{
+ gnm_float x;
+
+ do {
+ x = random_01 ();
+ } while (x == 0.0);
+
+ return gnm_pow (-b / gnm_log (x), 1 / a);
+}
+
+/*
+ * Generate a stable Levy-distributed random number. From the GNU
+ * Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
+ *
+ * The stable Levy probability distributions have the form
+ *
+ * p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
+ *
+ * with 0 < alpha <= 2.
+ *
+ * For alpha = 1, we get the Cauchy distribution
+ * For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.
+ *
+ * Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
+ * Simulation". The original reference given there is,
+ *
+ * J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
+ * simulating stable random variates". Journal of the American
+ * Statistical Association, JASA 71 340-344 (1976).
+ */
+gnm_float
+random_levy (gnm_float c, gnm_float alpha)
+{
+ gnm_float u, v, t, s;
+
+ do {
+ u = random_01 ();
+ } while (u == 0.0);
+
+ u = M_PIgnum * (u - 0.5);
+
+ if (alpha == 1) { /* cauchy case */
+ t = gnm_tan (u);
+ return c * t;
+ }
+
+ do {
+ v = random_exponential (1.0);
+ } while (v == 0);
+
+ if (alpha == 2) { /* gaussian case */
+ t = 2 * gnm_sin (u) * gnm_sqrt (v);
+ return c * t;
+ }
+
+ /* general case */
+
+ t = gnm_sin (alpha * u) / gnm_pow (gnm_cos (u), 1 / alpha);
+ s = gnm_pow (gnm_cos ((1 - alpha) * u) / v, (1 - alpha) / alpha);
+
+ return c * t * s;
+}
+
+/*
+ * The following routine for the skew-symmetric case was provided by
+ * Keith Briggs.
+ *
+ * The stable Levy probability distributions have the form
+ *
+ * 2*pi* p(x) dx
+ *
+ * = int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for
+ * alpha != 1
+ * = int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|))) for
+ alpha == 1
+ *
+ * with 0<alpha<=2, -1<=beta<=1, sigma>0.
+ *
+ * For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
+ *
+ * For alpha = 1, beta=0, we get the Lorentz distribution
+ * For alpha = 2, beta=0, we get the Gaussian distribution
+ *
+ * See A. Weron and R. Weron: Computer simulation of Lévy alpha-stable
+ * variables and processes, preprint Technical University of Wroclaw.
+ * http://www.im.pwr.wroc.pl/~hugo/Publications.html
+ */
+gnm_float
+random_levy_skew (gnm_float c, gnm_float alpha, gnm_float beta)
+{
+ gnm_float V, W, X;
+
+ if (beta == 0) /* symmetric case */
+ return random_levy (c, alpha);
+
+ do {
+ V = random_01 ();
+ } while (V == 0.0);
+
+ V = M_PIgnum * (V - 0.5);
+
+ do {
+ W = random_exponential (1.0);
+ } while (W == 0);
+
+ if (alpha == 1) {
+ X = ((M_PI_2gnum + beta * V) * gnm_tan (V) -
+ beta * gnm_log (M_PI_2gnum * W * gnm_cos (V) /
+ (M_PI_2gnum + beta * V))) / M_PI_2gnum;
+ return c * (X + beta * gnm_log (c) / M_PI_2gnum);
+ } else {
+ gnm_float t = beta * gnm_tan (M_PI_2gnum * alpha);
+ gnm_float B = gnm_atan (t) / alpha;
+ gnm_float S = pow1p (t * t, 1 / (2 * alpha));
+
+ X = S * gnm_sin (alpha * (V + B)) / gnm_pow (gnm_cos (V),
+ 1 / alpha)
+ * gnm_pow (gnm_cos (V - alpha * (V + B)) / W,
+ (1 - alpha) / alpha);
+ return c * X;
+ }
+}
+
+gnm_float
+random_exppow_pdf (gnm_float x, gnm_float a, gnm_float b)
+{
+ gnm_float lngamma = lgamma1p (1 / b);
+
+ return (1 / (2 * a)) * gnm_exp (-gnm_pow (gnm_abs (x / a), b) - lngamma);
+}
+
+/*
+ * The exponential power probability distribution is
+ *
+ * p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx
+ *
+ * for -infty < x < infty. For b = 1 it reduces to the Laplace
+ * distribution.
+ *
+ * The exponential power distribution is related to the gamma
+ * distribution by E = a * pow(G(1/b),1/b), where E is an exponential
+ * power variate and G is a gamma variate.
+ *
+ * We use this relation for b < 1. For b >=1 we use rejection methods
+ * based on the laplace and gaussian distributions which should be
+ * faster.
+ *
+ * See P. R. Tadikamalla, "Random Sampling from the Exponential Power
+ * Distribution", Journal of the American Statistical Association,
+ * September 1980, Volume 75, Number 371, pages 683-686.
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
+ */
+
+gnm_float
+random_exppow (gnm_float a, gnm_float b)
+{
+ /* See http://www.mcgill.ca/files/economics/propertiesandestimation.pdf */
+ if (!(a > 0) || gnm_isnan (b))
+ return gnm_nan;
+
+ if (b < 1) {
+ gnm_float u = random_01 ();
+ gnm_float v = random_gamma (1 / b, 1.0);
+ gnm_float z = a * gnm_pow (v, 1 / b) ;
+
+ if (u > 0.5)
+ return z;
+ else
+ return -z;
+ } else if (b == 1)
+ return random_laplace (a); /* Laplace distribution */
+ else if (b < 2) {
+ /* Use laplace distribution for rejection method */
+ gnm_float x, y, h, ratio, u;
+
+ /* Scale factor chosen by upper bound on ratio at b = 2 */
+ gnm_float s = 1.4489;
+ do {
+ x = random_laplace (a);
+ y = random_laplace_pdf (x, a);
+ h = random_exppow_pdf (x, a, b);
+ ratio = h / (s * y);
+ u = random_01 ();
+ } while (u > ratio);
+
+ return x ;
+ } else if (b == 2) /* Gaussian distribution */
+ return random_gaussian (a / gnm_sqrt (2.0));
+ else {
+ /* Use gaussian for rejection method */
+ gnm_float x, y, h, ratio, u;
+ const gnm_float sigma = a / gnm_sqrt (2.0);
+
+ /* Scale factor chosen by upper bound on ratio at b = infinity.
+ * This could be improved by using a rational function
+ * approximation to the bounding curve. */
+
+ gnm_float s = 2.4091 ; /* this is sqrt(pi) e / 2 */
+
+ do {
+ x = random_gaussian (sigma) ;
+ y = dnorm (x, 0.0, gnm_abs (sigma), FALSE) ;
+ h = random_exppow_pdf (x, a, b) ;
+ ratio = h / (s * y) ;
+ u = random_01 ();
+ } while (u > ratio);
+
+ return x;
+ }
+}
+
+/*
+ * Generate a Gaussian tail-distributed random number. From the GNU
+ * Scientific library 1.1.1.
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
+ */
+gnm_float
+random_gaussian_tail (gnm_float a, gnm_float sigma)
+{
+ /*
+ * Returns a gaussian random variable larger than a
+ * This implementation does one-sided upper-tailed deviates.
+ */
+
+ gnm_float s = a / sigma;
+
+ if (s < 1) {
+ /* For small s, use a direct rejection method. The limit s < 1
+ * can be adjusted to optimise the overall efficiency */
+
+ gnm_float x;
+
+ do {
+ x = random_gaussian (1.0);
+ } while (x < s);
+ return x * sigma;
+ } else {
+ /* Use the "supertail" deviates from the last two steps
+ * of Marsaglia's rectangle-wedge-tail method, as described
+ * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11,
+ * p139, and the solution, p586.)
+ */
+
+ gnm_float u, v, x;
+
+ do {
+ u = random_01 ();
+ do {
+ v = random_01 ();
+ } while (v == 0.0);
+ x = gnm_sqrt (s * s - 2 * gnm_log (v));
+ } while (x * u > s);
+ return x * sigma;
+ }
+}
+
+/*
+ * Generate a Landau-distributed random number. From the GNU Scientific
+ * library 1.1.1.
+ *
+ * Copyright (C) 2001 David Morrison
+ *
+ * Adapted from the CERN library routines DENLAN, RANLAN, and DISLAN
+ * as described in http://consult.cern.ch/shortwrups/g110/top.html.
+ * Original author: K.S. K\"olbig.
+ *
+ * The distribution is given by the complex path integral,
+ *
+ * p(x) = (1/(2 pi i)) \int_{c-i\inf}^{c+i\inf} ds exp(s log(s) + x s)
+ *
+ * which can be converted into a real integral over [0,+\inf]
+ *
+ * p(x) = (1/pi) \int_0^\inf dt \exp(-t log(t) - x t) sin(pi t)
+ */
+
+gnm_float
+random_landau (void)
+{
+ static gnm_float F[983] = {
+ 0.0000000, /*
+ * Add empty element [0] to account for difference
+ * between C and Fortran convention for lower bound.
+ */
+ 00.000000, 00.000000, 00.000000, 00.000000, 00.000000,
+ -2.244733, -2.204365, -2.168163, -2.135219, -2.104898,
+ -2.076740, -2.050397, -2.025605, -2.002150, -1.979866,
+ -1.958612, -1.938275, -1.918760, -1.899984, -1.881879,
+ -1.864385, -1.847451, -1.831030, -1.815083, -1.799574,
+ -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
+ -1.714187, -1.701029, -1.688130, -1.675477, -1.663057,
+ -1.650858, -1.638868, -1.627078, -1.615477, -1.604058,
+ -1.592811, -1.581729, -1.570806, -1.560034, -1.549407,
+ -1.538919, -1.528565, -1.518339, -1.508237, -1.498254,
+ -1.488386, -1.478628, -1.468976, -1.459428, -1.449979,
+ -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
+ -1.395194, -1.386356, -1.377594, -1.368906, -1.360291,
+ -1.351746, -1.343269, -1.334859, -1.326512, -1.318229,
+ -1.310006, -1.301843, -1.293737, -1.285688, -1.277693,
+ -1.269752, -1.261863, -1.254024, -1.246235, -1.238494,
+ -1.230800, -1.223153, -1.215550, -1.207990, -1.200474,
+ -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
+ -1.156220, -1.148977, -1.141770, -1.134598, -1.127459,
+ -1.120354, -1.113282, -1.106242, -1.099233, -1.092255,
+ -1.085306, -1.078388, -1.071498, -1.064636, -1.057802,
+ -1.050996, -1.044215, -1.037461, -1.030733, -1.024029,
+ -1.017350, -1.010695, -1.004064, -0.997456, -0.990871,
+ -0.984308, -0.977767, -0.971247, -0.964749, -0.958271,
+ -0.951813, -0.945375, -0.938957, -0.932558, -0.926178,
+ -0.919816, -0.913472, -0.907146, -0.900838, -0.894547,
+ -0.888272, -0.882014, -0.875773, -0.869547, -0.863337,
+ -0.857142, -0.850963, -0.844798, -0.838648, -0.832512,
+ -0.826390, -0.820282, -0.814187, -0.808106, -0.802038,
+ -0.795982, -0.789940, -0.783909, -0.777891, -0.771884,
+ -0.765889, -0.759906, -0.753934, -0.747973, -0.742023,
+ -0.736084, -0.730155, -0.724237, -0.718328, -0.712429,
+ -0.706541, -0.700661, -0.694791, -0.688931, -0.683079,
+ -0.677236, -0.671402, -0.665576, -0.659759, -0.653950,
+ -0.648149, -0.642356, -0.636570, -0.630793, -0.625022,
+ -0.619259, -0.613503, -0.607754, -0.602012, -0.596276,
+ -0.590548, -0.584825, -0.579109, -0.573399, -0.567695,
+ -0.561997, -0.556305, -0.550618, -0.544937, -0.539262,
+ -0.533592, -0.527926, -0.522266, -0.516611, -0.510961,
+ -0.505315, -0.499674, -0.494037, -0.488405, -0.482777,
+ -0.477153, -0.471533, -0.465917, -0.460305, -0.454697,
+ -0.449092, -0.443491, -0.437893, -0.432299, -0.426707,
+ -0.421119, -0.415534, -0.409951, -0.404372, -0.398795,
+ -0.393221, -0.387649, -0.382080, -0.376513, -0.370949,
+ -0.365387, -0.359826, -0.354268, -0.348712, -0.343157,
+ -0.337604, -0.332053, -0.326503, -0.320955, -0.315408,
+ -0.309863, -0.304318, -0.298775, -0.293233, -0.287692,
+ -0.282152, -0.276613, -0.271074, -0.265536, -0.259999,
+ -0.254462, -0.248926, -0.243389, -0.237854, -0.232318,
+ -0.226783, -0.221247, -0.215712, -0.210176, -0.204641,
+ -0.199105, -0.193568, -0.188032, -0.182495, -0.176957,
+ -0.171419, -0.165880, -0.160341, -0.154800, -0.149259,
+ -0.143717, -0.138173, -0.132629, -0.127083, -0.121537,
+ -0.115989, -0.110439, -0.104889, -0.099336, -0.093782,
+ -0.088227, -0.082670, -0.077111, -0.071550, -0.065987,
+ -0.060423, -0.054856, -0.049288, -0.043717, -0.038144,
+ -0.032569, -0.026991, -0.021411, -0.015828, -0.010243,
+ -0.004656, 00.000934, 00.006527, 00.012123, 00.017722,
+ 00.023323, 00.028928, 00.034535, 00.040146, 00.045759,
+ 00.051376, 00.056997, 00.062620, 00.068247, 00.073877,
+ 00.079511, 00.085149, 00.090790, 00.096435, 00.102083,
+ 00.107736, 00.113392, 00.119052, 00.124716, 00.130385,
+ 00.136057, 00.141734, 00.147414, 00.153100, 00.158789,
+ 00.164483, 00.170181, 00.175884, 00.181592, 00.187304,
+ 00.193021, 00.198743, 00.204469, 00.210201, 00.215937,
+ 00.221678, 00.227425, 00.233177, 00.238933, 00.244696,
+ 00.250463, 00.256236, 00.262014, 00.267798, 00.273587,
+ 00.279382, 00.285183, 00.290989, 00.296801, 00.302619,
+ 00.308443, 00.314273, 00.320109, 00.325951, 00.331799,
+ 00.337654, 00.343515, 00.349382, 00.355255, 00.361135,
+ 00.367022, 00.372915, 00.378815, 00.384721, 00.390634,
+ 00.396554, 00.402481, 00.408415, 00.414356, 00.420304,
+ 00.426260, 00.432222, 00.438192, 00.444169, 00.450153,
+ 00.456145, 00.462144, 00.468151, 00.474166, 00.480188,
+ 00.486218, 00.492256, 00.498302, 00.504356, 00.510418,
+ 00.516488, 00.522566, 00.528653, 00.534747, 00.540850,
+ 00.546962, 00.553082, 00.559210, 00.565347, 00.571493,
+ 00.577648, 00.583811, 00.589983, 00.596164, 00.602355,
+ 00.608554, 00.614762, 00.620980, 00.627207, 00.633444,
+ 00.639689, 00.645945, 00.652210, 00.658484, 00.664768,
+ 00.671062, 00.677366, 00.683680, 00.690004, 00.696338,
+ 00.702682, 00.709036, 00.715400, 00.721775, 00.728160,
+ 00.734556, 00.740963, 00.747379, 00.753807, 00.760246,
+ 00.766695, 00.773155, 00.779627, 00.786109, 00.792603,
+ 00.799107, 00.805624, 00.812151, 00.818690, 00.825241,
+ 00.831803, 00.838377, 00.844962, 00.851560, 00.858170,
+ 00.864791, 00.871425, 00.878071, 00.884729, 00.891399,
+ 00.898082, 00.904778, 00.911486, 00.918206, 00.924940,
+ 00.931686, 00.938446, 00.945218, 00.952003, 00.958802,
+ 00.965614, 00.972439, 00.979278, 00.986130, 00.992996,
+ 00.999875, 01.006769, 01.013676, 01.020597, 01.027533,
+ 01.034482, 01.041446, 01.048424, 01.055417, 01.062424,
+ 01.069446, 01.076482, 01.083534, 01.090600, 01.097681,
+ 01.104778, 01.111889, 01.119016, 01.126159, 01.133316,
+ 01.140490, 01.147679, 01.154884, 01.162105, 01.169342,
+ 01.176595, 01.183864, 01.191149, 01.198451, 01.205770,
+ 01.213105, 01.220457, 01.227826, 01.235211, 01.242614,
+ 01.250034, 01.257471, 01.264926, 01.272398, 01.279888,
+ 01.287395, 01.294921, 01.302464, 01.310026, 01.317605,
+ 01.325203, 01.332819, 01.340454, 01.348108, 01.355780,
+ 01.363472, 01.371182, 01.378912, 01.386660, 01.394429,
+ 01.402216, 01.410024, 01.417851, 01.425698, 01.433565,
+ 01.441453, 01.449360, 01.457288, 01.465237, 01.473206,
+ 01.481196, 01.489208, 01.497240, 01.505293, 01.513368,
+ 01.521465, 01.529583, 01.537723, 01.545885, 01.554068,
+ 01.562275, 01.570503, 01.578754, 01.587028, 01.595325,
+ 01.603644, 01.611987, 01.620353, 01.628743, 01.637156,
+ 01.645593, 01.654053, 01.662538, 01.671047, 01.679581,
+ 01.688139, 01.696721, 01.705329, 01.713961, 01.722619,
+ 01.731303, 01.740011, 01.748746, 01.757506, 01.766293,
+ 01.775106, 01.783945, 01.792810, 01.801703, 01.810623,
+ 01.819569, 01.828543, 01.837545, 01.846574, 01.855631,
+ 01.864717, 01.873830, 01.882972, 01.892143, 01.901343,
+ 01.910572, 01.919830, 01.929117, 01.938434, 01.947781,
+ 01.957158, 01.966566, 01.976004, 01.985473, 01.994972,
+ 02.004503, 02.014065, 02.023659, 02.033285, 02.042943,
+ 02.052633, 02.062355, 02.072110, 02.081899, 02.091720,
+ 02.101575, 02.111464, 02.121386, 02.131343, 02.141334,
+ 02.151360, 02.161421, 02.171517, 02.181648, 02.191815,
+ 02.202018, 02.212257, 02.222533, 02.232845, 02.243195,
+ 02.253582, 02.264006, 02.274468, 02.284968, 02.295507,
+ 02.306084, 02.316701, 02.327356, 02.338051, 02.348786,
+ 02.359562, 02.370377, 02.381234, 02.392131, 02.403070,
+ 02.414051, 02.425073, 02.436138, 02.447246, 02.458397,
+ 02.469591, 02.480828, 02.492110, 02.503436, 02.514807,
+ 02.526222, 02.537684, 02.549190, 02.560743, 02.572343,
+ 02.583989, 02.595682, 02.607423, 02.619212, 02.631050,
+ 02.642936, 02.654871, 02.666855, 02.678890, 02.690975,
+ 02.703110, 02.715297, 02.727535, 02.739825, 02.752168,
+ 02.764563, 02.777012, 02.789514, 02.802070, 02.814681,
+ 02.827347, 02.840069, 02.852846, 02.865680, 02.878570,
+ 02.891518, 02.904524, 02.917588, 02.930712, 02.943894,
+ 02.957136, 02.970439, 02.983802, 02.997227, 03.010714,
+ 03.024263, 03.037875, 03.051551, 03.065290, 03.079095,
+ 03.092965, 03.106900, 03.120902, 03.134971, 03.149107,
+ 03.163312, 03.177585, 03.191928, 03.206340, 03.220824,
+ 03.235378, 03.250005, 03.264704, 03.279477, 03.294323,
+ 03.309244, 03.324240, 03.339312, 03.354461, 03.369687,
+ 03.384992, 03.400375, 03.415838, 03.431381, 03.447005,
+ 03.462711, 03.478500, 03.494372, 03.510328, 03.526370,
+ 03.542497, 03.558711, 03.575012, 03.591402, 03.607881,
+ 03.624450, 03.641111, 03.657863, 03.674708, 03.691646,
+ 03.708680, 03.725809, 03.743034, 03.760357, 03.777779,
+ 03.795300, 03.812921, 03.830645, 03.848470, 03.866400,
+ 03.884434, 03.902574, 03.920821, 03.939176, 03.957640,
+ 03.976215, 03.994901, 04.013699, 04.032612, 04.051639,
+ 04.070783, 04.090045, 04.109425, 04.128925, 04.148547,
+ 04.168292, 04.188160, 04.208154, 04.228275, 04.248524,
+ 04.268903, 04.289413, 04.310056, 04.330832, 04.351745,
+ 04.372794, 04.393982, 04.415310, 04.436781, 04.458395,
+ 04.480154, 04.502060, 04.524114, 04.546319, 04.568676,
+ 04.591187, 04.613854, 04.636678, 04.659662, 04.682807,
+ 04.706116, 04.729590, 04.753231, 04.777041, 04.801024,
+ 04.825179, 04.849511, 04.874020, 04.898710, 04.923582,
+ 04.948639, 04.973883, 04.999316, 05.024942, 05.050761,
+ 05.076778, 05.102993, 05.129411, 05.156034, 05.182864,
+ 05.209903, 05.237156, 05.264625, 05.292312, 05.320220,
+ 05.348354, 05.376714, 05.405306, 05.434131, 05.463193,
+ 05.492496, 05.522042, 05.551836, 05.581880, 05.612178,
+ 05.642734, 05.673552, 05.704634, 05.735986, 05.767610,
+ 05.799512, 05.831694, 05.864161, 05.896918, 05.929968,
+ 05.963316, 05.996967, 06.030925, 06.065194, 06.099780,
+ 06.134687, 06.169921, 06.205486, 06.241387, 06.277630,
+ 06.314220, 06.351163, 06.388465, 06.426130, 06.464166,
+ 06.502578, 06.541371, 06.580553, 06.620130, 06.660109,
+ 06.700495, 06.741297, 06.782520, 06.824173, 06.866262,
+ 06.908795, 06.951780, 06.995225, 07.039137, 07.083525,
+ 07.128398, 07.173764, 07.219632, 07.266011, 07.312910,
+ 07.360339, 07.408308, 07.456827, 07.505905, 07.555554,
+ 07.605785, 07.656608, 07.708035, 07.760077, 07.812747,
+ 07.866057, 07.920019, 07.974647, 08.029953, 08.085952,
+ 08.142657, 08.200083, 08.258245, 08.317158, 08.376837,
+ 08.437300, 08.498562, 08.560641, 08.623554, 08.687319,
+ 08.751955, 08.817481, 08.883916, 08.951282, 09.019600,
+ 09.088889, 09.159174, 09.230477, 09.302822, 09.376233,
+ 09.450735, 09.526355, 09.603118, 09.681054, 09.760191,
+ 09.840558, 09.922186, 10.005107, 10.089353, 10.174959,
+ 10.261958, 10.350389, 10.440287, 10.531693, 10.624646,
+ 10.719188, 10.815362, 10.913214, 11.012789, 11.114137,
+ 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
+ 11.762390, 11.877664, 11.995170, 12.114979, 12.237161,
+ 12.361791, 12.488946, 12.618708, 12.751161, 12.886394,
+ 13.024498, 13.165570, 13.309711, 13.457026, 13.607625,
+ 13.761625, 13.919145, 14.080314, 14.245263, 14.414134,
+ 14.587072, 14.764233, 14.945778, 15.131877, 15.322712,
+ 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
+ 16.578520, 16.808433, 17.044929, 17.288305, 17.538873,
+ 17.796967, 18.062943, 18.337176, 18.620068, 18.912049,
+ 19.213574, 19.525133, 19.847249, 20.180480, 20.525429,
+ 20.882738, 21.253102, 21.637266, 22.036036, 22.450278,
+ 22.880933, 23.329017, 23.795634, 24.281981, 24.789364,
+ 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
+ 28.365274, 29.068370, 29.808638, 30.589157, 31.413354,
+ 32.285060, 33.208568, 34.188705, 35.230920, 36.341388,
+ 37.527131, 38.796172, 40.157721, 41.622399, 43.202525,
+ 44.912465, 46.769077, 48.792279, 51.005773, 53.437996,
+ 56.123356, 59.103894
+ };
+ gnm_float X, U, V, RANLAN;
+ int I;
+
+ do {
+ X = random_01 ();
+ } while (X == 0.0);
+ U = 1000.0 * X;
+ I = U;
+ U = U - I;
+
+ if (I >= 70 && I <= 800)
+ RANLAN = F[I] + U * (F[I + 1] - F[I]);
+ else if (I >= 7 && I <= 980)
+ RANLAN = F[I] + U * (F[I + 1] - F[I] -
+ 0.25 * (1 - U) * (F[I + 2] - F[I + 1] -
+ F[I] + F[I - 1]));
+ else if (I < 7) {
+ V = gnm_log (X);
+ U = 1 / V;
+ RANLAN = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * U) * U) /
+ (1 + (3.41760202E1 + 4.01244582 * U) * U)) *
+ ( -gnm_log ( -0.91893853 - V) - 1);
+ } else {
+ U = 1 - X;
+ V = U * U;
+ if (X <= 0.999)
+ RANLAN = (1.00060006 + 2.63991156E2 *
+ U + 4.37320068E3 * V) /
+ ((1 + 2.57368075E2 * U + 3.41448018E3 * V) * U);
+ else
+ RANLAN = (1.00001538 + 6.07514119E3 * U +
+ 7.34266409E5 * V) /
+ ((1 + 6.06511919E3 * U + 6.94021044E5 * V) * U);
+ }
+
+ return RANLAN;
+}
+
+
+/* ------------------------------------------------------------------------ */
+
+/*
+ * Generate a skew-normal distributed random number.
+ *
+ * based on the information provided at
+ * http://azzalini.stat.unipd.it/SN/faq.html
+ *
+ */
+
+gnm_float
+random_skew_normal (gnm_float a)
+{
+ gnm_float result;
+ gnm_float delta = a / gnm_sqrt(1 + a * a);
+ gnm_float u = random_normal ();
+ gnm_float v = random_normal ();
+
+ result = delta * u + gnm_sqrt (1-delta*delta) * v;
+
+ return ((u < 0.) ? -result : result);
+}
+
+
+/* ------------------------------------------------------------------------ */
+
+/*
+ * Generate a skew-t distributed random number.
+ *
+ * based on the information provided at
+ * http://azzalini.stat.unipd.it/SN/faq.html
+ *
+ */
+
+gnm_float
+random_skew_tdist (gnm_float nu, gnm_float a)
+{
+ gnm_float chi = random_chisq (nu);
+ gnm_float z = random_skew_normal (a);;
+
+ return (z * gnm_sqrt(nu/chi));
+}
+
+/* ------------------------------------------------------------------------ */
diff --git a/src/gnm-random.h b/src/gnm-random.h
new file mode 100644
index 0000000..547fa05
--- /dev/null
+++ b/src/gnm-random.h
@@ -0,0 +1,44 @@
+#ifndef _GNM_RANDOM_H_
+#define _GNM_RANDOM_H_
+
+#include "numbers.h"
+
+G_BEGIN_DECLS
+
+gnm_float random_01 (void);
+gnm_float random_poisson (gnm_float lambda);
+gnm_float random_binomial (gnm_float p, gnm_float trials);
+gnm_float random_negbinom (gnm_float p, gnm_float f);
+gnm_float random_exponential (gnm_float b);
+gnm_float random_bernoulli (gnm_float p);
+gnm_float random_normal (void);
+gnm_float random_cauchy (gnm_float a);
+gnm_float random_lognormal (gnm_float zeta, gnm_float sigma);
+gnm_float random_weibull (gnm_float a, gnm_float b);
+gnm_float random_laplace (gnm_float a);
+gnm_float random_rayleigh (gnm_float sigma);
+gnm_float random_rayleigh_tail (gnm_float a, gnm_float sigma);
+gnm_float random_gamma (gnm_float a, gnm_float b);
+gnm_float random_pareto (gnm_float a, gnm_float b);
+gnm_float random_fdist (gnm_float nu1, gnm_float nu2);
+gnm_float random_beta (gnm_float a, gnm_float b);
+gnm_float random_logistic (gnm_float a);
+gnm_float random_geometric (gnm_float p);
+gnm_float random_hypergeometric (gnm_float n1, gnm_float n2, gnm_float t);
+gnm_float random_logarithmic (gnm_float p);
+gnm_float random_chisq (gnm_float nu);
+gnm_float random_tdist (gnm_float nu);
+gnm_float random_gumbel1 (gnm_float a, gnm_float b);
+gnm_float random_gumbel2 (gnm_float a, gnm_float b);
+gnm_float random_levy (gnm_float c, gnm_float alpha);
+gnm_float random_levy_skew (gnm_float c, gnm_float alpha,
+ gnm_float beta);
+gnm_float random_exppow (gnm_float a, gnm_float b);
+gnm_float random_landau (void);
+gnm_float random_gaussian_tail (gnm_float a, gnm_float sigma);
+gnm_float random_skew_normal (gnm_float a);
+gnm_float random_skew_tdist (gnm_float nu, gnm_float a);
+
+G_END_DECLS
+
+#endif
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 4224366..bfc15cd 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -7,10 +7,7 @@
* Morten Welinder <terra gnome org>
* Miguel de Icaza (miguel gnu org)
* Jukka-Pekka Iivonen (iivonen iki fi)
- * James Theiler. (See note 2.)
- * Brian Gough. (See note 2.)
- * Makoto Matsumoto and Takuji Nishimura (Mersenne Twister, see note in code)
- * Ian Smith (iandjmsmith aol com). (See note 3.)
+ * Ian Smith (iandjmsmith aol com). (See note 2.)
*/
/*
@@ -24,16 +21,7 @@
*/
/*
- * NOTE 2: most of the random distribution code comes from the GNU Scientific
- * Library (GSL), notably version 1.1.1. GSL is distributed under GPL licence,
- * see COPYING. The relevant parts are copyright (C) 1996, 1997, 1998, 1999,
- * 2000 James Theiler and Brian Gough.
- *
- * Thank you!
- */
-
-/*
- * NOTE 3: the pbeta (and support) code comes from Ian Smith. (Translated
+ * NOTE 2: the pbeta (and support) code comes from Ian Smith. (Translated
* into C, adapted to Gnumeric naming convensions, and R's API conventions
* by Morten Welinder. Blame me for problems.)
*
@@ -45,10 +33,6 @@
*/
-/* for random() */
-#define _SVID_SOURCE 1
-#define _BSD_SOURCE 1
-
#include <gnumeric-config.h>
#include "gnumeric.h"
#include "mathfunc.h"
@@ -58,7 +42,6 @@
#include <errno.h>
#include <stdlib.h>
#include <float.h>
-#include <fcntl.h>
#include <unistd.h>
#include <locale.h>
#include <string.h>
@@ -78,7 +61,6 @@
#define M_1_SQRT_2PI GNM_const(0.398942280401432677939946059934) /* 1/sqrt(2pi) */
#define M_SQRT_2dPI GNM_const(0.797884560802865355879892119869) /* sqrt(2/pi) */
#define M_2PIgnum (2 * M_PIgnum)
-#define M_Egnum GNM_const(2.718281828459045235360287471352662497757247)
#define ML_UNDERFLOW (GNM_EPSILON * GNM_EPSILON)
#define ML_ERROR(cause) /* Nothing */
@@ -6344,1355 +6326,7 @@ pbinom2 (gnm_float x0, gnm_float x1, gnm_float n, gnm_float p)
}
}
-
/* ------------------------------------------------------------------------- */
-/* http://www.math.keio.ac.jp/matumoto/CODES/MT2002/mt19937ar.c */
-/* Imported by hand -- MW. */
-
-/*
- A C-program for MT19937, with initialization improved 2002/1/26.
- Coded by Takuji Nishimura and Makoto Matsumoto.
-
- Before using, initialize the state by using init_genrand(seed)
- or init_by_array(init_key, key_length).
-
- Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
- All rights reserved.
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions
- are met:
-
- 1. Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
-
- 2. Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
-
- 3. The names of its contributors may not be used to endorse or promote
- products derived from this software without specific prior written
- permission.
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
- A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
- CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-
- Any feedback is very welcome.
- http://www.math.keio.ac.jp/matumoto/emt.html
- email: matumoto math keio ac jp
-*/
-
-#if 0
-#include <stdio.h>
-#endif
-
-/* Period parameters */
-#define N 624
-#define M 397
-#define MATRIX_A 0x9908b0dfUL /* constant vector a */
-#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
-#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
-
-static unsigned long mt[N]; /* the array for the state vector */
-static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
-
-/* initializes mt[N] with a seed */
-static void init_genrand(unsigned long s)
-{
- mt[0]= s & 0xffffffffUL;
- for (mti=1; mti<N; mti++) {
- mt[mti] =
- (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
- /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
- /* In the previous versions, MSBs of the seed affect */
- /* only MSBs of the array mt[]. */
- /* 2002/01/09 modified by Makoto Matsumoto */
- mt[mti] &= 0xffffffffUL;
- /* for >32 bit machines */
- }
-}
-
-/* initialize by an array with array-length */
-/* init_key is the array for initializing keys */
-/* key_length is its length */
-static void mt_init_by_array(unsigned long init_key[], int key_length)
-{
- int i, j, k;
- init_genrand(19650218UL);
- i=1; j=0;
- k = (N>key_length ? N : key_length);
- for (; k; k--) {
- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
- + init_key[j] + j; /* non linear */
- mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
- i++; j++;
- if (i>=N) { mt[0] = mt[N-1]; i=1; }
- if (j>=key_length) j=0;
- }
- for (k=N-1; k; k--) {
- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
- - i; /* non linear */
- mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
- i++;
- if (i>=N) { mt[0] = mt[N-1]; i=1; }
- }
-
- mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
-}
-
-/* generates a random number on [0,0xffffffff]-interval */
-static unsigned long genrand_int32(void)
-{
- unsigned long y;
- static unsigned long mag01[2]={0x0UL, MATRIX_A};
- /* mag01[x] = x * MATRIX_A for x=0,1 */
-
- if (mti >= N) { /* generate N words at one time */
- int kk;
-
- if (mti == N+1) /* if init_genrand() has not been called, */
- init_genrand(5489UL); /* a default initial seed is used */
-
- for (kk=0;kk<N-M;kk++) {
- y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
- mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
- }
- for (;kk<N-1;kk++) {
- y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
- mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
- }
- y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
- mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
-
- mti = 0;
- }
-
- y = mt[mti++];
-
- /* Tempering */
- y ^= (y >> 11);
- y ^= (y << 7) & 0x9d2c5680UL;
- y ^= (y << 15) & 0xefc60000UL;
- y ^= (y >> 18);
-
- return y;
-}
-
-#if 0
-/* generates a random number on [0,0x7fffffff]-interval */
-long genrand_int31(void)
-{
- return (long)(genrand_int32()>>1);
-}
-
-/* generates a random number on [0,1]-real-interval */
-double genrand_real1(void)
-{
- return genrand_int32()*(1.0/4294967295.0);
- /* divided by 2^32-1 */
-}
-
-/* generates a random number on [0,1)-real-interval */
-double genrand_real2(void)
-{
- return genrand_int32()*(1.0/4294967296.0);
- /* divided by 2^32 */
-}
-
-/* generates a random number on (0,1)-real-interval */
-double genrand_real3(void)
-{
- return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
- /* divided by 2^32 */
-}
-#endif
-
-/* generates a random number on [0,1) with 53-bit resolution*/
-static double genrand_res53(void)
-{
- unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
- return(a*67108864.0+b)*(1.0/9007199254740992.0);
-}
-/* These real versions are due to Isaku Wada, 2002/01/09 added */
-
-#if 0
-int main(void)
-{
- int i;
- unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
- init_by_array(init, length);
- printf("1000 outputs of genrand_int32()\n");
- for (i=0; i<1000; i++) {
- printf("%10lu ", genrand_int32());
- if (i%5==4) printf("\n");
- }
- printf("\n1000 outputs of genrand_real2()\n");
- for (i=0; i<1000; i++) {
- printf("%10.8f ", genrand_real2());
- if (i%5==4) printf("\n");
- }
- return 0;
-}
-#endif
-
-
-#undef N
-#undef M
-#undef MATRIX_A
-#undef UPPER_MASK
-#undef LOWER_MASK
-
-/* ------------------------------------------------------------------------ */
-
-/* FIXME: we need something that catches partials and EAGAIN. */
-#define fullread read
-
-#define RANDOM_DEVICE "/dev/urandom"
-
-/*
- * Conservative random number generator. The result is (supposedly) uniform
- * and between 0 and 1. (0 possible, 1 not.) The result should have about
- * 64 bits randomness.
- */
-gnm_float
-random_01 (void)
-{
- static int device_fd = -2;
- static int seeded = -2;
-
- while (seeded) {
- if (seeded == -2) {
- char const *seed = g_getenv ("GNUMERIC_PRNG_SEED");
- if (seed) {
- int len = strlen (seed);
- int i;
- unsigned long *longs = g_new (unsigned long, len + 1);
-
- /* We drop only one character into each long. */
- for (i = 0; i < len; i++)
- longs[i] = (unsigned char)seed[i];
- mt_init_by_array (longs, len);
- g_free (longs);
- seeded = 1;
-
- g_warning ("Using pseudo-random numbers.");
- } else {
- seeded = 0;
- break;
- }
- }
-
- /*
- * Only 52-bit precision. But hey, if you are using pseudo
- * random numbers that ought to be good enough to you.
- */
- return genrand_res53 ();
- }
-
- if (device_fd == -2) {
- device_fd = g_open (RANDOM_DEVICE, O_RDONLY, 0);
- /*
- * We could check that we really have a device, but it hard
- * to come up with a non-paranoid reason to.
- */
- }
-
- if (device_fd >= 0) {
- static ssize_t bytes_left = 0;
- static unsigned char data[32 * sizeof (gnm_float)];
- gnm_float res = 0;
- size_t i;
-
- if (bytes_left < (ssize_t)sizeof (gnm_float)) {
- ssize_t got = fullread (device_fd, &data, sizeof (data));
- if (got < (ssize_t)sizeof (gnm_float))
- goto failure;
- bytes_left += got;
- }
-
- bytes_left -= sizeof (gnm_float);
- for (i = 0; i < sizeof (gnm_float); i++)
- res = (res + data[bytes_left + i]) / 256;
- return res;
-
- failure:
- /* It failed when it shouldn't. Disable. */
- g_warning ("Reading from %s failed; reverting to pseudo-random.",
- RANDOM_DEVICE);
- close (device_fd);
- device_fd = -1;
- }
-
- return genrand_res53 ();
-}
-
-/*
- * Generate a N(0,1) distributed number.
- */
-gnm_float
-random_normal (void)
-{
- static gboolean has_saved = FALSE;
- static gnm_float saved;
-
- if (has_saved) {
- has_saved = FALSE;
- return saved;
- } else {
- gnm_float u, v, r2, rsq;
- do {
- u = 2 * random_01 () - 1;
- v = 2 * random_01 () - 1;
- r2 = u * u + v * v;
- } while (r2 > 1 || r2 == 0);
-
- rsq = gnm_sqrt (-2 * gnm_log (r2) / r2);
-
- has_saved = TRUE;
- saved = v * rsq;
-
- return u * rsq;
- }
-}
-
-gnm_float
-random_lognormal (gnm_float zeta, gnm_float sigma)
-{
- return gnm_exp (sigma * random_normal () + zeta);
-}
-
-static gnm_float
-random_gaussian (gnm_float sigma)
-{
- return sigma * random_normal ();
-}
-
-/*
- * Generate a poisson distributed number.
- */
-gnm_float
-random_poisson (gnm_float lambda)
-{
- /*
- * This may not be optimal code, but it sure is easy to
- * understand compared to R's code.
- */
- return qpois (random_01 (), lambda, TRUE, FALSE);
-}
-
-/*
- * Generate a binomial distributed number.
- */
-gnm_float
-random_binomial (gnm_float p, gnm_float trials)
-{
- return qbinom (random_01 (), trials, p, TRUE, FALSE);
-}
-
-/*
- * Generate a negative binomial distributed number.
- */
-gnm_float
-random_negbinom (gnm_float p, gnm_float f)
-{
- return qnbinom (random_01 (), f, p, TRUE, FALSE);
-}
-
-/*
- * Generate an exponential distributed number.
- */
-gnm_float
-random_exponential (gnm_float b)
-{
- return -b * gnm_log (random_01 ());
-}
-
-/*
- * Generate a bernoulli distributed number.
- */
-gnm_float
-random_bernoulli (gnm_float p)
-{
- gnm_float r = random_01 ();
-
- return (r <= p) ? 1.0 : 0.0;
-}
-
-/*
- * Generate a cauchy distributed number. From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_cauchy (gnm_float a)
-{
- gnm_float u;
-
- do {
- u = random_01 ();
- } while (u == 0.5);
-
- return a * gnm_tan (M_PIgnum * u);
-}
-
-/*
- * Generate a Weibull distributed number. From the GNU Scientific
- * library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_weibull (gnm_float a, gnm_float b)
-{
- gnm_float x, z;
-
- do {
- x = random_01 ();
- } while (x == 0.0);
-
- z = gnm_pow (-gnm_log (x), 1 / b);
-
- return a * z;
-}
-
-/*
- * Generate a Laplace (two-sided exponential probability) distributed number.
- * From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_laplace (gnm_float a)
-{
- gnm_float u;
-
- do {
- u = 2 * random_01 () - 1.0;
- } while (u == 0.0);
-
- if (u < 0)
- return a * gnm_log (-u);
- else
- return -a * gnm_log (u);
-}
-
-gnm_float
-random_laplace_pdf (gnm_float x, gnm_float a)
-{
- return (1 / (2 * a)) * gnm_exp (-gnm_abs (x) / a);
-}
-
-/*
- * Generate a Rayleigh distributed number. From the GNU Scientific library
- * 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_rayleigh (gnm_float sigma)
-{
- gnm_float u;
-
- do {
- u = random_01 ();
- } while (u == 0.0);
-
- return sigma * gnm_sqrt (-2.0 * gnm_log (u));
-}
-
-/*
- * Generate a Rayleigh tail distributed number. From the GNU Scientific library
- * 1.1.1. The Rayleigh tail distribution has the form
- * p(x) dx = (x / sigma^2) exp((a^2 - x^2)/(2 sigma^2)) dx
- *
- * for x = a ... +infty
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_rayleigh_tail (gnm_float a, gnm_float sigma)
-{
- gnm_float u;
-
- do {
- u = random_01 ();
- } while (u == 0.0);
-
- return gnm_sqrt (a * a - 2.0 * sigma * sigma * gnm_log (u));
-}
-
-/* The Gamma distribution of order a>0 is defined by:
- *
- * p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx
- *
- * for x>0. If X and Y are independent gamma-distributed random
- * variables of order a1 and a2 with the same scale parameter b, then
- * X+Y has gamma distribution of order a1+a2.
- *
- * The algorithms below are from Knuth, vol 2, 2nd ed, p. 129.
- */
-
-static gnm_float
-gamma_frac (gnm_float a)
-{
- /* This is exercise 16 from Knuth; see page 135, and the solution is
- * on page 551. */
-
- gnm_float x, q;
- gnm_float p = M_Egnum / (a + M_Egnum);
- do {
- gnm_float v;
- gnm_float u = random_01 ();
- do {
- v = random_01 ();
- } while (v == 0.0);
-
- if (u < p) {
- x = gnm_pow (v, 1 / a);
- q = gnm_exp (-x);
- } else {
- x = 1 - gnm_log (v);
- q = gnm_pow (x, a - 1);
- }
- } while (random_01 () >= q);
-
- return x;
-}
-
-static gnm_float
-gamma_large (gnm_float a)
-{
- /*
- * Works only if a > 1, and is most efficient if a is large
- *
- * This algorithm, reported in Knuth, is attributed to Ahrens. A
- * faster one, we are told, can be found in: J. H. Ahrens and
- * U. Dieter, Computing 12 (1974) 223-246.
- */
-
- gnm_float sqa, x, y, v;
- sqa = gnm_sqrt (2 * a - 1);
- do {
- do {
- y = gnm_tan (M_PIgnum * random_01 ());
- x = sqa * y + a - 1;
- } while (x <= 0);
- v = random_01 ();
- } while (v > (1 + y * y) * gnm_exp ((a - 1) * gnm_log (x / (a - 1)) -
- sqa * y));
-
- return x;
-}
-
-static gnm_float
-ran_gamma_int (gnm_float a)
-{
- if (a < 12) {
- gnm_float prod;
-
- do {
- unsigned int i, ua;
- prod = 1;
- ua = (unsigned int)a;
-
- for (i = 0; i < ua; i++)
- prod *= random_01 ();
-
- /*
- * This handles the 0-probability event of getting
- * an actual zero as well as the possibility of
- * underflow.
- */
- } while (prod == 0);
-
- return -gnm_log (prod);
- } else
- return gamma_large (a);
-}
-
-/*
- * Generate a Gamma distributed number. From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_gamma (gnm_float a, gnm_float b)
-{
- gnm_float na;
-
- if (gnm_isnan (a) || gnm_isnan (b) || a <= 0 || b <= 0)
- return gnm_nan;
-
- na = gnm_floor (a);
-
- if (a == na)
- return b * ran_gamma_int (na);
- else if (na == 0)
- return b * gamma_frac (a);
- else
- return b * (ran_gamma_int (na) + gamma_frac (a - na));
-}
-
-/*
- * Generate a Pareto distributed number. From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_pareto (gnm_float a, gnm_float b)
-{
- gnm_float x;
-
- do {
- x = random_01 ();
- } while (x == 0.0);
-
- return b * gnm_pow (x, -1 / a);
-}
-
-/*
- * Generate a F-distributed number. From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_fdist (gnm_float nu1, gnm_float nu2)
-{
- gnm_float Y1 = random_gamma (nu1 / 2, 2.0);
- gnm_float Y2 = random_gamma (nu2 / 2, 2.0);
-
- return (Y1 * nu2) / (Y2 * nu1);
-}
-
-/*
- * Generate a Beta-distributed number. From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_beta (gnm_float a, gnm_float b)
-{
- gnm_float x1 = random_gamma (a, 1.0);
- gnm_float x2 = random_gamma (b, 1.0);
-
- return x1 / (x1 + x2);
-}
-
-/*
- * Generate a Chi-Square-distributed number. From the GNU Scientific library
- * 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_chisq (gnm_float nu)
-{
- return 2 * random_gamma (nu / 2, 1.0);
-}
-
-/*
- * Generate a logistic-distributed number. From the GNU Scientific library
- * 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_logistic (gnm_float a)
-{
- gnm_float x;
-
- do {
- x = random_01 ();
- } while (x == 0 || x == 1);
-
- return a * gnm_log (x / (1 - x));
-}
-
-/*
- * Generate a geometric-distributed number. From the GNU Scientific library
- * 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_geometric (gnm_float p)
-{
- gnm_float u;
-
- if (p == 1)
- return 1;
- do {
- u = random_01 ();
- } while (u == 0);
-
- return gnm_floor (gnm_log (u) / gnm_log1p (-p) + 1);
-}
-
-gnm_float
-random_hypergeometric (gnm_float n1, gnm_float n2, gnm_float t)
-{
- return qhyper (random_01 (), n1, n2, t, TRUE, FALSE);
-}
-
-
-/*
- * Generate a logarithmic-distributed number. From the GNU Scientific library
- * 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_logarithmic (gnm_float p)
-{
- gnm_float c, v;
-
- c = gnm_log1p (-p);
- do {
- v = random_01 ();
- } while (v == 0);
-
- if (v >= p)
- return 1;
- else {
- gnm_float u, q;
-
- do {
- u = random_01 ();
- } while (u == 0);
- q = gnm_expm1 (c * u);
-
- if (v <= q * q)
- return gnm_floor (1 + gnm_log (v) / gnm_log (q));
- else if (v <= q)
- return 2;
- else
- return 1;
- }
-}
-
-/*
- * Generate a T-distributed number. From the GNU Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_tdist (gnm_float nu)
-{
- if (nu <= 2) {
- gnm_float Y1 = random_normal ();
- gnm_float Y2 = random_chisq (nu);
-
- gnm_float t = Y1 / gnm_sqrt (Y2 / nu);
-
- return t;
- } else {
- gnm_float Y1, Y2, Z, t;
- do {
- Y1 = random_normal ();
- Y2 = random_exponential (1 / (nu / 2 - 1));
-
- Z = Y1 * Y1 / (nu - 2);
- } while (1 - Z < 0 || gnm_exp (-Y2 - Z) > (1 - Z));
-
- /* Note that there is a typo in Knuth's formula, the line below
- * is taken from the original paper of Marsaglia, Mathematics
- * of Computation, 34 (1980), p 234-256. */
-
- t = Y1 / gnm_sqrt ((1 - 2 / nu) * (1 - Z));
- return t;
- }
-}
-
-/*
- * Generate a Type I Gumbel-distributed random number. From the GNU
- * Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_gumbel1 (gnm_float a, gnm_float b)
-{
- gnm_float x;
-
- do {
- x = random_01 ();
- } while (x == 0.0);
-
- return (gnm_log (b) - gnm_log (-gnm_log (x))) / a;
-}
-
-/*
- * Generate a Type II Gumbel-distributed random number. From the GNU
- * Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- */
-gnm_float
-random_gumbel2 (gnm_float a, gnm_float b)
-{
- gnm_float x;
-
- do {
- x = random_01 ();
- } while (x == 0.0);
-
- return gnm_pow (-b / gnm_log (x), 1 / a);
-}
-
-/*
- * Generate a stable Levy-distributed random number. From the GNU
- * Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
- *
- * The stable Levy probability distributions have the form
- *
- * p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
- *
- * with 0 < alpha <= 2.
- *
- * For alpha = 1, we get the Cauchy distribution
- * For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.
- *
- * Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
- * Simulation". The original reference given there is,
- *
- * J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
- * simulating stable random variates". Journal of the American
- * Statistical Association, JASA 71 340-344 (1976).
- */
-gnm_float
-random_levy (gnm_float c, gnm_float alpha)
-{
- gnm_float u, v, t, s;
-
- do {
- u = random_01 ();
- } while (u == 0.0);
-
- u = M_PIgnum * (u - 0.5);
-
- if (alpha == 1) { /* cauchy case */
- t = gnm_tan (u);
- return c * t;
- }
-
- do {
- v = random_exponential (1.0);
- } while (v == 0);
-
- if (alpha == 2) { /* gaussian case */
- t = 2 * gnm_sin (u) * gnm_sqrt (v);
- return c * t;
- }
-
- /* general case */
-
- t = gnm_sin (alpha * u) / gnm_pow (gnm_cos (u), 1 / alpha);
- s = gnm_pow (gnm_cos ((1 - alpha) * u) / v, (1 - alpha) / alpha);
-
- return c * t * s;
-}
-
-/*
- * The following routine for the skew-symmetric case was provided by
- * Keith Briggs.
- *
- * The stable Levy probability distributions have the form
- *
- * 2*pi* p(x) dx
- *
- * = int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for
- * alpha != 1
- * = int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|))) for
- alpha == 1
- *
- * with 0<alpha<=2, -1<=beta<=1, sigma>0.
- *
- * For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
- *
- * For alpha = 1, beta=0, we get the Lorentz distribution
- * For alpha = 2, beta=0, we get the Gaussian distribution
- *
- * See A. Weron and R. Weron: Computer simulation of Lévy alpha-stable
- * variables and processes, preprint Technical University of Wroclaw.
- * http://www.im.pwr.wroc.pl/~hugo/Publications.html
- */
-gnm_float
-random_levy_skew (gnm_float c, gnm_float alpha, gnm_float beta)
-{
- gnm_float V, W, X;
-
- if (beta == 0) /* symmetric case */
- return random_levy (c, alpha);
-
- do {
- V = random_01 ();
- } while (V == 0.0);
-
- V = M_PIgnum * (V - 0.5);
-
- do {
- W = random_exponential (1.0);
- } while (W == 0);
-
- if (alpha == 1) {
- X = ((M_PI_2gnum + beta * V) * gnm_tan (V) -
- beta * gnm_log (M_PI_2gnum * W * gnm_cos (V) /
- (M_PI_2gnum + beta * V))) / M_PI_2gnum;
- return c * (X + beta * gnm_log (c) / M_PI_2gnum);
- } else {
- gnm_float t = beta * gnm_tan (M_PI_2gnum * alpha);
- gnm_float B = gnm_atan (t) / alpha;
- gnm_float S = pow1p (t * t, 1 / (2 * alpha));
-
- X = S * gnm_sin (alpha * (V + B)) / gnm_pow (gnm_cos (V),
- 1 / alpha)
- * gnm_pow (gnm_cos (V - alpha * (V + B)) / W,
- (1 - alpha) / alpha);
- return c * X;
- }
-}
-
-gnm_float
-random_exppow_pdf (gnm_float x, gnm_float a, gnm_float b)
-{
- gnm_float lngamma = lgamma1p (1 / b);
-
- return (1 / (2 * a)) * gnm_exp (-gnm_pow (gnm_abs (x / a), b) - lngamma);
-}
-
-/*
- * The exponential power probability distribution is
- *
- * p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx
- *
- * for -infty < x < infty. For b = 1 it reduces to the Laplace
- * distribution.
- *
- * The exponential power distribution is related to the gamma
- * distribution by E = a * pow(G(1/b),1/b), where E is an exponential
- * power variate and G is a gamma variate.
- *
- * We use this relation for b < 1. For b >=1 we use rejection methods
- * based on the laplace and gaussian distributions which should be
- * faster.
- *
- * See P. R. Tadikamalla, "Random Sampling from the Exponential Power
- * Distribution", Journal of the American Statistical Association,
- * September 1980, Volume 75, Number 371, pages 683-686.
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
- */
-
-gnm_float
-random_exppow (gnm_float a, gnm_float b)
-{
- /* See http://www.mcgill.ca/files/economics/propertiesandestimation.pdf */
- if (!(a > 0) || gnm_isnan (b))
- return gnm_nan;
-
- if (b < 1) {
- gnm_float u = random_01 ();
- gnm_float v = random_gamma (1 / b, 1.0);
- gnm_float z = a * gnm_pow (v, 1 / b) ;
-
- if (u > 0.5)
- return z;
- else
- return -z;
- } else if (b == 1)
- return random_laplace (a); /* Laplace distribution */
- else if (b < 2) {
- /* Use laplace distribution for rejection method */
- gnm_float x, y, h, ratio, u;
-
- /* Scale factor chosen by upper bound on ratio at b = 2 */
- gnm_float s = 1.4489;
- do {
- x = random_laplace (a);
- y = random_laplace_pdf (x, a);
- h = random_exppow_pdf (x, a, b);
- ratio = h / (s * y);
- u = random_01 ();
- } while (u > ratio);
-
- return x ;
- } else if (b == 2) /* Gaussian distribution */
- return random_gaussian (a / gnm_sqrt (2.0));
- else {
- /* Use gaussian for rejection method */
- gnm_float x, y, h, ratio, u;
- const gnm_float sigma = a / gnm_sqrt (2.0);
-
- /* Scale factor chosen by upper bound on ratio at b = infinity.
- * This could be improved by using a rational function
- * approximation to the bounding curve. */
-
- gnm_float s = 2.4091 ; /* this is sqrt(pi) e / 2 */
-
- do {
- x = random_gaussian (sigma) ;
- y = dnorm (x, 0.0, gnm_abs (sigma), FALSE) ;
- h = random_exppow_pdf (x, a, b) ;
- ratio = h / (s * y) ;
- u = random_01 ();
- } while (u > ratio);
-
- return x;
- }
-}
-
-/*
- * Generate a Gaussian tail-distributed random number. From the GNU
- * Scientific library 1.1.1.
- * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
- */
-gnm_float
-random_gaussian_tail (gnm_float a, gnm_float sigma)
-{
- /*
- * Returns a gaussian random variable larger than a
- * This implementation does one-sided upper-tailed deviates.
- */
-
- gnm_float s = a / sigma;
-
- if (s < 1) {
- /* For small s, use a direct rejection method. The limit s < 1
- * can be adjusted to optimise the overall efficiency */
-
- gnm_float x;
-
- do {
- x = random_gaussian (1.0);
- } while (x < s);
- return x * sigma;
- } else {
- /* Use the "supertail" deviates from the last two steps
- * of Marsaglia's rectangle-wedge-tail method, as described
- * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11,
- * p139, and the solution, p586.)
- */
-
- gnm_float u, v, x;
-
- do {
- u = random_01 ();
- do {
- v = random_01 ();
- } while (v == 0.0);
- x = gnm_sqrt (s * s - 2 * gnm_log (v));
- } while (x * u > s);
- return x * sigma;
- }
-}
-
-/*
- * Generate a Landau-distributed random number. From the GNU Scientific
- * library 1.1.1.
- *
- * Copyright (C) 2001 David Morrison
- *
- * Adapted from the CERN library routines DENLAN, RANLAN, and DISLAN
- * as described in http://consult.cern.ch/shortwrups/g110/top.html.
- * Original author: K.S. K\"olbig.
- *
- * The distribution is given by the complex path integral,
- *
- * p(x) = (1/(2 pi i)) \int_{c-i\inf}^{c+i\inf} ds exp(s log(s) + x s)
- *
- * which can be converted into a real integral over [0,+\inf]
- *
- * p(x) = (1/pi) \int_0^\inf dt \exp(-t log(t) - x t) sin(pi t)
- */
-
-gnm_float
-random_landau (void)
-{
- static gnm_float F[983] = {
- 0.0000000, /*
- * Add empty element [0] to account for difference
- * between C and Fortran convention for lower bound.
- */
- 00.000000, 00.000000, 00.000000, 00.000000, 00.000000,
- -2.244733, -2.204365, -2.168163, -2.135219, -2.104898,
- -2.076740, -2.050397, -2.025605, -2.002150, -1.979866,
- -1.958612, -1.938275, -1.918760, -1.899984, -1.881879,
- -1.864385, -1.847451, -1.831030, -1.815083, -1.799574,
- -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
- -1.714187, -1.701029, -1.688130, -1.675477, -1.663057,
- -1.650858, -1.638868, -1.627078, -1.615477, -1.604058,
- -1.592811, -1.581729, -1.570806, -1.560034, -1.549407,
- -1.538919, -1.528565, -1.518339, -1.508237, -1.498254,
- -1.488386, -1.478628, -1.468976, -1.459428, -1.449979,
- -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
- -1.395194, -1.386356, -1.377594, -1.368906, -1.360291,
- -1.351746, -1.343269, -1.334859, -1.326512, -1.318229,
- -1.310006, -1.301843, -1.293737, -1.285688, -1.277693,
- -1.269752, -1.261863, -1.254024, -1.246235, -1.238494,
- -1.230800, -1.223153, -1.215550, -1.207990, -1.200474,
- -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
- -1.156220, -1.148977, -1.141770, -1.134598, -1.127459,
- -1.120354, -1.113282, -1.106242, -1.099233, -1.092255,
- -1.085306, -1.078388, -1.071498, -1.064636, -1.057802,
- -1.050996, -1.044215, -1.037461, -1.030733, -1.024029,
- -1.017350, -1.010695, -1.004064, -0.997456, -0.990871,
- -0.984308, -0.977767, -0.971247, -0.964749, -0.958271,
- -0.951813, -0.945375, -0.938957, -0.932558, -0.926178,
- -0.919816, -0.913472, -0.907146, -0.900838, -0.894547,
- -0.888272, -0.882014, -0.875773, -0.869547, -0.863337,
- -0.857142, -0.850963, -0.844798, -0.838648, -0.832512,
- -0.826390, -0.820282, -0.814187, -0.808106, -0.802038,
- -0.795982, -0.789940, -0.783909, -0.777891, -0.771884,
- -0.765889, -0.759906, -0.753934, -0.747973, -0.742023,
- -0.736084, -0.730155, -0.724237, -0.718328, -0.712429,
- -0.706541, -0.700661, -0.694791, -0.688931, -0.683079,
- -0.677236, -0.671402, -0.665576, -0.659759, -0.653950,
- -0.648149, -0.642356, -0.636570, -0.630793, -0.625022,
- -0.619259, -0.613503, -0.607754, -0.602012, -0.596276,
- -0.590548, -0.584825, -0.579109, -0.573399, -0.567695,
- -0.561997, -0.556305, -0.550618, -0.544937, -0.539262,
- -0.533592, -0.527926, -0.522266, -0.516611, -0.510961,
- -0.505315, -0.499674, -0.494037, -0.488405, -0.482777,
- -0.477153, -0.471533, -0.465917, -0.460305, -0.454697,
- -0.449092, -0.443491, -0.437893, -0.432299, -0.426707,
- -0.421119, -0.415534, -0.409951, -0.404372, -0.398795,
- -0.393221, -0.387649, -0.382080, -0.376513, -0.370949,
- -0.365387, -0.359826, -0.354268, -0.348712, -0.343157,
- -0.337604, -0.332053, -0.326503, -0.320955, -0.315408,
- -0.309863, -0.304318, -0.298775, -0.293233, -0.287692,
- -0.282152, -0.276613, -0.271074, -0.265536, -0.259999,
- -0.254462, -0.248926, -0.243389, -0.237854, -0.232318,
- -0.226783, -0.221247, -0.215712, -0.210176, -0.204641,
- -0.199105, -0.193568, -0.188032, -0.182495, -0.176957,
- -0.171419, -0.165880, -0.160341, -0.154800, -0.149259,
- -0.143717, -0.138173, -0.132629, -0.127083, -0.121537,
- -0.115989, -0.110439, -0.104889, -0.099336, -0.093782,
- -0.088227, -0.082670, -0.077111, -0.071550, -0.065987,
- -0.060423, -0.054856, -0.049288, -0.043717, -0.038144,
- -0.032569, -0.026991, -0.021411, -0.015828, -0.010243,
- -0.004656, 00.000934, 00.006527, 00.012123, 00.017722,
- 00.023323, 00.028928, 00.034535, 00.040146, 00.045759,
- 00.051376, 00.056997, 00.062620, 00.068247, 00.073877,
- 00.079511, 00.085149, 00.090790, 00.096435, 00.102083,
- 00.107736, 00.113392, 00.119052, 00.124716, 00.130385,
- 00.136057, 00.141734, 00.147414, 00.153100, 00.158789,
- 00.164483, 00.170181, 00.175884, 00.181592, 00.187304,
- 00.193021, 00.198743, 00.204469, 00.210201, 00.215937,
- 00.221678, 00.227425, 00.233177, 00.238933, 00.244696,
- 00.250463, 00.256236, 00.262014, 00.267798, 00.273587,
- 00.279382, 00.285183, 00.290989, 00.296801, 00.302619,
- 00.308443, 00.314273, 00.320109, 00.325951, 00.331799,
- 00.337654, 00.343515, 00.349382, 00.355255, 00.361135,
- 00.367022, 00.372915, 00.378815, 00.384721, 00.390634,
- 00.396554, 00.402481, 00.408415, 00.414356, 00.420304,
- 00.426260, 00.432222, 00.438192, 00.444169, 00.450153,
- 00.456145, 00.462144, 00.468151, 00.474166, 00.480188,
- 00.486218, 00.492256, 00.498302, 00.504356, 00.510418,
- 00.516488, 00.522566, 00.528653, 00.534747, 00.540850,
- 00.546962, 00.553082, 00.559210, 00.565347, 00.571493,
- 00.577648, 00.583811, 00.589983, 00.596164, 00.602355,
- 00.608554, 00.614762, 00.620980, 00.627207, 00.633444,
- 00.639689, 00.645945, 00.652210, 00.658484, 00.664768,
- 00.671062, 00.677366, 00.683680, 00.690004, 00.696338,
- 00.702682, 00.709036, 00.715400, 00.721775, 00.728160,
- 00.734556, 00.740963, 00.747379, 00.753807, 00.760246,
- 00.766695, 00.773155, 00.779627, 00.786109, 00.792603,
- 00.799107, 00.805624, 00.812151, 00.818690, 00.825241,
- 00.831803, 00.838377, 00.844962, 00.851560, 00.858170,
- 00.864791, 00.871425, 00.878071, 00.884729, 00.891399,
- 00.898082, 00.904778, 00.911486, 00.918206, 00.924940,
- 00.931686, 00.938446, 00.945218, 00.952003, 00.958802,
- 00.965614, 00.972439, 00.979278, 00.986130, 00.992996,
- 00.999875, 01.006769, 01.013676, 01.020597, 01.027533,
- 01.034482, 01.041446, 01.048424, 01.055417, 01.062424,
- 01.069446, 01.076482, 01.083534, 01.090600, 01.097681,
- 01.104778, 01.111889, 01.119016, 01.126159, 01.133316,
- 01.140490, 01.147679, 01.154884, 01.162105, 01.169342,
- 01.176595, 01.183864, 01.191149, 01.198451, 01.205770,
- 01.213105, 01.220457, 01.227826, 01.235211, 01.242614,
- 01.250034, 01.257471, 01.264926, 01.272398, 01.279888,
- 01.287395, 01.294921, 01.302464, 01.310026, 01.317605,
- 01.325203, 01.332819, 01.340454, 01.348108, 01.355780,
- 01.363472, 01.371182, 01.378912, 01.386660, 01.394429,
- 01.402216, 01.410024, 01.417851, 01.425698, 01.433565,
- 01.441453, 01.449360, 01.457288, 01.465237, 01.473206,
- 01.481196, 01.489208, 01.497240, 01.505293, 01.513368,
- 01.521465, 01.529583, 01.537723, 01.545885, 01.554068,
- 01.562275, 01.570503, 01.578754, 01.587028, 01.595325,
- 01.603644, 01.611987, 01.620353, 01.628743, 01.637156,
- 01.645593, 01.654053, 01.662538, 01.671047, 01.679581,
- 01.688139, 01.696721, 01.705329, 01.713961, 01.722619,
- 01.731303, 01.740011, 01.748746, 01.757506, 01.766293,
- 01.775106, 01.783945, 01.792810, 01.801703, 01.810623,
- 01.819569, 01.828543, 01.837545, 01.846574, 01.855631,
- 01.864717, 01.873830, 01.882972, 01.892143, 01.901343,
- 01.910572, 01.919830, 01.929117, 01.938434, 01.947781,
- 01.957158, 01.966566, 01.976004, 01.985473, 01.994972,
- 02.004503, 02.014065, 02.023659, 02.033285, 02.042943,
- 02.052633, 02.062355, 02.072110, 02.081899, 02.091720,
- 02.101575, 02.111464, 02.121386, 02.131343, 02.141334,
- 02.151360, 02.161421, 02.171517, 02.181648, 02.191815,
- 02.202018, 02.212257, 02.222533, 02.232845, 02.243195,
- 02.253582, 02.264006, 02.274468, 02.284968, 02.295507,
- 02.306084, 02.316701, 02.327356, 02.338051, 02.348786,
- 02.359562, 02.370377, 02.381234, 02.392131, 02.403070,
- 02.414051, 02.425073, 02.436138, 02.447246, 02.458397,
- 02.469591, 02.480828, 02.492110, 02.503436, 02.514807,
- 02.526222, 02.537684, 02.549190, 02.560743, 02.572343,
- 02.583989, 02.595682, 02.607423, 02.619212, 02.631050,
- 02.642936, 02.654871, 02.666855, 02.678890, 02.690975,
- 02.703110, 02.715297, 02.727535, 02.739825, 02.752168,
- 02.764563, 02.777012, 02.789514, 02.802070, 02.814681,
- 02.827347, 02.840069, 02.852846, 02.865680, 02.878570,
- 02.891518, 02.904524, 02.917588, 02.930712, 02.943894,
- 02.957136, 02.970439, 02.983802, 02.997227, 03.010714,
- 03.024263, 03.037875, 03.051551, 03.065290, 03.079095,
- 03.092965, 03.106900, 03.120902, 03.134971, 03.149107,
- 03.163312, 03.177585, 03.191928, 03.206340, 03.220824,
- 03.235378, 03.250005, 03.264704, 03.279477, 03.294323,
- 03.309244, 03.324240, 03.339312, 03.354461, 03.369687,
- 03.384992, 03.400375, 03.415838, 03.431381, 03.447005,
- 03.462711, 03.478500, 03.494372, 03.510328, 03.526370,
- 03.542497, 03.558711, 03.575012, 03.591402, 03.607881,
- 03.624450, 03.641111, 03.657863, 03.674708, 03.691646,
- 03.708680, 03.725809, 03.743034, 03.760357, 03.777779,
- 03.795300, 03.812921, 03.830645, 03.848470, 03.866400,
- 03.884434, 03.902574, 03.920821, 03.939176, 03.957640,
- 03.976215, 03.994901, 04.013699, 04.032612, 04.051639,
- 04.070783, 04.090045, 04.109425, 04.128925, 04.148547,
- 04.168292, 04.188160, 04.208154, 04.228275, 04.248524,
- 04.268903, 04.289413, 04.310056, 04.330832, 04.351745,
- 04.372794, 04.393982, 04.415310, 04.436781, 04.458395,
- 04.480154, 04.502060, 04.524114, 04.546319, 04.568676,
- 04.591187, 04.613854, 04.636678, 04.659662, 04.682807,
- 04.706116, 04.729590, 04.753231, 04.777041, 04.801024,
- 04.825179, 04.849511, 04.874020, 04.898710, 04.923582,
- 04.948639, 04.973883, 04.999316, 05.024942, 05.050761,
- 05.076778, 05.102993, 05.129411, 05.156034, 05.182864,
- 05.209903, 05.237156, 05.264625, 05.292312, 05.320220,
- 05.348354, 05.376714, 05.405306, 05.434131, 05.463193,
- 05.492496, 05.522042, 05.551836, 05.581880, 05.612178,
- 05.642734, 05.673552, 05.704634, 05.735986, 05.767610,
- 05.799512, 05.831694, 05.864161, 05.896918, 05.929968,
- 05.963316, 05.996967, 06.030925, 06.065194, 06.099780,
- 06.134687, 06.169921, 06.205486, 06.241387, 06.277630,
- 06.314220, 06.351163, 06.388465, 06.426130, 06.464166,
- 06.502578, 06.541371, 06.580553, 06.620130, 06.660109,
- 06.700495, 06.741297, 06.782520, 06.824173, 06.866262,
- 06.908795, 06.951780, 06.995225, 07.039137, 07.083525,
- 07.128398, 07.173764, 07.219632, 07.266011, 07.312910,
- 07.360339, 07.408308, 07.456827, 07.505905, 07.555554,
- 07.605785, 07.656608, 07.708035, 07.760077, 07.812747,
- 07.866057, 07.920019, 07.974647, 08.029953, 08.085952,
- 08.142657, 08.200083, 08.258245, 08.317158, 08.376837,
- 08.437300, 08.498562, 08.560641, 08.623554, 08.687319,
- 08.751955, 08.817481, 08.883916, 08.951282, 09.019600,
- 09.088889, 09.159174, 09.230477, 09.302822, 09.376233,
- 09.450735, 09.526355, 09.603118, 09.681054, 09.760191,
- 09.840558, 09.922186, 10.005107, 10.089353, 10.174959,
- 10.261958, 10.350389, 10.440287, 10.531693, 10.624646,
- 10.719188, 10.815362, 10.913214, 11.012789, 11.114137,
- 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
- 11.762390, 11.877664, 11.995170, 12.114979, 12.237161,
- 12.361791, 12.488946, 12.618708, 12.751161, 12.886394,
- 13.024498, 13.165570, 13.309711, 13.457026, 13.607625,
- 13.761625, 13.919145, 14.080314, 14.245263, 14.414134,
- 14.587072, 14.764233, 14.945778, 15.131877, 15.322712,
- 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
- 16.578520, 16.808433, 17.044929, 17.288305, 17.538873,
- 17.796967, 18.062943, 18.337176, 18.620068, 18.912049,
- 19.213574, 19.525133, 19.847249, 20.180480, 20.525429,
- 20.882738, 21.253102, 21.637266, 22.036036, 22.450278,
- 22.880933, 23.329017, 23.795634, 24.281981, 24.789364,
- 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
- 28.365274, 29.068370, 29.808638, 30.589157, 31.413354,
- 32.285060, 33.208568, 34.188705, 35.230920, 36.341388,
- 37.527131, 38.796172, 40.157721, 41.622399, 43.202525,
- 44.912465, 46.769077, 48.792279, 51.005773, 53.437996,
- 56.123356, 59.103894
- };
- gnm_float X, U, V, RANLAN;
- int I;
-
- do {
- X = random_01 ();
- } while (X == 0.0);
- U = 1000.0 * X;
- I = U;
- U = U - I;
-
- if (I >= 70 && I <= 800)
- RANLAN = F[I] + U * (F[I + 1] - F[I]);
- else if (I >= 7 && I <= 980)
- RANLAN = F[I] + U * (F[I + 1] - F[I] -
- 0.25 * (1 - U) * (F[I + 2] - F[I + 1] -
- F[I] + F[I - 1]));
- else if (I < 7) {
- V = gnm_log (X);
- U = 1 / V;
- RANLAN = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * U) * U) /
- (1 + (3.41760202E1 + 4.01244582 * U) * U)) *
- ( -gnm_log ( -0.91893853 - V) - 1);
- } else {
- U = 1 - X;
- V = U * U;
- if (X <= 0.999)
- RANLAN = (1.00060006 + 2.63991156E2 *
- U + 4.37320068E3 * V) /
- ((1 + 2.57368075E2 * U + 3.41448018E3 * V) * U);
- else
- RANLAN = (1.00001538 + 6.07514119E3 * U +
- 7.34266409E5 * V) /
- ((1 + 6.06511919E3 * U + 6.94021044E5 * V) * U);
- }
-
- return RANLAN;
-}
-
-
-/* ------------------------------------------------------------------------ */
-
-/*
- * Generate a skew-normal distributed random number.
- *
- * based on the information provided at
- * http://azzalini.stat.unipd.it/SN/faq.html
- *
- */
-
-gnm_float
-random_skew_normal (gnm_float a)
-{
- gnm_float result;
- gnm_float delta = a / gnm_sqrt(1 + a * a);
- gnm_float u = random_normal ();
- gnm_float v = random_normal ();
-
- result = delta * u + gnm_sqrt (1-delta*delta) * v;
-
- return ((u < 0.) ? -result : result);
-}
-
-
-/* ------------------------------------------------------------------------ */
-
-/*
- * Generate a skew-t distributed random number.
- *
- * based on the information provided at
- * http://azzalini.stat.unipd.it/SN/faq.html
- *
- */
-
-gnm_float
-random_skew_tdist (gnm_float nu, gnm_float a)
-{
- gnm_float chi = random_chisq (nu);
- gnm_float z = random_skew_normal (a);;
-
- return (z * gnm_sqrt(nu/chi));
-}
-
-
-/* ------------------------------------------------------------------------ */
-
-
gnm_float
combin (gnm_float n, gnm_float k)
diff --git a/src/mathfunc.h b/src/mathfunc.h
index e73837e..d306a7c 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -21,6 +21,7 @@ G_BEGIN_DECLS
#define M_LN2gnum GNM_const(0.693147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996419)
#define M_LN10gnum GNM_const(2.302585092994045684017991454684364207601101488628772976033327900967572609677352480235997205089598)
#define M_SQRT2gnum GNM_const(1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327642)
+#define M_Egnum GNM_const(2.718281828459045235360287471352662497757247)
/* ------------------------------------------------------------------------- */
@@ -124,41 +125,6 @@ gnm_float qgeom (gnm_float p, gnm_float prob, gboolean lower_tail, gboolean log_
gnm_float dcauchy (gnm_float x, gnm_float location, gnm_float scale, gboolean give_log);
gnm_float pcauchy (gnm_float x, gnm_float location, gnm_float scale, gboolean lower_tail, gboolean log_p);
-/* Random number generation. */
-gnm_float random_01 (void);
-gnm_float random_poisson (gnm_float lambda);
-gnm_float random_binomial (gnm_float p, gnm_float trials);
-gnm_float random_negbinom (gnm_float p, gnm_float f);
-gnm_float random_exponential (gnm_float b);
-gnm_float random_bernoulli (gnm_float p);
-gnm_float random_normal (void);
-gnm_float random_cauchy (gnm_float a);
-gnm_float random_lognormal (gnm_float zeta, gnm_float sigma);
-gnm_float random_weibull (gnm_float a, gnm_float b);
-gnm_float random_laplace (gnm_float a);
-gnm_float random_rayleigh (gnm_float sigma);
-gnm_float random_rayleigh_tail (gnm_float a, gnm_float sigma);
-gnm_float random_gamma (gnm_float a, gnm_float b);
-gnm_float random_pareto (gnm_float a, gnm_float b);
-gnm_float random_fdist (gnm_float nu1, gnm_float nu2);
-gnm_float random_beta (gnm_float a, gnm_float b);
-gnm_float random_logistic (gnm_float a);
-gnm_float random_geometric (gnm_float p);
-gnm_float random_hypergeometric (gnm_float n1, gnm_float n2, gnm_float t);
-gnm_float random_logarithmic (gnm_float p);
-gnm_float random_chisq (gnm_float nu);
-gnm_float random_tdist (gnm_float nu);
-gnm_float random_gumbel1 (gnm_float a, gnm_float b);
-gnm_float random_gumbel2 (gnm_float a, gnm_float b);
-gnm_float random_levy (gnm_float c, gnm_float alpha);
-gnm_float random_levy_skew (gnm_float c, gnm_float alpha,
- gnm_float beta);
-gnm_float random_exppow (gnm_float a, gnm_float b);
-gnm_float random_landau (void);
-gnm_float random_gaussian_tail (gnm_float a, gnm_float sigma);
-gnm_float random_skew_normal (gnm_float a);
-gnm_float random_skew_tdist (gnm_float nu, gnm_float a);
-
/* The probability density functions. */
gnm_float random_exppow_pdf (gnm_float x, gnm_float a, gnm_float b);
gnm_float random_laplace_pdf (gnm_float x, gnm_float a);
diff --git a/src/tools/data-shuffling.c b/src/tools/data-shuffling.c
index 90bfd0e..aff6fd1 100644
--- a/src/tools/data-shuffling.c
+++ b/src/tools/data-shuffling.c
@@ -35,7 +35,7 @@
#include <command-context.h>
#include <goffice/goffice.h>
-#include "mathfunc.h"
+#include "gnm-random.h"
#include "data-shuffling.h"
#include "dao.h"
#include "expr.h"
diff --git a/src/tools/goal-seek.c b/src/tools/goal-seek.c
index 4314beb..c626ba4 100644
--- a/src/tools/goal-seek.c
+++ b/src/tools/goal-seek.c
@@ -15,8 +15,8 @@
#include "numbers.h"
#include "gnumeric.h"
#include "goal-seek.h"
+#include "gnm-random.h"
-#include "mathfunc.h"
#include <stdlib.h>
#include <math.h>
#include <limits.h>
diff --git a/src/tools/random-generator.c b/src/tools/random-generator.c
index 5ca7819..e1d462a 100644
--- a/src/tools/random-generator.c
+++ b/src/tools/random-generator.c
@@ -29,7 +29,7 @@
#include "gnumeric.h"
#include "random-generator.h"
-#include "mathfunc.h"
+#include "gnm-random.h"
#include "rangefunc.h"
#include "tools.h"
#include "value.h"
diff --git a/src/wbc-gtk-actions.c b/src/wbc-gtk-actions.c
index 8f97648..3c1f111 100644
--- a/src/wbc-gtk-actions.c
+++ b/src/wbc-gtk-actions.c
@@ -103,7 +103,7 @@ static GNM_ACTION_DEF (cb_file_save) { gui_file_save (wbcg, wb_control_view (WOR
static GNM_ACTION_DEF (cb_file_save_as) { gui_file_save_as (wbcg, wb_control_view (WORKBOOK_CONTROL (wbcg))); }
#ifndef HAVE_MKDTEMP
-#include <mathfunc.h>
+#include "gnm-random.h"
#ifdef G_OS_WIN32
#include <process.h>
#include <errno.h>
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]