random.cc

一个用MATLAB语言编写的摄像机标定工具箱,内容丰富

//
// random.cc
//  Original C source: D. P. Mitchell  92/10/02.
//  C++ class wrapping: Dan Wallach
//
// $Id: random.cc,v 1.1.1.1 2001/02/28 00:28:38 cstolte Exp $
//

#include <sgl/random.h>
#include <sgl/vecmath2.h>
#include <sgl/vecmath3.h>
#include <time.h>
#include <math.h>

#define M	(1<<31)
#define R	48
#define S	8
#define LOGR	6

/*
 *	physically-random initial state
 */
static int randomData[48] = {
    1947146662, 255872804, 611859725, 1698768494, 1168331706, 903266320, 156960988, 
    2108094856, 1962943837, 190733878, 1611490366, 1064800627, 750133665, 1388252873, 
    1292226713, 2024731946, 1741633152, 1048965367, 1242081949, 648944060, 2146570054, 
    1918590288, 1154878382, 1659238901, 1517088874, 1343505038, 702442230, 223319626, 
    2112415214, 1198445789, 1667134997, 81636584, 1477696870, 1396665379, 296474908, 
    221800638, 285530547, 506210295, 322166150, 242833116, 542659480, 1371231536, 
    164404762, 580757470, 1564914124, 2004579430, 389459662, 1039937051
};

Random::Random() {
    init((int)time((time_t *)0));    // initialize from the system clock
}

Random::Random(int seed) {
    init(seed);
}

Random::~Random()
{
}

void Random::init(int seed) {
    register int i;

    feed = 0;
    tap = R - S;
    borrow = 0;
    /*
     *	the seed pseudorandomly perturbes a
     *	physically random state.
     */
    seed = 1103515245*seed + 12345;
    seed = 1103515245*seed + 12345;
    seed = 1103515245*seed + 12345;
    for (i = 0; i < R; i++) {
	seed = 1103515245*seed + 12345;
	vec[i] = seed ^ randomData[i];
    }
    for (i = 0; i < R*LOGR; i++)
	(void) simpleInt31();
    for (i = 0; i < 571; i++)
	V[i] = simpleInt31();
    Y = simpleInt31();
    for (i = 0; i < R*LOGR; i++)
	(void) int31();
}

/*
 *      Marsaglia's subtract-with-borrow generator
 *      (The Annals of Applied Probability, 1991, 1(3), 462-480)
 *	  from code by J. Reeds 91/11/19.
 *
 *	Augmented by Bays & Durham's output shuffling.
 */
int Random::int31() {
    register int tmp;
    register int j;

    tmp = (vec[tap] - vec[feed]) - borrow;
    if (tmp < 0) {
	borrow = 1;
	tmp += M;
    } else
	borrow = 0;
    vec[feed] = tmp;
    if (++feed >= R)
	feed = 0;
    if (++tap >= R)
	tap = 0;
    j = Y & 511;			/* Bays & Durham */
    Y = V[j];
    V[j] = tmp;
    return Y;
}

int Random::simpleInt31() {
    register int tmp;

    tmp = (vec[tap] - vec[feed]) - borrow;
    if (tmp < 0) {
	borrow = 1;
	tmp += M;
    } else
	borrow = 0;
    vec[feed] = tmp;
    if (++feed >= R)
	feed = 0;
    if (++tap >= R)
	tap = 0;
    return tmp;			/* no Bays & Durham */
}

/*
 *	uniform random integer from 0 to n-1
 */
int Random::integer(int n) {
    int slop, v;

    slop = int(0x7fffffffL % n);
    do {
	v = int31();
    } while (v <= slop);
    return v % n;
}

/*
 *	uniform distribution in [0, 1.0)
 */
double Random::uniform() {
    double x;

    do {
	x = (double)int31() * (1.0/2147483648.0);
	x = (x + (double)int31()) * (1.0/2147483648.0);
    } while (x >= 1.0);
    return x;
}

/*
 *	uniform points in unit disk
 */
double Random::disk(double *xp, double *yp) {
    double x, y, s;

    do {
	x = 2.0*uniform() - 1.0;
	y = 2.0*uniform() - 1.0;
	s = x*x + y*y;
    } while (s >= 1.0);
    *xp = x;
    *yp = y;
    return s;
}

Vector2 Random::disk() {
    double x, y, s;

    do {
	x = 2.0*uniform() - 1.0;
	y = 2.0*uniform() - 1.0;
	s = x*x + y*y;
    } while (s >= 1.0);
    
    return Vector2(x,y);
}

/*
 *	normal distribution
 */
double
Random::normal(double sigma, double mean) {
    double v1, x;
    static double v2, s = 1.1;

    if (s > 1.0) {
	s = disk(&v1, &v2);
	x = v1*sqrt(-2.0*log(s)/s);
    } else {
	x = v2*sqrt(-2.0*log(s)/s);
	s = 1.1;
    }
    return mean + sigma*x;
}

void
Random::vector2(double *x, double *y) {
    double s, r1, r2;

    s = disk(&r1, &r2);
    *x = (r1*r1 - r2*r2)/s;
    *y = 2.0*r1*r2/s;
}

Vector2
Random::vector2() {
    double s, r1, r2;

    s = disk(&r1, &r2);
    return Vector2((r1*r1 - r2*r2)/s, 2.0*r1*r2/s);
}

/*
 *	random 3D unit vector
 */
void Random::vector3(double *x, double *y, double *z) {
    double w;

    *z = 2.0*uniform() - 1.0;
    vector2(x, y);
    w = sqrt(1.0 - *z * *z);
    *x *= w;
    *y *= w;
}

Vector3
Random::vector3() {
    double x, y, z, w;

    z = 2.0*uniform() - 1.0;
    vector2(&x, &y);
    w = sqrt(1.0 - z * z);
    return Vector3(x * w, y * w, z);
}

double Random::exponential(double mu) {
    return -mu * log(uniform());
}

double Random::gamma(double a) {
    double u, x, y, v, v1, v2, p, q, t;
    int k;

    k = (int)floor(a);
    if (a < 10.0 && a == (double)k) {
	for (u = 1.0; k > 0; --k)
	    u *= uniform();
	return - log(u);
    }
    if (a > 1.0) {
	do {
	    do {
		(void) disk(&v1, &v2);
		y = v2/v1;
		x = sqrt(2.0 * a - 1.0) * y + a - 1.0;
	    } while (x <= 0.0);
	    t = (1.0 + y*y) * exp((a - 1.0)*log(x/(a - 1.0))
				  - x + a - 1.0);
	} while (uniform() > t);
    } else {
	p = M_E/(a + M_E);
	do {
	    u = uniform();
	    v = uniform();
	    if (u < p) {
		x = pow(v, 1.0/a);
		q = exp(-x);
	    } else {
		x = 1 - log(v);
		q = pow(x, a - 1.0);
	    }
	} while (uniform() >= q);
    }
    return x;
}

double Random::beta(double a, double b) {
    double x1, x2;

    x1 = gamma(a);
    x2 = gamma(b);
    return x1 / (x1 + x2);
}

double Random::chi(double nu) {
    return 2.0*gamma(0.5*nu);
}

double Random::Fdist(double nu1, double nu2) {
    return chi(nu1)*nu2 / (chi(nu2)*nu1);
}

double Random::tdist(double nu) {
    return normal(1.0, 0.0)/sqrt(chi(nu)/nu);
}

int Random::geometric(double p) {
    return (int)ceil(log(uniform()) / log(1.0 - p));
}

int Random::binomial(double p, int trials) {
    register int i, n;
    int a, b;
    double x;

    if (trials < 10) {
	n = 0;
	for (i = 0; i < trials; i++)
	    n += uniform() < p;
	return n;
    } else {
	a = (int)(1.0 + floor(0.5*(double)trials));
	b = trials - a + 1;
	x = beta(a, b);
	if (x >= p)
	    return binomial(p/x, a - 1);
	else
	    return a + binomial((p - x)/(1.0 - x), b - 1);
    }
}

int Random::poisson(double mu) {
    register int n;
    double m;
    double x, u;

    if (mu < 10) {
	u = exp(-mu);
	x = 1.0;
	n = 0;
	do {
	    n++;
	    x *= uniform();
	} while (x > u);
	return n - 1;
    }
    m = floor(7.0*mu/8.0);
    x = gamma(m);
    if (x < mu)
	return (int)m + poisson(mu - x);
    else
	return binomial(mu/x, (int)(m - 1.0));
}