mch3.c

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#define NR 10001

double rnd(void)
{
  return rand()/(RAND_MAX+1.0);
}
int main(void)
{
        int n;
        double p[NR] = {0};
        double xo = 0;
        double xs = 0;
        double xx = 0.0;
        double k[NR] = {0};
        for(n = 0; n <= NR; n++) {
                xx = ((n > 1) ? p[n-1] : 0) + 0.5 - rnd() ;
                if(rnd() < exp(-0.5*(xx*xx - p[n] * p[n])))
                        p[n] = xx;
                else if (n > 1)
                        p[n] = p[n-1]; 
                xs += p[n-1]/NR; 
        }
        for(n = 0; n < NR; n++) {
                xo += (p[n] - xs)*(p[n-1] - xs);
        }
        for (n=0; n <= NR; n++) {
                k[n] = (p[2]*p[2+n])/(p[2]*p[2]);
        }

        for(n = 0; n < NR; n++) {
                printf ("%d %g %lf %lf %g\n", n, p[n], xs, sqrt(xo/(NR - 1)), k[n]);
        } 

        return 0;
}