📄 random.c
字号:
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Filename: random.c
//
// Purpose: Random utility procedures for BayeSys3.
//
// History: Random.c 17 Nov 1994 - 13 Sep 2003
//
// Acknowledgments:
// "Numerical Recipes", Press et al, for ideas
// "Handbook of Mathematical Functions", Abramowitz and Stegun, for formulas
// Peter J Acklam website, for inverse normal code
//-----------------------------------------------------------------------------
/*
Copyright (c) 1994-2003 Maximum Entropy Data Consultants Ltd,
114c Milton Road, Cambridge CB4 1XE, England
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#include "license.txt"
*/
#include <stddef.h>
#include <time.h>
#include <math.h>
#include <float.h>
#include "random.h"
static const double SHIFT32 = 1.0 / 4294967296.0; // 2^-32
static const double HalfPi = 1.57079632679489661922; // pi/2
static const double TwoxPi = 6.28318530717958647688; // 2*pi
static const double SqrPi2 = 1.25331413731550025120; // sqrt(pi/2)
static const double Sqr2Pi = 2.50662827463100050240; // sqrt(2*pi)
static const double LogSqr2Pi = 0.91893853320467274177; // log(sqrt(2*pi))
// Internal prototypes
static void Positive2 (Rand_t, double, double, double, double*, double*);
static void Positive2A (double, double, double, double,
double*, double*, double*);
static void Positive2B (double, double, double, double,
double*, double*, double*);
static void Positive2C (double, double, double, double,
double*, double*, double*);
static void Positive2D (double, double, double*, double*, double*);
static void Posneg2 (Rand_t, double, double, double, double, double,
double*, double*);
static void Posneg2A (double, double, double, double, double, double,
double*, double*, double*);
static void Posneg2B (double, double, double, double, double, double,
double*, double*, double*);
static void Posneg2C (double, double, double, double, double, double,
double*, double*, double*);
static void Posneg2D (double, double, double, double, double, double,
double*, double*, double*);
static double eigenerf2 (double, double, double, double, double);
static double erfint2 (double, double);
static double wedge (double, double, double);
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: RanInit
//
// Purpose: Initialise generator of 4 unsigneds for other Random procedures
//
// Format: Rand[0] + Rand[1] * 2^32 = # calls to Ranint since initialisation
// Rand[2] = generator offset, re-randomised every 2^32 calls
// Rand[3] = extra random integer available after each call
//
// History: John Skilling 15 Jan 2002, 31 Oct 2002
//-----------------------------------------------------------------------------
//
int RanInit( // O Seed, either from input or time
Rand_t Rand, // O Random generator state [4]
int seed) // I Seed: +ve = value, -ve = time seed
{
unsigned j, k;
k = 1;
for( j = 0; k; ++j )
k += k;
if( j != 32 ) // Check 32-bit arithmetic
return E_RAN_ARITH;
if( seed < 0 )
{
seed = (int)time(NULL);
if( seed < 0 )
seed = ~seed; // still OK after A.D.2030
}
Rand[0] = Rand[1] = 0; // 64-bit counter
Rand[2] = 1013904223 + 1664525 * seed; // sticky offset
Rand[3] = 1013904223 + 1664525 * Rand[2]; // extra random integer
return seed;
}
#if 0
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative initialiser #1, crude congruential.
//-----------------------------------------------------------------------------
void RanInit1(
Rand_t Rand, // O Random generator state [1]
int seed) // I Non-negative seed
{
Rand[0] = seed;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative initialiser #2, double congruential with shuffling.
//-----------------------------------------------------------------------------
void RanInit2(
Rand_t Rand, // O Random generator state [35]
int seed) // I Non-negative seed
{
int i, j;
i = (seed > 2147483397) ? seed - 2147483397 : seed + 1; // [1...2147483398]
Rand[33] = i;
for( j = 39; j >= 0; j-- )
{
i = 40014 * i - 2147483563 * (i / 53668);
if( i < 0 )
i += 2147483563; // [1...2147483562]
if( j < 32 )
Rand[j] = i;
}
Rand[34] = Rand[0];
Rand[32] = i;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative initialiser #3, Knuth's subtractive.
//-----------------------------------------------------------------------------
void RanInit3(
Rand_t Rand, // O Random generator state [112]
int seed) // I Non-negative seed
{
unsigned i, j, k, m;
m = (unsigned)(161803398 - seed);
Rand[54] = m;
Rand[109] = 0;
k = 1;
for( i = 1; i < 55; ++i )
{
j = (21 * i - 1) % 55;
Rand[j] = k;
k = m - k;
m = Rand[j];
Rand[54+i] = i;
}
for( k = 0; k < 4; ++k )
for( i = 0; i < 55; ++i )
Rand[i] -= Rand[(i + 31) % 55];
Rand[110] = 0;
Rand[111] = 31;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative initialiser #4, 64-bit hashing.
//-----------------------------------------------------------------------------
void RanInit4(
Rand_t Rand, // O Random generator state [3]
int seed) // I Non-negative seed
{
Rand[0] = Rand[1] = 0;
Rand[2] = seed;
}
#endif
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Ranint
//
// Purpose: Random integer sample, in [-2^31,2^31).
//
// Method: 64-bit hashing is on #calls + offset to produce 32-bit candidate
// with algorithm from ran4 in "Numerical Recipes".
// After each 2^32 calls, the 64-bit offset is incremented by
// congruential generators on its 32-bit halves.
//
// Notes: (1) Period is just less than 2^95, extremely large.
// (2) All random calls are directed through this procedure.
// (3) After a call, Rand[3] is available as another random integer.
//
// History: John Skilling 28 Jan 2002, 31 Oct 2002, 17 Dec 2002
//-----------------------------------------------------------------------------
//
int Ranint( // O 32-bit integer
Rand_t Rand) // I O Random generator state [4]
{
unsigned m, n; // 64-bit register
unsigned u, v, w;
int i;
// 64-bit counter, for hashing
if( !(Rand[0] ++) )
{
Rand[1] ++;
// occasional update offset in Rand[2]
i = Rand[2] >> 1; // [0...2147483647]
i -= i / 24683721; // [0...2147483561]
i++; // [1...2147483562]
i = 40014 * i - 2147483563 * (i / 53668);
if( i < 0 ) // (40014 * i) % 2147483563
i += 2147483563; // [1...2147483562]
i--; // [0...2147483561]
i += i / 24683720; // [0...2147483647] (with holes)
Rand[2] = i << 1; // even
if( 1013904223 + 1664525 * i < 0 ) // chance
Rand[2]++; // even or odd
}
// Two double steps of 64-bit hash
n = Rand[0] + Rand[2];
m = Rand[1] + Rand[2];
w = n ^ 0xbaa96887;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
m ^= (((u >> 16) | (u << 16)) ^ 0xb4f0c4a7) + w * v;
w = m ^ 0x1e17d32c;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
n ^= (((u >> 16) | (u << 16)) ^ 0x178b0f3c) + w * v;
w = n ^ 0x03bcdc3c;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
m ^= (((u >> 16) | (u << 16)) ^ 0x96aa3a59) + w * v;
w = m ^ 0x0f33d1b2;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
n ^= (((u >> 16) | (u << 16)) ^ 0xaa5835b9) + w * v;
Rand[3] = m;
return n;
}
#if 0
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative generator #1, crude congruential.
//-----------------------------------------------------------------------------
int Ranint1( // O 32-bit integer
Rand_t Rand) // I O Random generator state [1]
{
Rand[0] = Rand[0] * 1664525 + 1013904223;
return Rand[0];
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative generator #2, double congruential with shuffling.
//-----------------------------------------------------------------------------
int Ranint2( // O Odd integer with 85 holes: ONLY (2^31 - 85) VALUES
Rand_t Rand) // I O Random generator state [35]
{
int i, j, k;
// First congruential sequence
i = Rand[32];
i = 40014 * i - 2147483563 * (i / 53668);
if( i < 0 )
i += 2147483563;
Rand[32] = i; // [1...2147483562]
// Second congruential sequence
k = Rand[33];
k = 40692 * k - 2147483399 * (k / 52774);
if( k < 0 )
k += 2147483399;
Rand[33] = k; // [1...2147483398]
// Subtract from shuffled table
j = Rand[34] / 67108862; // [0...31]
k = Rand[j] - k; // [-2147483397...2147483561]
if( k < 0 )
k += 2147483563; // (2^31 - 85) values
Rand[34] = k; // [0...2147483562]
Rand[j] = i;
// Expand result k
k += k / 24970740; // [0...2^31) with holes at 24970740*[1...85]
return (k << 1) | 1;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative generator #3, Knuth's subtractive.
//-----------------------------------------------------------------------------
int Ranint3( // O 32-bit integer
Rand_t Rand) // I O Random generator state [112]
{
unsigned j, k;
j = Rand[110]; Rand[110] = Rand[55+j];
k = Rand[111]; Rand[111] = Rand[55+k];
return (Rand[j] -= Rand[k]);
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Kernel of alternative generator #4, 64-bit hashing.
//-----------------------------------------------------------------------------
int Ranint4( // O 32-bit integer
Rand_t Rand) // I O Random generator state [3]
{
unsigned u, v, w, m, n;
// 64-bit counter, for hashing
if( !(Rand[0] ++) )
Rand[1] ++;
// 64-bit hash
n = Rand[0] + Rand[2];
m = Rand[1] + Rand[2];
w = n ^ 0xbaa96887;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
m ^= (((u >> 16) | (u << 16)) ^ 0xb4f0c4a7) + w * v;
w = m ^ 0x1e17d32c;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
n ^= (((u >> 16) | (u << 16)) ^ 0x178b0f3c) + w * v;
w = n ^ 0x03bcdc3c;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
m ^= (((u >> 16) | (u << 16)) ^ 0x96aa3a59) + w * v;
w = m ^ 0x0f33d1b2;
v = w >> 16;
w &= 0xffff;
u = (v - w) * (v + w);
n ^= (((u >> 16) | (u << 16)) ^ 0xaa5835b9) + w * v;
return n;
}
#endif
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Ranfloat
//
// Purpose: Random single-precision floating-point sample.
// The 2^23 allowed values are odd multiples of 2^-24,
// symmetrically placed in strict interior of (0,1).
//
// Notes: (1) Tuned to 23-bit mantissa in "float" format.
// (2) Uses Ranint.
//
// History: John Skilling 20 Oct 2002
//-----------------------------------------------------------------------------
//
float Ranfloat( // O Value
Rand_t Rand) // I O Random generator state
{
unsigned u = ((unsigned)Ranint(Rand)
& 0xfffffe00) ^ 0x00000100; // switch lowest (2^-24) bit ON
return (float)u * (float)SHIFT32;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Randouble
//
// Purpose: Random double-precision floating-point sample.
// The 2^52 allowed values are odd multiples of 2^-53,
// symmetrically placed in strict interior of (0,1).
//
// Notes: (1) Tuned to 52-bit mantissa in "double" format.
// (2) Uses one call to Ranint to get 64 random bits, with extra
// random integer available in Rand[3].
// (3) All floating-point random calls are directed through this code,
// except Rangauss which uses the extra random integer in Rand[3].
//
// History: John Skilling 6 May 1995, 3 Dec 1995, 24 Aug 1996
// 20 Oct 2002, 17 Dec 2002
//-----------------------------------------------------------------------------
//
double Randouble( // O Value within (0,1)
Rand_t Rand) // I O Random generator state
{
unsigned hi, lo;
hi = (unsigned)Ranint(Rand); // top 32 bits
lo = (Rand[3] // low bits
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -