kolibrios/programs/develop/tinypy/modules/random/random.c

#include <stdlib.h>
#include <time.h>
#include <errno.h>

#include <math.h>
#ifndef M_E
 #define M_E     2.7182818284590452354
#endif
#ifndef M_PI
 #define M_PI    3.14159265358979323846
#endif

/*
 * following random generators are mutual exclusive.
 */
#define __USE_MERSENNE_TWISTER_RANDOM_GENERATOR
/*#define __USE_WICHMANN_HILL_RANDOM_GENERATOR*/
/*#define __USE_POSIX_RANDOM_GENERATOR*/


#if defined(__USE_MERSENNE_TWISTER_RANDOM_GENERATOR)

/***********************************************************
 * following code is borrowed from Python's _randommodule.c
 ***********************************************************/
/* Random objects */

/* ------------------------------------------------------------------
   The code in this module was based on a download from:
      http://www.math.keio.ac.jp/~matumoto/MT2002/emt19937ar.html

   It was modified in 2002 by Raymond Hettinger as follows:

    * the principal computational lines untouched except for tabbing.

    * renamed genrand_res53() to random_random() and wrapped
      in python calling/return code.

    * genrand_int32() and the helper functions, init_genrand()
      and init_by_array(), were declared static, wrapped in
      Python calling/return code.  also, their global data
      references were replaced with structure references.

    * unused functions from the original were deleted.
      new, original C python code was added to implement the
      Random() interface.

   The following are the verbatim comments from the original code:

   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
*/

/* Period parameters -- These are all magic.  Don't change. */
#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 */

typedef struct {
    unsigned long state[N];
    int index;
    int has_seed;
} RandomObject;

/* generates a random number on [0,0xffffffff]-interval */
static unsigned long
genrand_int32(RandomObject *self)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */
    unsigned long *mt;

    mt = self->state;
    if (self->index >= N) { /* generate N words at one time */
        int kk;

        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];

        self->index = 0;
    }

    y = mt[self->index++];
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);
    return y;
}

/* initializes mt[N] with a seed */
static void
init_genrand(RandomObject *self, unsigned long s)
{
    int mti;
    unsigned long *mt;

    mt = self->state;
    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 */
    }
    self->index = mti;
    return;
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
static void
init_by_array(RandomObject *self, unsigned long init_key[], unsigned long key_length)
{
    unsigned int i, j, k;   /* was signed in the original code. RDH 12/16/2002 */
    unsigned long *mt;

    mt = self->state;
    init_genrand(self, 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 */
    return;
}

/*----------------------end of Mersenne Twister Algorithm----------------*/

/*************************************************************************
 * following are tinypy related stuffs.
 *************************************************************************/

static RandomObject _gRandom;       /* random object */

static tp_obj random_seed(TP)
{
    tp_obj arg = TP_DEFAULT(tp_None);

    if (arg.type == TP_NONE) {
        time_t now;

        (void)time(&now);
        init_genrand(&_gRandom, (unsigned long)now);
        _gRandom.has_seed = 1;
    } else if (arg.type == TP_NUMBER) {
        init_genrand(&_gRandom, (unsigned long)arg.number.val);
        _gRandom.has_seed = 1;
    } else if (arg.type == TP_STRING) {
        unsigned long seed;
        
        seed = (unsigned long)tp_hash(tp, arg);
        init_genrand(&_gRandom, seed);
        _gRandom.has_seed = 1;
    } else {
        tp_raise(tp_None,tp_printf(tp, "%s", "invalid argument for seed()"));
    }
    
    return (tp_None);
}

static tp_obj random_getstate(TP)
{
    tp_obj state_list = tp_list(tp);
    int i;

    for (i = 0; i < N; i++) {
        _tp_list_append(tp, state_list.list.val, tp_number(_gRandom.state[i]));
    }
    _tp_list_append(tp, state_list.list.val, tp_number(_gRandom.index));

    return (state_list);
}

/*
 * @state_list  must contain exactly N+1(625) integer.
 */
static tp_obj random_setstate(TP)
{
    tp_obj state_list = TP_OBJ();
    tp_obj state_elem;
    tp_obj len;
    int i;

    len = tp_len(tp, state_list);
    if (len.number.val != N+1) {
        tp_raise(tp_None,tp_printf(tp, "%s: state vector's size invalid(should be %d)", 
                __func__, N+1));
    }

    for (i = 0; i < N; i++) {
        state_elem = tp_get(tp, state_list, tp_number(i));
        _gRandom.state[i] = (unsigned long)state_elem.number.val;
    }
    state_elem = tp_get(tp, state_list, tp_number(i));
    _gRandom.index = (int)state_elem.number.val;

    return (tp_None);
}

/*
 * Jumpahead should be a fast way advance the generator n-steps ahead, but
 * lacking a formula for that, the next best is to use n and the existing
 * state to create a new state far away from the original.
 *
 * The generator uses constant spaced additive feedback, so shuffling the
 * state elements ought to produce a state which would not be encountered
 * (in the near term) by calls to random().  Shuffling is normally
 * implemented by swapping the ith element with another element ranging
 * from 0 to i inclusive.  That allows the element to have the possibility
 * of not being moved.  Since the goal is to produce a new, different
 * state, the swap element is ranged from 0 to i-1 inclusive.  This assures
 * that each element gets moved at least once.
 *
 * To make sure that consecutive calls to jumpahead(n) produce different
 * states (even in the rare case of involutory shuffles), i+1 is added to
 * each element at position i.  Successive calls are then guaranteed to
 * have changing (growing) values as well as shuffled positions.
 *
 * Finally, the self->index value is set to N so that the generator itself
 * kicks in on the next call to random().  This assures that all results
 * have been through the generator and do not just reflect alterations to
 * the underlying state.
 */
static tp_obj random_jumpahead(TP)
{
    long n = (long)TP_NUM();
    long i, j;
    unsigned long *mt;
    unsigned long tmp;

    mt = _gRandom.state;

    for (i = N-1; i > 1; i--) {
        j = n % i;
        if (j == -1L) {
            tp_raise(tp_None,tp_printf(tp, "error: %s: j = %ld", __func__, j));
        }
        tmp   = mt[i];
        mt[i] = mt[j];
        mt[j] = tmp;
    }

    for (i = 0; i < N; i++)
        mt[i] += i + 1;

    _gRandom.index = N;

    return (tp_None);
}

/* random_random is the function named genrand_res53 in the original code;
 * generates a random number on [0,1) with 53-bit resolution; note that
 * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
 * multiply-by-reciprocal in the (likely vain) hope that the compiler will
 * optimize the division away at compile-time.  67108864 is 2**26.  In
 * effect, a contains 27 random bits shifted left 26, and b fills in the
 * lower 26 bits of the 53-bit numerator.
 * The orginal code credited Isaku Wada for this algorithm, 2002/01/09.
 */
static tp_obj random_random(TP)
{
    RandomObject *self = &_gRandom;
    unsigned long a, b;

    if (! self->has_seed)
        random_seed(tp);

    a = genrand_int32(self)>>5;
    b = genrand_int32(self)>>6;
    
    return tp_number((a*67108864.0+b)*(1.0/9007199254740992.0));
}


#elif defined(__USE_WICHMANN_HILL_RANDOM_GENERATOR)

/************************************************* 
 * divmod(x, y)
 * a helper function borrowed from Python
 *************************************************/
/* Compute Python divmod(x, y), returning the quotient and storing the
 * remainder into *r.  The quotient is the floor of x/y, and that's
 * the real point of this.  C will probably truncate instead (C99
 * requires truncation; C89 left it implementation-defined).
 * Simplification:  we *require* that y > 0 here.  That's appropriate
 * for all the uses made of it.  This simplifies the code and makes
 * the overflow case impossible (divmod(LONG_MIN, -1) is the only
 * overflow case).
 */
#include <assert.h>

static int
divmod(int x, int y, int *r)
{
    int quo;

    assert(y > 0);
    quo = x / y;
    *r = x - quo * y;
    if (*r < 0) {
        --quo;
        *r += y;
    }
    assert(0 <= *r && *r < y);
    return quo;
}


typedef struct WH_RandomObject {
    struct seed {
        unsigned long   x;
        unsigned long   y;
        unsigned long   z;
    } seed;
    int has_seed;
} WH_RandomObject;

static WH_RandomObject _gWhRandom;

static tp_obj random_seed(TP)
{
    long a;
    int x, y, z;
    tp_obj arg = TP_DEFAULT(tp_None);

    if (arg.type == TP_NONE) {
        time_t now;

        (void)time(&now);
        a = (long)now * 256;
    } else if (arg.type == TP_NUMBER) {
        a = (long)arg.number.val;
    } else {
        tp_raise(tp_None,tp_printf(tp, "%s", "invalid argument for seed()"));
    }
    
    a = divmod(a, 30268, &x);
    a = divmod(a, 30306, &y);
    a = divmod(a, 30322, &z);
    _gWhRandom.seed.x = (int)x + 1;
    _gWhRandom.seed.y = (int)y + 1;
    _gWhRandom.seed.z = (int)z + 1;
    _gWhRandom.has_seed = 1;

    return (tp_None);
}

/*
 * following comments are from Python's random.py
 *---------------------------------------
 * Wichman-Hill random number generator.
 *
 * Wichmann, B. A. & Hill, I. D. (1982)
 * Algorithm AS 183:
 * An efficient and portable pseudo-random number generator
 * Applied Statistics 31 (1982) 188-190
 *
 * see also:
 *        Correction to Algorithm AS 183
 *        Applied Statistics 33 (1984) 123
 *
 *        McLeod, A. I. (1985)
 *        A remark on Algorithm AS 183
 *        Applied Statistics 34 (1985),198-200
 * This part is thread-unsafe:
 * BEGIN CRITICAL SECTION
 */
static tp_obj random_random(TP)
{
    long x, y, z;
    double r;

    if (! _gWhRandom.has_seed)
        random_seed(tp);

    x = _gWhRandom.seed.x;
    y = _gWhRandom.seed.y;
    z = _gWhRandom.seed.z;

    x = (171 * x) % 30269;
    y = (172 * y) % 30307;
    z = (170 * z) % 30323;

    _gWhRandom.seed.x = x;
    _gWhRandom.seed.y = y;
    _gWhRandom.seed.z = z;

    /*
     * Note:  on a platform using IEEE-754 double arithmetic, this can
     * never return 0.0 (asserted by Tim; proof too long for a comment).
     */
    errno = 0;
    r = fmod(((double)x/30269.0+(double)y/30307.0+(double)z/30323.0), 1.0);
    if (errno == EDOM)
        tp_raise(tp_None,tp_printf(tp, "%s", "fmod(): denominator can't be zero"));

    return tp_number(r);
}

/*
 * for compatibility
 */
static tp_obj random_setstate(TP)
{
    return (tp_None);
}

/*
 * for compatibility
 */
static tp_obj random_getstate(TP)
{
    return (tp_None);
}

/*
 * FIXME: risk of overflow.
 * following comments are from Python's random.py
 * --------------------------------------------
 * Act as if n calls to random() were made, but quickly.
 *
 * n is an int, greater than or equal to 0.
 *
 * Example use:  If you have 2 threads and know that each will
 * consume no more than a million random numbers, create two Random
 * objects r1 and r2, then do
 *      r2.setstate(r1.getstate())
 *      r2.jumpahead(1000000)
 * Then r1 and r2 will use guaranteed-disjoint segments of the full
 * period.
 */
static tp_obj random_jumpahead(TP)
{
    int n = (int)TP_NUM();
    long x, y, z;

    if (n < 0)
        tp_raise(tp_None,tp_printf(tp, "%s: n = %d invalid, should >= 0", __func__, n));

    x = _gWhRandom.seed.x;
    y = _gWhRandom.seed.y;
    z = _gWhRandom.seed.z;

    x = (int)(x * ((long)pow(171, n) % 30269)) % 30269;
    y = (int)(y * ((long)pow(172, n) % 30307)) % 30307;
    z = (int)(z * ((long)pow(170, n) % 30323)) % 30323;

    _gWhRandom.seed.x = x;
    _gWhRandom.seed.y = y;
    _gWhRandom.seed.z = z;

    return (tp_None);
}


#elif defined(__USE_POSIX_RANDOM_GENERATOR)

/*
 * judge whether seeded
 */
static int has_seed = 0;

static tp_obj random_seed(TP)
{
    tp_obj arg = TP_DEFAULT(tp_None);

    if (arg.type == TP_NONE) {
        time_t now;

        (void)time(&now);
        srandom((unsigned int)now);
        has_seed = 1;
    } else if (arg.type == TP_NUMBER) {
        srandom((unsigned long)arg.number.val);
        has_seed = 1;
    } else {
        tp_raise(tp_None,tp_printf(tp, "%s", "invalid argument for seed()"));
    }
    
    return (tp_None);
}

/*
 * random()
 *
 * generate successive pseudo random number ranging from [0.0, 1.0).
 * usually RAND_MAX is huge number, thus the periods of the success-
 * ive random number is very long, about 16*((2**31)-1).
 * NOTE: if seed() not called before random(), random() will
 * automatically call seed() with current time.
 */
tp_obj random_random(TP)
{
    double r = 0.0;

    if (! has_seed)
        random_seed(tp);

    r = (tp_num)random()/(tp_num)RAND_MAX;

    return (tp_number(r));
}

/*
 * setstate(state)
 *
 * for compatibility.
 */
tp_obj random_setstate(TP)
{
    return (tp_None);
}

/*
 * getstate()
 *
 * for compatibility.
 */
tp_obj random_getstate(TP)
{
    return (tp_None);
}

/*
 * jumpahead()
 *
 * for compatibility.
 */
tp_obj random_jumpahead(TP)
{
    return (tp_None);
}

#else

#error no underlying random generator is specified

#endif

/************************************************************
 * some usual distributions
 ************************************************************/

/*
 * return real number in range [a, b)
 * a and b can be negtive, but a must be less than b.
 */
tp_obj random_uniform(TP)
{
    double a = TP_NUM();
    double b = TP_NUM();
    double r = 0.0;
    tp_obj rvo;         /* random variable object */

    if (a >= b)
        tp_raise(tp_None,tp_printf(tp, "%s: a(%f) must be less than b(%f)", a, b));

    rvo = random_random(tp);
    r = a + (b - a) * rvo.number.val;
    
    return (tp_number(r));
}

/*
 * Normal distribution
 * @mu      mean
 * @sigma   standard deviation
 *-----------------------------
 * Uses Kinderman and Monahan method. Reference: Kinderman,
 * A.J. and Monahan, J.F., "Computer generation of random
 * variables using the ratio of uniform deviates", ACM Trans
 * Math Software, 3, (1977), pp257-260.
 */
tp_obj random_normalvariate(TP)
{
    double mu = TP_NUM();
    double sigma = TP_NUM();
    double NV_MAGICCONST;
    double u1, u2;
    double z, zz;
    double r = 0.0;
    tp_obj rvo;         /* random variable object */

    NV_MAGICCONST = 4.0 * exp(-0.5) / sqrt(2.0);
    while (1) {
        rvo = random_random(tp);
        u1  = rvo.number.val;
        rvo = random_random(tp);
        u2  = 1.0 - rvo.number.val;
        z   = NV_MAGICCONST * (u1 - 0.5) / u2;
        zz  = z * z / 4.0;
        if (zz <= - log(u2))
            break;
    }

    r = mu + z * sigma;

    return (tp_number(r));
}

/*
 * Log normal distribution
 *
 * If take natural logarithm on log normal distribution, normal
 * distribution with mean mu and standard deviation sigma will
 * return.
 * @mu      mean, can be any value
 * @sigma   standard deviation, must be > 0.
 */
tp_obj random_lognormvariate(TP)
{
    double mu = TP_NUM();
    double sigma = TP_NUM();
    tp_obj params;
    tp_obj normvar;     /* normal distribution variate */
    double r = 0.0;

    /*
     * call random_normalvariate() actually
     */
    params = tp_params_v(tp, 2, tp_number(mu), tp_number(sigma));
    normvar = tp_ez_call(tp, "random", "normalvariate", params);
    r = exp(normvar.number.val);

    return (tp_number(r));
}

/*
 * Exponential distribution
 *
 * @lambda      reciprocal of mean.
 * return value range (0, +inf)
 */
tp_obj random_expovariate(TP)
{
    double lambda = TP_NUM();
    double u, r;
    tp_obj rvo;

    do {
        rvo = random_random(tp);
        u = rvo.number.val;
    } while (u <= 0.0000001);

    r = -log(u) / lambda;

    return (tp_number(r));
}

/*
 * Circular data distribution.
 *
 * mu is the mean angle, expressed in radians between 0 and 2*pi, and
 * kappa is the concentration parameter, which must be greater than or
 * equal to zero.  If kappa is equal to zero, this distribution reduces
 * to a uniform random angle over the range 0 to 2*pi.
 *
 * mu:    mean angle (in radians between 0 and 2*pi)
 * kappa: concentration parameter kappa (>= 0)
 * if kappa = 0 generate uniform random angle
 *
 * Based upon an algorithm published in: Fisher, N.I.,
 * "Statistical Analysis of Circular Data", Cambridge
 * University Press, 1993.
 *
 * Thanks to Magnus Kessler for a correction to the
 * implementation of step 4.
 */
tp_obj random_vonmisesvariate(TP)
{
    double mu = TP_NUM();
    double kappa = TP_NUM();
    tp_obj rvo;
    double a, b, c, r;
    double u1, u2, u3, z, f;
    double theta;
    double TWOPI = 2.0 * M_PI;

    if (kappa <= 1e-6) {
        rvo = random_random(tp);
        theta = TWOPI * rvo.number.val;
        return (tp_number(theta));
    }

    a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa);
    b = (a - sqrt(2.0 * a))/(2.0 * kappa);
    r = (1.0 + b * b)/(2.0 * b);

    while (1) {
        rvo = random_random(tp);
        u1 = rvo.number.val;

        z = cos(M_PI * u1);
        f = (1.0 + r * z)/(r + z);
        c = kappa * (r - f);

        rvo = random_random(tp);
        u2 = rvo.number.val;

        if ((u2 < (c * (2.0 - c))) ||
            (u2 <= (c * exp(1.0 - c))))
            break;
    }

    rvo = random_random(tp);
    u3 = rvo.number.val;
    if (u3 > 0.5)
        theta = fmod(mu, TWOPI) + acos(f);
    else
        theta = fmod(mu, TWOPI) - acos(f);

    return (tp_number(theta));
}

/*
 * Gamma distribution.  Not the gamma function!
 *
 * Conditions on the parameters are alpha > 0 and beta > 0.
 */
tp_obj random_gammavariate(TP)
{
    double alpha = TP_NUM();
    double beta  = TP_NUM();
    tp_obj rvo;
    double res;
    double LOG4 = log(4.0);
    double SG_MAGICCONST = 1.0 + log(4.5);

    /*
     * alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
     * Warning: a few older sources define the gamma distribution in terms
     * of alpha > -1.0
     */
    if ((alpha <= 0.0) || (beta <= 0.0))
        tp_raise(tp_None,tp_printf(tp, "%s: alpha(%f) and beta(%f) must be > 0.0",
                __func__, alpha, beta));

    if (alpha > 1.0) {

        /*
         * Uses R.C.H. Cheng, "The generation of Gamma
         * variables with non-integral shape parameters",
         * Applied Statistics, (1977), 26, No. 1, p71-74
         */

        double ainv;
        double bbb, ccc;
        double u1, u2;
        double v, x, z, r;

        ainv = sqrt(2.0 * alpha - 1.0);
        bbb = alpha - LOG4;
        ccc = alpha + ainv;

        while (1) {
            rvo = random_random(tp);
            u1 = rvo.number.val;
            if (! ((1e-7 < u1) && (u1 < 0.9999999)))
                continue;
            rvo = random_random(tp);
            u2 = 1.0 - rvo.number.val;
            v = log(u1 / (1.0 - u1)) / ainv;
            x = alpha * exp(v);
            z = u1 * u1 * u2;
            r = bbb + ccc * v - x;
            if ((r + SG_MAGICCONST - 4.5 * z >= 0.0) ||
                (r >= log(z))) {
                res = x * beta;
                return (tp_number(res));
            }
        }
    }
    else if (alpha == 1.0) {

        /*
         * expovariate(1)
         */

        double u;

        do {
            rvo = random_random(tp);
            u = rvo.number.val;
        } while (u <= 1e-7);

        res = - log(u) * beta;
        return (tp_number(res));
    } else {

        /*
         * alpha is between 0 and 1 (exclusive)
         *
         * Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
         */

        double b, p, u, u1, x;

        while (1) {
            rvo = random_random(tp);
            u = rvo.number.val;
            b = (M_E + alpha) / M_E;
            p = b * u;
            if (p <= 1.0)
                /*FIXME: x = p ** (1.0/alpha)*/
                x = pow(p, 1.0/alpha);
            else
                x = - log((b - p) / alpha);
            rvo = random_random(tp);
            u1 = rvo.number.val;
            if (p > 1.0) {
                /*FIXME: if u1 <= x ** (alpha - 1.0):*/
                if (u1 <= pow(x, alpha - 1.0))
                    break;
            }
            else if (u1 <= exp(-x))
                break;
        }

        res = x * beta;
        return (tp_number(res));
    }
}

/*
 * Beta distribution.
 *
 *        Conditions on the parameters are alpha > 0 and beta > 0.
 *        Returned values range between 0 and 1.
 *        
 *        # This version due to Janne Sinkkonen, and matches all the std
 *        # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
 *        
 * See also:
 * http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
 * for Ivan Frohne's insightful analysis of why the original implementation:
 *
 *    def betavariate(self, alpha, beta):
 *        # Discrete Event Simulation in C, pp 87-88.
 *
 *        y = self.expovariate(alpha)
 *        z = self.expovariate(1.0/beta)
 *        return z/(y+z)
 *
 * was dead wrong, and how it probably got that way.
 *
 */
tp_obj random_betavariate(TP)
{
    double alpha = TP_NUM();
    double beta  = TP_NUM();
    double t;
    double r = 0.0;    
    tp_obj y;
    tp_obj params;

    params = tp_params_v(tp, 2, tp_number(alpha), tp_number(1.0));
    y = tp_ez_call(tp, "random", "gammavariate", params);
    if (y.number.val == 0) {
        return (y);
    } else {
        params = tp_params_v(tp, 2, tp_number(beta), tp_number(1.0));
        t = y.number.val;
        y = tp_ez_call(tp, "random", "gammavariate", params);
        r = t / (t + y.number.val);
        return (tp_number(r));
    }
}

/*
 * Pareto distribution.  alpha is the shape parameter.
 * # Jain, pg. 495
 */
tp_obj random_paretovariate(TP)
{
    double alpha = TP_NUM();
    double u;
    double r;
    tp_obj rvo;
    
    rvo = random_random(tp);
    u = 1.0 - rvo.number.val;
    r = 1.0 / pow(u, 1.0/alpha);
    
    return (tp_number(r));
}

/*
 * Weibull distribution.
 *
 * alpha is the scale parameter and beta is the shape parameter.
 *
 * Jain, pg. 499; bug fix courtesy Bill Arms
 */
tp_obj random_weibullvariate(TP)
{
    double alpha = TP_NUM();
    double beta  = TP_NUM();
    double u;
    double r;
    tp_obj rvo;
    
    rvo = random_random(tp);
    u = 1.0 - rvo.number.val;
    r = alpha * pow(-log(u), 1.0/beta);
    return (tp_number(r));
}

/*
 * randomly select an element from range ([start,] stop[, step])
 *
 * 'stop' must be larger than 'start', both can be negative;
 * 'step' must be integer larger than zero.
 */
tp_obj random_randrange(TP)
{
    tp_obj start = TP_OBJ();
    tp_obj stop = TP_DEFAULT(tp_None);
    tp_obj step = TP_DEFAULT(tp_number(1));
    tp_obj rvo = random_random(tp);
    int istart = (int)start.number.val;
    int istep = (int)step.number.val;
    int istop;
    int iwidth;
    double res;
    
    if (stop.type == TP_NONE) {
        /*
                        * if only one argument, then start just means stop
                        */
        istop = istart;
        res = (rvo.number.val * istop);
        return (tp_number(res));
    } else if (stop.type == TP_NUMBER) {
        istop = (int)stop.number.val;
        iwidth = istop - istart;
        if (iwidth < 0)
            tp_raise(tp_None,tp_printf(tp, "%s", "stop must be > start"));
        if (istep <= 0)
            tp_raise(tp_None,tp_printf(tp, "%s", "step must be integer larger than 0"));
            
        if (istep == 1) {
            res = (int)(istart + (int)(rvo.number.val * iwidth));
            return (tp_number(res));
        } else {
            int n = (iwidth + istep - 1) / istep;
            res = (int)(istart + istep * (int)(n * rvo.number.val));
            return (tp_number(res));
        }
    } else {
        tp_raise(tp_None,tp_printf(tp, "%s", "wrong type of stop"));
    }
}

/*
 * return random integer between [a, b]
 */
tp_obj random_randint(TP)
{
    double a = TP_NUM();
    double b = TP_NUM();
    tp_obj r;
    tp_obj params;

    params = tp_params_v(tp, 2, tp_number(a), tp_number(b + 1));
    r = tp_ez_call(tp, "random", "randrange", params);
    return (r);
}

/*
 * return a random element of sequence 'seq'. 'seq' mustn't be empty.
 */
tp_obj random_choice(TP)
{
    tp_obj seq = TP_OBJ();
    tp_obj len;
    tp_obj rvo;
    tp_obj r;
    int i;
    
    len = tp_len(tp, seq);
    if (len.number.val <= 0)
        tp_raise(tp_None,tp_printf(tp, "%s", "seq mustn't be empty"));
    
    rvo = random_random(tp);
    i = (int)(len.number.val * rvo.number.val);
    r = tp_get(tp, seq, tp_number(i));
    
    return (r);
}

/*
 * shuffle sequence 'seq' in place, return None
 */
tp_obj random_shuffle(TP)
{
    tp_obj seq = TP_OBJ();
    tp_obj elmi;
    tp_obj elmj;
    tp_obj params;
    tp_obj rvo;
    tp_obj len = tp_len(tp, seq);
    int i, j;
    
    if (len.number.val <= 0)
        return (tp_None);
    
    for (i = len.number.val - 1; i > len.number.val / 2; i--) {
        /*
                       * randomly exchange elment i and elment j, element i from the behind end of 'seq', while
                       * element j from the front end of 'seq'.
                       */
        params = tp_params_v(tp, 2, tp_number(0), tp_number(len.number.val / 2));
        rvo = tp_ez_call(tp, "random", "randint", params);
        j = (int)rvo.number.val;
        elmi = tp_get(tp, seq, tp_number(i));
        elmj = tp_get(tp, seq, tp_number(j));        
        
        tp_set(tp, seq, tp_number(i), elmj);
        tp_set(tp, seq, tp_number(j), elmi);
    }
    
    for (i = len.number.val / 2; i >= 0; i--) {
        /*
                       * randomly exchange elment i and elment j, element i from the front end of 'seq', while
                       * element j from the behind end of 'seq'.
                       */
        params = tp_params_v(tp, 2, tp_number(len.number.val / 2), tp_number(len.number.val - 1));
        rvo = tp_ez_call(tp, "random", "randint", params);
        j = (int)rvo.number.val;
        elmi = tp_get(tp, seq, tp_number(i));
        elmj = tp_get(tp, seq, tp_number(j));
        
        tp_set(tp, seq, tp_number(i), elmj);
        tp_set(tp, seq, tp_number(j), elmi);
    }
    
    return (tp_None);
}