乱数生成プログラム

/*	random.c	擬似乱数(Pseudo-Random numbers)	*/
#include <stdio.h>
#include <stdlib.h>
#include "random.h"

/*#define		DEBUG		/* debug routine */

ulong _idum_u0, _idum_u1, _idum_u2;		/* urnd0()〜urnd6()の乱数核 */
ulong _idum_u3, _idum_u4, _idum_u5;		/*      〃        */
ulong _idum_u6;							/*      〃        */
ulong *_prx, *_prx_end;					/* mrnd()の乱数核ほか変数 */
ulong _rx[522];							/*      〃        */
unsigned int r_s1, r_s2, r_s3;			/* rrnd()の乱数核ほか変数 */
long _idum_o0, _idum_o1, _idum_o2;		/* ornd0()〜ornd2()の乱数核 */
long _idum_b, b_nd, b_iy;				/* brnd()の乱数核ほか変数 */
long b_iv[NTAB];						/*      〃        */
long _idum_e1, _idum_e2, e_nd, e_iy;	/* ernd()の乱数核ほか変数 */
long e_iv[NTAB];						/*      〃        */
long *k_ma_end, *pk1, *pk2;				/* krnd()の乱数核ほか変数 */
long k_ma[57];							/*      〃        */
ulong _idum_d;							/* drnd()の乱数核 */
ulong _idum_q;							/* qrnd()の乱数核 */
ulong _idum_p0, _idum_p1, _idum_p2;		/* prnd0()〜prnd33()の乱数核 */
ulong _idum_p3, _idum_p4, _idum_p5;		/*      〃        */
ulong _idum_p6, _idum_p7, _idum_p8;		/*      〃        */
ulong _idum_p9, _idum_p10, _idum_p11;	/*      〃        */
ulong _idum_p12, _idum_p13, _idum_p14;	/*      〃        */
ulong _idum_p15, _idum_p16, _idum_p17;	/*      〃        */
ulong _idum_p18, _idum_p19, _idum_p20;	/*      〃        */
ulong _idum_p21, _idum_p22, _idum_p23;	/*      〃        */
ulong _idum_p24, _idum_p25, _idum_p26;	/*      〃        */
ulong _idum_p27, _idum_p28, _idum_p29;	/*      〃        */
ulong _idum_p30, _idum_p31, _idum_p32;	/*      〃        */
ulong _idum_p33;						/*      〃        */
double (*_rand)(void);					/* 一様乱数以外の乱数で呼び出す一様分布乱数 */ 
double _tri;

void _error(char *s, int r)
{
	fprintf(stderr, "Error : %s\n", s);
	exit(r);
}

void init_urnd0(long seed)
{
	_idum_u0 = seed ^ MASK;
}

void init_urnd1(long seed)
{
	_idum_u1 = seed ^ MASK;
}

void init_urnd2(long seed)
{
	_idum_u2 = seed ^ MASK;
}

void init_urnd3(long seed)
{
	_idum_u3 = seed ^ MASK;
}

void init_urnd4(long seed)
{
	_idum_u4 = seed ^ MASK;
}

void init_urnd5(long seed)
{
	_idum_u5 = seed ^ MASK;
}

void init_urnd6(long seed)
{
	_idum_u6 = seed ^ MASK;
}

void init_mrnd(long seed)
{
	int i, j;
	ulong *p, *p1, *p2, *p3, s, u;

	_prx_end = &_rx[521];

	s = seed ^ MASK;
	u = 0;
	for(i = 0, _prx = _rx; i <= 16; i++)
	{
		for(j = 0; j < 32; j++)
		{
			s = s * 1566083941L + 1;
			u = (u >> 1) | (s & (1L << 31));
		}
		*_prx++ = u;
	}
	_prx--;
	*_prx = (0xffffffff & (*_prx << 23)) ^ (*_rx >> 9) ^ *(_rx + 15);
	p1 = _rx;
	p2 = &_rx[1];
	p3 = _prx++;
	while(_prx != _prx_end)
		*_prx++ = (0xffffffff & (*p1++ << 23)) ^ (*p2++ >> 9) ^ *p3++;
	rnd521();		/* wram up */
	rnd521();
	rnd521();
	_prx = &_rx[520];
}

void init_rrnd(long seed)
{
	r_s1 = (seed ^ MASK) & 0x7fff;
	r_s2 = ((r_s1 * r_s1) ^ MASK) & 0x7fff;
	r_s3 = ((r_s1 * r_s2) ^ MASK) & 0x7fff;
#ifdef DEBUG
printf("r_s1 = %d, r_s2 = %d, r_s3 = %d\n", r_s1, r_s2, r_s3);
#endif
}

void init_ornd0(long seed)
{
	_idum_o0 = seed ^ MASK;
}

void init_ornd1(long seed)
{
	_idum_o1 = seed ^ MASK;
}

void init_ornd2(long seed)
{
	_idum_o2 = seed ^ MASK;
}

long b_iy;

void init_brnd(long seed)
{
	int i;
	long *p;

	_idum_b = seed ^ MASK;
	for(i = 0; i < 8; i++)	bcal();
	for(i = NTAB - 1; i >= 0; i--)	b_iv[i] = bcal();
	b_nd = 1 + 2147483646L / NTAB;
	b_iy = b_iv[0];
} 

void init_ernd(long seed)
{
	int i;

	_idum_e1 = seed ^ MASK;
	_idum_e2 = (_idum_e1 * _idum_e1) ^ MASK;
	for(i = 0; i < 8; i++)	ecal1();
	for(i = NTAB - 1; i >= 0; i--)	e_iv[i] = ecal1();
	e_nd = 1 + 2147483562L / NTAB;
	e_iy = e_iv[0];
} 

void init_krnd(long seed)
{
	long mj, mk;
	int i, k;

	k_ma_end = &k_ma[56];
	mj = (MSEED - (seed ^ MASK)) % MBIG;
	k_ma[55] = mj;
	mk = 1L;
	for(i = 1; i <= 54; i++)
	{
		k = (21 * i) % 55;
		k_ma[k] = mk;
		mk = mj - mk;
		if(mk < 0)	mk += MBIG;
		mj = k_ma[k];
	}
	for(k = 1; k <= 4; k++)
	{
		for(i = 1; i <= 55; i++)
		{
			k_ma[i] -= k_ma[1 + (i + 30) % 55];
			if(k_ma[i] < 0)	k_ma[i] += MBIG;
		}
	}
	pk1 = &k_ma[0];
	pk2 = &k_ma[31];
}

void init_drnd(long seed)
{
	_idum_d = seed ^ MASK;
}

void init_qrnd(long seed)
{
	_idum_q = seed ^ MASK;
}

long init_prnd0(long seed)
{
	_idum_p0 = (seed ^ MASK) & 0x00000fff;
	return 6075;
}

long init_prnd1(long seed)
{
	_idum_p1 = (seed ^ MASK) & 0x00000fff;
	return 7875;
}

long init_prnd2(long seed)
{
	_idum_p2 = (seed ^ MASK) & 0x00000fff;
	return 7875;
}

long init_prnd3(long seed)
{
	_idum_p3 = (seed ^ MASK) & 0x00000fff;
	return 6075;
}

long init_prnd4(long seed)
{
	_idum_p4 = (seed ^ MASK) & 0x00000fff;
	return 6655;
}

long init_prnd5(long seed)
{
	_idum_p5 = (seed ^ MASK) & 0x00000fff;
	return 11979;
}

long init_prnd6(long seed)
{
	_idum_p6 = (seed ^ MASK) & 0x00000fff;
	return 14406;
}

long init_prnd7(long seed)
{
	_idum_p7 = (seed ^ MASK) & 0x00000fff;
	return 29282;
}

long init_prnd8(long seed)
{
	_idum_p8 = (seed ^ MASK) & 0x00000fff;
	return 53125;
}

long init_prnd9(long seed)
{
	_idum_p9 = (seed ^ MASK) & 0x00000fff;
	return 12960;
}

long init_prnd10(long seed)
{
	_idum_p10 = (seed ^ MASK) & 0x00000fff;
	return 14000;
}

long init_prnd11(long seed)
{
	_idum_p11 = (seed ^ MASK) & 0x00000fff;
	return 21870;
}

long init_prnd12(long seed)
{
	_idum_p12 = (seed ^ MASK) & 0x00000fff;
	return 31104;
}

long init_prnd13(long seed)
{
	_idum_p13 = (seed ^ MASK) & 0x00000fff;
	return 139968;
}

long init_prnd14(long seed)
{
	_idum_p14 = (seed ^ MASK) & 0x00000fff;
	return 29282;
}

long init_prnd15(long seed)
{
	_idum_p15 = (seed ^ MASK) & 0x00000fff;
	return 81000;
}

long init_prnd16(long seed)
{
	_idum_p16 = (seed ^ MASK) & 0x00000fff;
	return 134456;
}

long init_prnd17(long seed)
{
	_idum_p17 = (seed ^ MASK) & 0x00000fff;
	return 86436;
}

long init_prnd18(long seed)
{
	_idum_p18 = (seed ^ MASK) & 0x00000fff;
	return 121500;
}

long init_prnd19(long seed)
{
	_idum_p19 = (seed ^ MASK) & 0x00000fff;
	return 259200;
}

long init_prnd20(long seed)
{
	_idum_p20 = (seed ^ MASK) & 0x00000fff;
	return 117128;
}

long init_prnd21(long seed)
{
	_idum_p21 = (seed ^ MASK) & 0x00000fff;
	return 121500;
}

long init_prnd22(long seed)
{
	_idum_p22 = (seed ^ MASK) & 0x00000fff;
	return 312500;
}

long init_prnd23(long seed)
{
	_idum_p23 = (seed ^ MASK) & 0x00000fff;
	return 145800;
}

long init_prnd24(long seed)
{
	_idum_p24 = (seed ^ MASK) & 0x00000fff;
	return 175000;
}

long init_prnd25(long seed)
{
	_idum_p25 = (seed ^ MASK) & 0x00000fff;
	return 233280;
}

long init_prnd26(long seed)
{
	_idum_p26 = (seed ^ MASK) & 0x00000fff;
	return 244944;
}

long init_prnd27(long seed)
{
	_idum_p27 = (seed ^ MASK) & 0x00000fff;
	return 139968;
}

long init_prnd28(long seed)
{
	_idum_p28 = (seed ^ MASK) & 0x00000fff;
	return 214326;
}

long init_prnd29(long seed)
{
	_idum_p29 = (seed ^ MASK) & 0x00000fff;
	return 714025;
}

long init_prnd30(long seed)
{
	_idum_p30 = (seed ^ MASK) & 0x00000fff;
	return 134456;
}

long init_prnd31(long seed)
{
	_idum_p31 = (seed ^ MASK) & 0x00000fff;
	return 259200;
}

long init_prnd32(long seed)
{
	_idum_p32 = (seed ^ MASK) & 0x00000fff;
	return 233280;
}

long init_prnd33(long seed)
{
	_idum_p33 = (seed ^ MASK) & 0x00000fff;
	return 714025;
}

double urnd0(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u0 *= 69069L) * d;
}

double urnd1(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u1 *= 1664525L) * d;
}

double urnd2(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u2 *= 39894229L) * d;
}

double urnd3(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u3 *= 48828125L) * d;
}

double urnd4(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u4 *= 1566083941L) * d;
}

double urnd5(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u5 *= 1812433253L) * d;
}

double urnd6(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_u6 *= 2110005341L) * d;
}

void rnd521(void)
{
	ulong *p, *p1, *p2;

	p1 = (_prx = p = p2 = _rx) + 489;
	while(p1 != _prx_end)	*p++ ^= *p1++;
	while(p  != _prx_end)	*p++ ^= *p2++;
}

double mrnd(void)
{
	double d = 1. / 4294967296.;
	if(++_prx == _prx_end)	rnd521();
	return (double)*_prx * d;
}

double rrnd(void)
{
	double d1 = 1. / 30269., d2 = 1. / 30307., d3 = 1. / 30323.;
	double r;

	if((r = (double)(r_s1 = (171 * r_s1) % 30269) * d1
		  + (double)(r_s2 = (172 * r_s2) % 30307) * d2
		  + (double)(r_s3 = (170 * r_s3) % 30323) * d3) < 1.)	return r;
	if(r < 2.)	return r - 1.;
	return r - 2.;
}

double ornd0(void)					/* Park and Miller の「最低基準」乱数 */
{
	long w;
	double d = 1. / 2147483647.;

	w = _idum_o0 / 127773L;
	if((_idum_o0 = (_idum_o0 - w * 127773L) * 16807 - w * 2836) <= 0)
		_idum_o0 += 2147483647L;
	return (double)_idum_o0 * d;
}

double ornd1(void)					/* Park and Miller の「最低基準」乱数 */
{
	long w;
	double d = 1. / 2147483647.;

	w = _idum_o1 / 44488L;
	if((_idum_o1 = (_idum_o1 - w * 44488L) * 48271L - w * 3399) <= 0)
		_idum_o1 += 2147483647L;
	return (double)_idum_o1 * d;
}

double ornd2(void)					/* Park and Miller の「最低基準」乱数 */
{
	long w;
	double d = 1. / 2147483647.;

	w = _idum_o2 / 30845;
	if((_idum_o2 = (_idum_o2 - w * 30845) * 69621L - w * 23902) <= 0)
		_idum_o2 += 2147483647L;
	return (double)_idum_o2 * d;
}

double brnd(void)
{
	double d = 1. / 2147483647.;
	int j;

	j = b_iy / b_nd;
	b_iy = b_iv[j];
	return (double)(b_iv[j] = bcal()) * d;
}

long bcal(void)
{
	long w;

	w = _idum_b / 127773L;
	if((_idum_b = (_idum_b - w * 127773L) * 16807 - w * 2836) <= 0)
		_idum_b += 2147483647L;
	return _idum_b;
}

double ernd(void)
{
	double d = 1. / 2147483563.;
	int j;

	j = e_iy / e_nd;
	if((e_iy = e_iv[j] - ecal2()) < 1)	e_iy += 2147483562L;
	return (double)(e_iv[j] = ecal1()) * d;
}

long ecal1(void)
{
	long w;

	w = _idum_e1 / 53668L;
	if((_idum_e1 = (_idum_e1 - w * 53668L) * 40014L - w * 12211L) <= 0)
		_idum_e1 += 2147483563L;
	return _idum_e1;
}

long ecal2(void)
{
	long w;

	w = _idum_e2 / 52774L;
	if((_idum_e2 = (_idum_e2 - w * 52774L) * 40692L - w * 3791) <= 0)
		_idum_e2 += 2147483399L;
	return _idum_e2;
}

double krnd(void)
{
	double d = 1. / (double)MBIG;
	long mj;

	if(++pk1 == k_ma_end)	pk1 = &k_ma[1];
	if(++pk2 == k_ma_end)	pk2 = &k_ma[1];
	if((mj = *pk1 - *pk2) < 0)	mj += MBIG;
	return (double)(*pk1 = mj) * d;
}

double drnd(void)
{
	double d = 1. / 4294967296.;
	ulong lw, rw, a, b, th, tl;

	a = _idum_d ^ 0xbaa96887L;
	b = (tl = a & 0xffff) * tl + ~((th = a >> 16) * th);
	a = (lw = ((((b >> 16)|(b << 16)) ^ 0x4b0f3b58L) + tl * th)) ^ 0x1e17d32cL;
	b = (tl = a & 0xffff) * tl + ~((th = a >> 16) * th);
	a = (rw = _idum_d++ ^ ((((b >> 16)|(b << 16)) ^ 0xe874f0c3L) + tl * th))
		^ 0x03bcdc3cL;
	b = (tl = a & 0xffff) * tl + ~((th = a >> 16) * th);
	a = lw ^ ((((b >> 16)|(b << 16)) ^ 0x6955c5a6L) + tl * th) ^ 0x0f33d1b2L;
	b = (tl = a & 0xffff) * tl + ~((th = a >> 16) * th);
	return (double)(rw ^ ((((b >> 16)|(b << 16)) ^ 0x55a7ca46L) + tl * th)) * d;
}

double qrnd(void)
{
	double d = 1. / 4294967296.;
	return (double)(_idum_q = _idum_q * 1664525L + 1013904223L) * d;
}

double prnd0(void)
{
	static double d = 1. / 6075.;
	return (double)(_idum_p0 = (_idum_p0 * 106 + 1283) % 6075) * d;
}

double prnd1(void)
{
	static double d = 1. / 7875.;
	return (double)(_idum_p1 = (_idum_p1 * 211 + 1663) % 7875) * d;
}

double prnd2(void)
{
	static double d = 1. / 7875.;
	return (double)(_idum_p2 = (_idum_p2 * 421 + 1663) % 7875) * d;
}

double prnd3(void)
{
	static double d = 1. / 6075.;
	return (double)(_idum_p3 = (_idum_p3 * 1366 + 1283) % 6075) * d;
}

double prnd4(void)
{
	static double d = 1. / 6655.;
	return (double)(_idum_p4 = (_idum_p4 * 936 + 1399) % 6655) * d;
}

double prnd5(void)
{
	static double d = 1. / 11979.;
	return (double)(_idum_p5 = (_idum_p5 * 430 + 2531) % 11979) * d;
}

double prnd6(void)
{
	static double d = 1. / 14406.;
	return (double)(_idum_p6 = (_idum_p6 * 967 + 3041) % 14406) * d;
}

double prnd7(void)
{
	static double d = 1. / 29282.;
	return (double)(_idum_p7 = (_idum_p7 * 419 + 6173) % 29282) * d;
}

double prnd8(void)
{
	static double d = 1. / 53125.;
	return (double)(_idum_p8 = (_idum_p8 * 171 + 11213) % 53125) * d;
}

double prnd9(void)
{
	static double d = 1. / 12960.;
	return (double)(_idum_p9 = (_idum_p9 * 1741 + 2731) % 12960) * d;
}

double prnd10(void)
{
	static double d = 1. / 14000.;
	return (double)(_idum_p10 = (_idum_p10 * 1541 + 2957) % 14000) * d;
}

double prnd11(void)
{
	static double d = 1. / 21870.;
	return (double)(_idum_p11 = (_idum_p11 * 1291 + 4621) % 21870) * d;
}

double prnd12(void)
{
	static double d = 1. / 31104.;
	return (double)(_idum_p12 = (_idum_p12 * 625 + 6571) % 31104) * d;
}

double prnd13(void)
{
	static double d = 1. / 139968.;
	return (double)(_idum_p13 = (_idum_p13 * 205 + 29573) % 139968) * d;
}

double prnd14(void)
{
	static double d = 1. / 29282.;
	return (double)(_idum_p14 = (_idum_p14 * 1255 + 6173) % 29282) * d;
}

double prnd15(void)
{
	static double d = 1. / 81000.;
	return (double)(_idum_p15 = (_idum_p15 * 421 + 17117) % 81000) * d;
}

double prnd16(void)
{
	static double d = 1. / 134456.;
	return (double)(_idum_p16 = (_idum_p16 * 281 + 28411) % 134456) * d;
}

double prnd17(void)
{
	static double d = 1. / 86436.;
	return (double)(_idum_p17 = (_idum_p17 * 1093 + 18257) % 86436) * d;
}

double prnd18(void)
{
	static double d = 1. / 121500.;
	return (double)(_idum_p18 = (_idum_p18 * 1021 + 25673) % 121500) * d;
}

double prnd19(void)
{
	static double d = 1. / 259200.;
	return (double)(_idum_p19 = (_idum_p19 * 421 + 54773) % 259200) * d;
}

double prnd20(void)
{
	static double d = 1. / 117128.;
	return (double)(_idum_p20 = (_idum_p20 * 1277 + 24749) % 117128) * d;
}

double prnd21(void)
{
	static double d = 1. / 121500.;
	return (double)(_idum_p21 = (_idum_p21 * 2041 + 25673) % 121500) * d;
}

double prnd22(void)
{
	static double d = 1. / 312500.;
	return (double)(_idum_p22 = (_idum_p22 * 741 + 66037) % 312500) * d;
}

double prnd23(void)
{
	static double d = 1. / 145800.;
	return (double)(_idum_p23 = (_idum_p23 * 3661 + 30809) % 145800) * d;
}

double prnd24(void)
{
	static double d = 1. / 175000.;
	return (double)(_idum_p24 = (_idum_p24 * 2661 + 36979) % 175000) * d;
}

double prnd25(void)
{
	static double d = 1. / 233280.;
	return (double)(_idum_p25 = (_idum_p25 * 1861 + 49297) % 233280) * d;
}

double prnd26(void)
{
	static double d = 1. / 244944.;
	return (double)(_idum_p26 = (_idum_p26 * 1597 + 51749) % 244944) * d;
}

double prnd27(void)
{
	static double d = 1. / 139968.;
	return (double)(_idum_p27 = (_idum_p27 * 3877 + 29573) % 139968) * d;
}

double prnd28(void)
{
	static double d = 1. / 214326.;
	return (double)(_idum_p28 = (_idum_p28 * 3613 + 45289) % 214326) * d;
}

double prnd29(void)
{
	static double d = 1. / 714025.;
	return (double)(_idum_p29 = (_idum_p29 * 1366 + 150889) % 714025) * d;
}

double prnd30(void)
{
	static double d = 1. / 134456.;
	return (double)(_idum_p30 = (_idum_p30 * 8121 + 28411) % 134456) * d;
}

double prnd31(void)
{
	static double d = 1. / 259200.;
	return (double)(_idum_p31 = (_idum_p31 * 7141 + 54773) % 259200) * d;
}

double prnd32(void)
{
	static double d = 1. / 233280.;
	return (double)(_idum_p32 = (_idum_p32 * 9301 + 49297) % 233280) * d;
}

double prnd33(void)
{
	static double d = 1. / 714025.;
	return (double)(_idum_p33 = (_idum_p33 * 4096 + 150889) % 714025) * d;
}

void init_irnd(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
}

int irnd(int max, int min)
{
	return min + (int)(_rand() * (double)(max - min + 1));
}

void init_gaussd(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
}

double gaussd(double av, double st)
{
	static double r1, r2, s;
	static int sw = 1;

	sw = 1 - sw;
	if(sw)	return av + st * r2 * s;
	do
	{
		r1 = 2. * _rand() - 1.;
		r2 = 2. * _rand() - 1.;
		s = r1 * r1 + r2 * r2;
	} while(s > 1.);
	s = sqrt(-2. * log(s) / s);
	return av + st * r1 * s;
}

void init_tri_rand(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
	_tri = _rand();
}

double tri_rand(void)
{
	double d;

	d = _tri;
	return (_tri = _rand()) - d;
}

void init_exp_rand(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
}

double exp_rand(void)
{
	return - log(_rand());
}

void init_binormal(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
}

void binormal(double r, double *x, double *y)
{
	double r1, r2, s;

	do
	{
		r1 = 2. * _rand() - 1.;
		r2 = 2. * _rand() - 1.;
		s = r1 * r1 + r2 * r2;
	} while(s > 1.);
	s = - log(s) / s;
	r1 *= sqrt((1 + r) * s);
	r2 *= sqrt((1 - r) * s);
	*x = r1 + r2;
	*y = r1 - r2;
	return;
}

void init_rnd_vector(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
}

void rnd_vector(double *x, double *y)
{
	double r;

	do
	{
		*x = 2. * _rand() - 1.;
		*y = 2. * _rand() - 1.;
		r = (*x * *x + *y * *y);
	} while(r > 1);
	r = sqrt(r);
	*x /= r;
	*y /= r;
	return;
}

void init_rnd_vector3(void (*_init)(long s), double (*_rnd)(void), long seed)
{
	_init(seed);
	_rand = _rnd;
}

void rnd_vector3(double *x, double *y, double *z)
{
	double t;

	do
	{
		*x = 2. * _rand() - 1.;
		*y = 2. * _rand() - 1.;
		*z = (*x * *x + *y * *y);
	} while(*z > 1);
	t = 2. * sqrt(1. - *z);
	*x *= t;
	*y *= t;
	*z = 2. * *z - 1.;
	return;
}