rand.c

随机数生成器

/***************************************************************************
* Rand.C  (Random Number Generators)                                       *
*                                                                          *
* Source file for pseudo-random number utilities.  Rand is the             *
* base class for several different algorithms for generating pseudo-random *
* numbers.  Any method can generate individual samples or arrays of        *
* samples using "Eval".  The random seed can be reset at any time by       *
* calling "Seed" with any integer.  Random permutations of the integers    *
* 0,1,...(n-1) are generated by "Perm(n,P)".                               *
*                                                                          *
*   HISTORY                                                                *
*      Name    Date        Description                                     *
*                                                                          *
*      arvo    08/04/97    Changed to virtual functions.                   *
*      arvo    06/06/93    Optimization, especially for array evaluators.  *
*      arvo    10/06/91    Converted to C++                                *
*      arvo    11/20/89    Added "gen_seed" function to handle.            *
*      arvo    10/30/89    "state" allocation now done in rand_alloc.      *
*      arvo    07/08/89    Initial coding.                                 *
*                                                                          *
*--------------------------------------------------------------------------*
* Copyright (C) 1989, James Arvo                                           *
*                                                                          *
* This program is free software; you can redistribute it and/or modify it  *
* under the terms of the GNU General Public License as published by the    *
* Free Software Foundation.  See http://www.fsf.org/copyleft/gpl.html      *
*                                                                          *
* This program is distributed in the hope that it will be useful, but      *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for    *
* any particular purpose.  See the GNU General Public License for more     *
* details.                                                                 *
*                                                                          *
***************************************************************************/
#include <stdio.h>
#include <math.h>
#include <util/Rand.h>

#ifndef ABS
#define ABS( x ) ((x) > 0 ? (x) : -(x))
#endif

/*-------------------------------------------------------------------------*
 * M E T H O D 1                                                           *
 *                                                                         *
 * From "Numerical Recipes," by William H. Press, Brian P. Flannery,       *
 * Saul A. Teukolsky, and William T. Vetterling, p. 197.                   *
 *                                                                         *
 *-------------------------------------------------------------------------*/
static const long   M1 = 714025;
static const long   IA =   1366;
static const long   IC = 150889;
static const double RM = 1.400512E-6;

float RandGen_1::Eval()
    {
    register long  *elem;
    register long  offset;
    register float rand;
    offset = 1 + ( 97 * index ) / M1;
    if( offset > 97 ) offset = 97;
    if( offset <  1 ) offset =  1;
    elem   = shuffle + offset;
    rand   = ( index = *elem ) * RM;
    *elem  = ( seed  = ( IA * seed + IC ) % M1 );
    return rand;
    }

void RandGen_1::Eval( int n, float *array )
    {
    register long *shfl = shuffle;
    register long *elem;
    register long offset;
    for( int i = 0; i < n; i++ ) 
        {
        offset   = 1 + ( 97 * index ) / M1;
        if( offset > 97 ) offset = 97;
        if( offset <  1 ) offset =  1;
        elem     = shfl + offset;
        *array++ = ( index = *elem ) * RM;
        *elem    = ( seed  = ( IA * seed + IC ) % M1 );
        }
    }

void RandGen_1::Seed( long seed )
    {
    long t = ( IC + ABS( seed ) + 1 ) % M1;
    for( register int k = 1; k <= 97; k++ )
        {
        t = ( IA * t + IC ) % M1;
        shuffle[k] = ABS( t );
        }
    t = ( IA * t + IC ) % M1;
    seed  = ABS( t );
    index = ABS( t );
    }

/*-------------------------------------------------------------------------*
 * M E T H O D 2                                                           *
 *                                                                         *
 * From "The Multiple Prime Random Number Generator," by Alexander Haas,   *
 * ACM Transactions on Mathematical Software, Vol. 13, No. 4, December     *
 * 1987, pp. 368-381.                                                      *
 *                                                                         *
 *-------------------------------------------------------------------------*/
float RandGen_2::Eval()
    {
    if( (m += 7    ) >=   9973 ) m -=  9871;
    if( (i += 1907 ) >=  99991 ) i -= 89989;
    if( (j += 73939) >= 224729 ) j -= 96233;
    r = ((r * m + i + j) % 100000) / 10;
    return r * 1.00010001E-4;
    }

void RandGen_2::Eval( int n, float *array )
    {
    for( register int k = 0; k < n; k++ ) 
        {
        if( (m += 7    ) >=   9973 ) m -=  9871;
        if( (i += 1907 ) >=  99991 ) i -= 89989;
        if( (j += 73939) >= 224729 ) j -= 96233;
        r = ((r * m + i + j) % 100000) / 10;
        *array++ = r * 1.00010001E-4;
        }
    }

void RandGen_2::Seed( long seed )
    {
    r = ABS( seed      );
    m = ABS( seed *  7 );
    i = ABS( seed * 11 );
    j = ABS( seed * 13 );
    if( m < 100    ) m += 100;
    if( i < 10000  ) i += 10000;
    if( j < 128000 ) j += 128000;
    }

/*-------------------------------------------------------------------------*
 * M E T H O D 3                                                           *
 *                                                                         *
 * From "A More Portable Fortran Random Number Generator," by Linus        *
 * Schrage, ACM Transactions on Mathematical Software, Vol. 5, No, 2,      *
 * June 1979, pp. 132-138.                                                 *
 *                                                                         *
 *-------------------------------------------------------------------------*/
static const long A3 = 16807;
static const long P3 = 2147483647;

float RandGen_3::Eval()
    {
    long xhi    = ix >> 16;
    long xalo   = ( ix & 0xFFFF ) * A3;
    long leftlo = xalo >> 16;
    long fhi    = xhi * A3 + leftlo;
    long k      = fhi >> 15;
    ix          = ( ((xalo - (leftlo << 16)) - P3) +
                    ((fhi - (k << 15)) << 16) ) + k;
    if( ix < 0 ) ix += P3;
    return ix * 4.656612875E-10;
    }

void RandGen_3::Eval( int n, float *array )
    {
    register long xhi, xalo, leftlo;
    register long fhi, k;
    for( register int i = 0; i < n; i++ ) 
        {
        xhi    = ix >> 16;
        xalo   = ( ix & 0xFFFF ) * A3;
        leftlo = xalo >> 16;
        fhi    = xhi * A3 + leftlo;
        k      = fhi >> 15;
        ix     = ( ((xalo - (leftlo << 16)) - P3) +
                   ((fhi - (k << 15)) << 16) ) + k;
        if( ix < 0 ) ix += P3;
        *array++ = ix * 4.656612875E-10;
        }
    }

void RandGen_3::Seed( long seed )
    {
    ix = ABS( seed );
    }

/*-------------------------------------------------------------------------*
 * R A N D : : P E R M        (Permutation)                                *
 *                                                                         *
 * This routine fills an integer array of length "len" with a random       *
 * permutation of the integers 0, 1, 2, ... (len-1).                       *
 *                                                                         *
 * For efficiency, the random numbers are generated in batches of up to    *
 * "Nmax" at a time.  The constant Nmax can be set to any value >= 1.      *
 *                                                                         *
 *-------------------------------------------------------------------------*/
static const int Nmax = 20;

void RandGen::Perm( int len, int perm[] )
    {
    float R[ Nmax ];    // A buffer for getting random numbers.
    int   L = len - 1;  // Total number of random numbers needed.
    int   N = 0;        // How many to generate when we call Eval.
    int   n = 0;        // The array index into R.

    // First initialize the array "perm" to the identity permutation.

    for( int j = 0; j < len; j++ ) perm[j] = j;

    // Now swap a random element in the front with the i'th element.
    // When i gets down to 0, we're done.

    for( int i = len - 1; i > 0; i-- )   // Element i is a swap candidate.
        {
        if( n == N )                     // Generate more random numbers.
            {
            N = ( L < Nmax ) ? L : Nmax; // Can't get more than "Nmax".
            Eval( N, R );                // Generate N random numbers.
            L -= N;                      // Decrement total counter.
            n  = 0;                      // Start index at beginning of R.
            }
        float r = ( i + 1 ) * R[ n++ ];  // Pick a float in [0,i+1].
        int   k = (int)r;                // Truncate r to an integer.
        if( k < i )                      // Disregard k == i and k == i+1.
            {
            int tmp = perm[i];           // Swap elements i and k.
            perm[i] = perm[k];
            perm[k] = tmp;
            }
        }
    }