not random enough for monte carlo

Question

I am trying to generate values from a normal distribution using a monte carlo method, as per the website http://math60082.blogspot.ca/2013/03/c-coding-random-numbers-and-monte-carlo.html

I modified the code a bit from the original so it calculates the variance and mean for the numbers generated directly to check if the method is working rather than do the tests separately (same difference really but just a heads up).

Question

Regardless of what I do, the variance is way above 1 and the mean is not zero. Is it possible the pseudo-random numbers generated aren't random enough?

Code

PLEASE NOTE THAT THE AUTHOR OF THE ABOVE GIVEN WEBSITE IS THE PERSON WHO WROTE THE CODE

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
// return a uniformly distributed random number
double uniformRandom()
{
    return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}

// return a normally distributed random number
double normalRandom()
{
    double u1=uniformRandom();
    double u2=uniformRandom();
    return cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1)); 
}

int main()
{
    double z;
    int N=1000;
    double array[N];
    double mean=0 ,variance=0;
    srand(time(NULL));

    for(int i=0;i<N;i++)
    {
        z=normalRandom();
        cout << i << "->"<< z<< endl;
        mean+=z;
        array[i]=z;
    }

    mean=mean/N ;
    cout << " mean = " << mean << endl;

    for(int i=0;i<N;i++)
    {
        variance = variance + (mean - array[i])*(mean - array[i]);
    }
    variance = variance/N;
    cout << " variance = " << variance << endl;

    return 0;
}

UPDATE

Apparently as pointed by users, I screwed up and the program was not working because of a very silly mistake.

Answer 1

You seems computed the mean in a wrong way. mean should be averaged over N, while you only sum over all array elements. current mean is actually sum.

mean = mean /N

Answer 2

rand() is a very low quality random numbers generator in most implementations. Some Linux versions would take value from kernel entropy pool, but it is not guaranteed across platforms (e.g. on Windows?) Use a Mersenne Twister instead. Boost libraries implement one.

EDIT: taocp answer highlights a coding problem, but the RNG issue still applies.