/*

==========================================================================

normalmc.c - a program to run a simple normal-distribution Markov chain.
            (Designed for use in conjunction with jpar.c)

Copyright (c) 1999 by Jeffrey S. Rosenthal (jeff@math.toronto.edu).

Licensed for general copying, distribution and modification according to
the GNU General Public License (http://www.gnu.org/copyleft/gpl.html).

----------------------------------------------------

Compile with "cc normalmc.c -lm -o normalmc".

----------------------------------------------------

USAGE:  normalmc [runningtime [seed]]

SYNOPSIS:

normalmc will run for <runningtime> seconds (or 3 seconds if
no runningtime is given), using the drand48() pseudo-random number
generator with intial seed <seed> (or the current time in microseconds
if no seed is given).

During this time, it will continually run a Markov chain of the form
X_0 ~ N(0,1); X_{n+1} = N(X_n/2, 0.75) for n=0,1,2,...
It will compute the average of the values of (X_n)^2, i.e.
of the SQUARES of the values.

When done, it will discard any simulation in progress to avoid bias, and
then output two numbers: the average of the squared values considered,
and the number of such values averaged.

==========================================================================

*/

#define DEFNUMSECS 3   /* Default number of seconds to run; may be modified. */

#define pi 3.1415926

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

main(int argc, char *argv[])
{
    int i, N, runningsecs, tmpseed;
    double sqtot = 0.0;
    double x, normvar();
    long presenttime();
    long starttime, runningtime;

    /* Determine running time. */
    if ( (argc >= 2) && ( (runningsecs=atoi(argv[1])) > 0 ) )
        ;
    else
        runningsecs = DEFNUMSECS;
    runningtime = runningsecs * 1000000;

    /* Seed the pseudo-random number generator. */
    if ( (argc >= 3) && ( (tmpseed=atoi(argv[2])) > 0 ) )
        srand48(tmpseed);
    else
        seedrand();

    /* Prepare to begin. */
    starttime = presenttime();
    N = 0;
    x = normvar(0.0, 1.0);

    /* Loop, adding up Unif[0,1] random variables. */
    while (1) {
	x = normvar(x/2.0, 0.75);
	/* Discard simulation in progress if time up. */
        if ( (N>0) && (presenttime() > starttime + runningtime) )
            break;
	sqtot = sqtot + x*x;
	N = N + 1;
    }
	printf("%f     %d\n", sqtot/N, N);
}


double expvar(double m)
{
    return( log(1.0/drand48())/m );
}


double normvar(double m, double v)
{
    double stdnormvar(), sqrt();

    return( m + sqrt(v)*stdnormvar() );
}


double stdnormvar()
{
    return( sqrt( 2.0*log(1.0/drand48()) ) * cos(2.0*pi*drand48()) );
}


seedrand()
{
    struct timeval tmptv;
    int seed;

    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) (tmptv.tv_usec - 1000000 * (int) 
	( ( (double) tmptv.tv_usec ) / ( (double) 1000000.0 ) ) );
    srand48(seed);
    return(0);
}


long presenttime()  /* Return present time in microseconds.  */
{
    struct timeval tmp;
    gettimeofday (&tmp, (struct timezone *)NULL);
    return (tmp.tv_usec + 1000000 * tmp.tv_sec);
}