🚀 go-pugleaf

Siew <siew@mx4.redestb.es> wrote:
: Hola a todos,
: �existe alg�n m�todo entre las librer�as de C incluidas en Linux que 
: permita obtener n�mero aleatorios, como el m�todo rand(), pero siguiendo 
: una distribuci�n normal? En el caso de que no existiera, �d�nde podr�a 
: encontrar informaci�n sobre como programarla?

/*
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <sys/time.h>
#include <math.h>

#define	sqr(a)	((a)*(a))
double f ( double x, double m, double s2 );

double	coso[1024];

int main( int argc, char *argv[] )
{
	int		i, n;
	double	s2, m, x, sum, med;
	struct timeval  tp;
	int plot;


	m = (double) atof( argv[1] );
	s2 = (double) atof( argv[2] );
	n = (int) atoi( argv[3] );
	plot = (int) atoi( argv[4] );

#if defined _POSIX_VERSION || defined _HPUX_SOURCE
	gettimeofday(&tp, NULL);
#else
	gettimeofday(&tp);
#endif

	srand48((long)tp.tv_usec);
	for( i = 0; (i<n) && (i<sizeof(coso)/sizeof(double)); )
	{
		x = (drand48() -0.5) *  1000.0;
		if ( drand48() < f(x, m, s2) )
			coso[i++] = x;
	}
	n = i;

	for( sum = 0, i=0; i<n ; i++ )
		sum += coso[i];
	med = sum / (double)n;

	if (plot)
	{
		for( i=0; i<n ; i++ )
			printf("%d %g\n",i, coso[i]);
	}

	fprintf( stderr, "m: %g\t", med);

	for( sum = 0.0, i=0; i<n ; i++ )
		sum += sqr (coso[i] - med);

	fprintf( stderr, "s: %g %d\n", sum / (double) n, n);
	
	return 0;
}

/**************************************************************************/

/*                            (x - m)^2
 *              1          - -----------
 * f(x) = -------------- e     2 s2
 *        s (2 M_PI)^0.5
 */
double f ( double x, double m, double s2 )
{
	x -= m;
	x *= x;
	x /= (2.0 * s2);
	x = exp(-x) / sqrt( 2.0 * M_PI * s2 );
	return x;
}

/*
 *    / x
 *    |       f(t)     dt
 *   / -inf       0,1
 */
double F ( double x )
{
	double t;

	t = 1.0 / (1 + 0.2316419 * x);

	return(
		1 - f(x, 0.0, 1.0) * (
		(((
		    1.330274429 * t
		  - 1.821255978) *
		t + 1.78147937) *
		t - 0.356563782) *
		t + 0.31938153) * t
		);
}

double u ( double x, double min, double max )
{
	if (( x <= min ) || ( x >= max )) return 0;
	else
		return ( 1 / ( max - min ) );
}

-- 
PGP and other useless info at      \
http://www.cps.unizar.es/~spd/      \
finger://daphne.cps.unizar.es/spd    \       Timeo Danaos et dona ferentes
ftp://ivo.cps.unizar.es/pub/          \                         (Virgilio)