/*
 * Using random numbers for evaluationg of an integral
 *
 * Compile:
 *   clang -Wall -Weverything -O3 -o pi random.c -lm
 *
 * Run:
 *   ./pi
 */

#include <stdio.h>
#include <math.h>

void rseed (unsigned int seed);
unsigned int rand (void);
double random (void);

int main (void)
{
    long N = 1000L;
    double area = 0.;
    double x;

    rseed (3);

    for (int j = 0; j < 8; j++)
    {
	for (long i = 0; i < N; i++)
	{
	    x = random();
	    area += sqrt(1. - x*x);
	}
       
	area = 4*area/N;
	printf("%15ld    %10.6f    %10.8f\n", N, area, fabs(M_PI - area));
        N *= 10;
    }

    return 0;
}
/* 
 * Linear congruental generator
 *
 * X_{n+1} = (a*X_n + c) mod m
 *
 * m, 0 < m          the "modulus"
 * a, 0 < a < m      the "multiplier"
 * c, 0 <= c < m     the "increment"
 * X_0, 0 <= X_0 < m the "seed"
 *
 */

/*
 * rand()      generates a pseudo random integer 
 *             between 0 and 2^31-1 inclusive
 * random()    generates a pseudo random double 
 *             between 0. and 1. inclusive
 * rseed(seed) initialises the generator
 */

static unsigned int xrand = 1;
void rseed (unsigned int seed)
{
    xrand = seed;
}

unsigned int rand (void)
{
    xrand = xrand * 1103515245 + 12345;
    return (unsigned int) (xrand % 2147483648);
}

double random (void)
{
    return rand () / 2147483647.;
}