frannor.c

来自「该程序是用vc开发的对动态数组进行管理的DLL」· C语言 代码 · 共 170 行

/* Copyright (c) Colorado School of Mines, 2003.*/
/* All rights reserved.                       */

/*********************** self documentation **********************/
/*****************************************************************************
FRANNOR - functions to generate a pseudo-random float normally distributed
		with N(0,1); i.e., with zero mean and unit variance.

frannor		return a normally distributed random float
srannor		seed random number generator for normal distribution

******************************************************************************
Function Prototypes:
float frannor (void);
void srannor (int seed);

******************************************************************************
frannor:
Input:		(none)
Returned:	normally distributed random float

srannor:
Input:
seed		different seeds yield different sequences of random numbers.

******************************************************************************
Notes:
Adapted from subroutine rnor in Kahaner,  Moler and Nash (1988)
which in turn was based on  an algorithm by 
Marsaglia and Tsang (1984).

******************************************************************************
References:
"Numerical Methods and Software", D.
Kahaner, C. Moler, S. Nash, Prentice Hall, 1988.

Marsaglia G. and Tsang, W. W., 1984,
A fast, easily implemented method for sampling from decreasing or symmetric
unimodal density functions:  SIAM J. Sci. Stat. Comput., v. 5, no. 2,
p. 349-359.

*****************************************************************************
Author:  Dave Hale, Colorado School of Mines, 01/21/89
*****************************************************************************/
/**************** end self doc ********************************/

#include "cwp.h"

/* constants for normal random number generator */
#define AA 12.37586F
#define B 0.4878992F
#define C 12.67706F
#define C1 0.9689279F
#define C2 1.301198F
#define PC 0.01958303F
#define XN 2.776994F
#define OXN 0.3601016F
#define NBITS 24
static float v[]={
	0.3409450f, 0.4573146f, 0.5397793f, 0.6062427f, 0.6631691f,
	0.7136975f, 0.7596125f, 0.8020356f, 0.8417227f, 0.8792102f, 0.9148948f,
	0.9490791f, 0.9820005f, 1.0138492f, 1.0447810f, 1.0749254f, 1.1043917f,
	1.1332738f, 1.1616530f, 1.1896010f, 1.2171815f, 1.2444516f, 1.2714635f,
	1.2982650f, 1.3249008f, 1.3514125f, 1.3778399f, 1.4042211f, 1.4305929f,
	1.4569915f, 1.4834526f, 1.5100121f, 1.5367061f, 1.5635712f, 1.5906454f,
	1.6179680f, 1.6455802f, 1.6735255f, 1.7018503f, 1.7306045f, 1.7598422f,
	1.7896223f, 1.8200099f, 1.8510770f, 1.8829044f, 1.9155830f, 1.9492166f,
	1.9839239f, 2.0198430f, 2.0571356f, 2.0959930f, 2.1366450f, 2.1793713f,
	2.2245175f, 2.2725185f, 2.3239338f, 2.3795007f, 2.4402218f, 2.5075117f,
	2.5834658f, 2.6713916f, 2.7769943f, 2.7769943f, 2.7769943f, 2.7769943f
};

/* internal state variables for uniform random number generator */
static int i=16,j=4;
static float u[]={
	0.8668672834288f,  0.3697986366357f,  0.8008968294805f,
	0.4173889774680f,  0.8254561579836f,  0.9640965269077f,
	0.4508667414265f,  0.6451309529668f,  0.1645456024730f,
	0.2787901807898f,  0.06761531340295f, 0.9663226330820f,
	0.01963343943798f, 0.02947398211399f, 0.1636231515294f,
	0.3976343250467f,  0.2631008574685f
};

/* macro to generate a random number uni uniform on [0,1) */
#define UNI(uni) \
	uni = u[i]-u[j]; if (uni<0.0) uni += 1.0; u[i] = uni; \
	if (--i<0) i = 16;  if (--j<0) j = 16

float
frannor(void)
/**************************************************************************
return a normally distributed random float
***************************************************************************
Returned:	normally distributed random float
**************************************************************************/
{
	int k;
	float uni,vni,rnor,x,y,s,bmbx,xnmx;
	
	/* uni is uniform on [0,1) */
	UNI(uni);
		
	/* vni is uniform on [-1,1) */
	vni = uni+uni-1.0F;
	
	/* k is in range [0,63] */
	k = ((int)(u[i]*128))%64;
	
	/* fast part */
	rnor = vni*v[k+1];
	if (ABS(rnor)<=v[k]) return rnor;
	
	/* slow part */
	x = (ABS(rnor)-v[k])/(v[k+1]-v[k]);
	UNI(y);
	s = x+y;
	if (s<=C2) {
		if (s<=C1) return rnor;
		bmbx = B-B*x;
		if (y<=C-AA*exp(-0.5*bmbx*bmbx)) {
			if (exp(-0.5*v[k+1]*v[k+1])+y*PC/v[k+1] <= 
				exp(-0.5*rnor*rnor)) return rnor;
			do {
				UNI(y);
				x = (float)(OXN*log(y));
				UNI(y);
			} while (-2.0*log(y)<=x*x);
			xnmx = XN-x;
			return (rnor>=0.0 ? ABS(xnmx) : -ABS(xnmx));
		}
	}
	bmbx = B-B*x;	
	return (rnor>=0.0 ? ABS(bmbx) : -ABS(bmbx));
}			

void
srannor (int seed)
/*****************************************************************************
seed random number generator
******************************************************************************
Input:
seed		different seeds yield different sequences of random numbers.
*****************************************************************************/
{
	int ii,jj,ia,ib,ic,id;
	float s,t;
	
	i = 16;
	j = 4;
	ia=ABS(seed)%32707;
	ib=1111;
	ic=1947;
	for (ii=0; ii<17; ii++) {
		s = 0.0;
		t = 0.5;
		for (jj=0; jj<64; jj++) {
			id = ic-ia;
			if (id<0) {
				id += 32707;
				s += t;
			}
			ia = ib;
			ib = ic;
			ic = id;
			t *= 0.5;
		}
		u[ii] = s;
	}
}