Monte Carlo Simulation of Percolation is not Giving Expected Results

Question

I've started working my way through the Princeton Algorithms course on Coursera. The course uses Java, but I decided to follow along with C as it is what I am most comfortable with. One of the assignments has you write a program to estimate the value of the percolation threshold via a Monte Carlo simulation (https://coursera.cs.princeton.edu/algs4/assignments/percolation/specification.php). I have written all the code, but my program's output is not the same as the one on the site (it is not completely off, but still incorrect.) e.g.

Their implementation:

~/Desktop/percolation> java-algs4 PercolationStats 200 100
mean                    = 0.5929934999999997
stddev                  = 0.00876990421552567
95% confidence interval = [0.5912745987737567, 0.5947124012262428]

~/Desktop/percolation> java-algs4 PercolationStats 2 100000
mean                    = 0.6669475
stddev                  = 0.11775205263262094
95% confidence interval = [0.666217665216461, 0.6676773347835391]

Mine:

~/percolation> ./percolation_stats 200 100
mean                    = 0.628564
stddev                  = 0.206286
95% confidence interval = [0.588132, 0.668996]

~/percolation> ./percolation_stats 2 100000
mean                    = 0.728548
stddev                  = 0.189745
95% confidence interval = [0.727371, 0.729724]

Here is my code:

percolation_stats.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "percolation.h"

double mean(double *, int);
double stddev(double *, int);
double confidencelo(double *, int);
double confidencehi(double *, int);

int main(int argc, char **argv) {
    int n, T, row, col;
    percolation *p;
    double *samp;

    srand((unsigned) time(NULL));
    sscanf(argv[1], "%d", &n);
    sscanf(argv[2], "%d", &T);
    samp = malloc(T * sizeof *samp);
    for (int i = 0; i < T; ++i) {
        p = creategrid(n);
        while (!percolates(p)) {
            do {
                row = rand() % n + 1;
                col = rand() % n + 1;
            } while (is_open(p, row, col));
            open(p, row, col);
        }
        samp[i] = (double) number_of_open_sites(p) / (n * n);
    }
    printf("mean                    = %g\n", mean(samp, T));
    printf("stddev                  = %g\n", stddev(samp, T));
    printf("95%% confidence interval = [%g, %g]\n", confidencelo(samp, T),
           confidencehi(samp, T));
    return 0;
}

// mean: sample mean of percolation threshold
double mean(double *a, int len) {
    double sum = 0.0;

    for (int i = 0; i < len; ++i) {
        sum += a[i];
    }
    return sum / len;
}

// stddev: sample standard deviation of percolation threshold
double stddev(double *a, int len) {
    double mean(double *, int);
    double sum = 0.0;
    double avg = mean(a, len);

    for (int i = 0; i < len; ++i) {
        sum += (a[i] - avg) * (a[i] - avg);
    }
    return sqrt(sum / (len - 1));
}

// confidencelo: low endpoint of 95% confidence interval
double confidencelo(double *a, int len) {
    double mean(double *, int);
    double stddev(double *, int);

    return mean(a, len) - (1.96 * stddev(a, len)) / sqrt(len);
}

// confidencehi: high endpoint of 95% confidence interval
double confidencehi(double *a, int len) {
    double mean(double *, int);
    double stddev(double *, int);

    return mean(a, len) + (1.96 * stddev(a, len)) / sqrt(len);
}

percolation.h

#include <stdbool.h>

typedef struct percolation percolation;

percolation *creategrid(int);
void open(percolation *, int, int);
bool is_open(percolation *, int, int);
bool is_full(percolation *, int, int);
int number_of_open_sites(percolation *);
bool percolates(percolation *);

percolation.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include "quick_union.h"

#define pos(p, x, y) (((x) - 1) * (p)->width + ((y) - 1))

typedef struct percolation {
    int width;
    int open_sites;
    UF *uf;
    bool *open;
} percolation;

// creategrid: creates n-by-n grid, with all sites initially blocked
percolation *creategrid(int n) {
    if (n <= 0) {
        fprintf(stderr, "open: illegal argument\n");
        exit(1);
    }
    percolation *p;

    p = malloc(sizeof *p);    
    p->width = n;
    p->open_sites = 0;
    p->uf = createuf(n * n + 2);
    p->open = malloc(n * n * sizeof *p->open);
    for (int i = 0; i < n * n; ++i) {
        p->open[i] = false;
    }
    for (int i = 0; i < n; ++i) {
        connect(p->uf, n * n, pos(p, 1, i));
        connect(p->uf, n * n + 1, pos(p, n, i));
    }
    return p;
}

// open: opens the site (row, col) if it is not open already
void open(percolation *p, int row, int col) {
    if (row < 1 || row > p->width || col < 1 || col > p->width) {
        fprintf(stderr, "open: illegal argument\n");
        exit(1);
    }
    if (p->open[pos(p, row, col)]) {
        return;
    }
    p->open[pos(p, row, col)] = true;
    ++p->open_sites;
    if (col > 1 && p->open[pos(p, row, col - 1)]) {
        connect(p->uf, pos(p, row, col), pos(p, row, col - 1));
    }
    if (col < p->width && p->open[pos(p, row, col + 1)]) {
        connect(p->uf, pos(p, row, col), pos(p, row, col + 1));
    }
    if (row > 1 && p->open[pos(p, row - 1, col)]) {
        connect(p->uf, pos(p, row, col), pos(p, row - 1, col));
    }
    if (row < p->width && p->open[pos(p, row + 1, col)]) {
        connect(p->uf, pos(p, row, col), pos(p, row + 1, col));
    }
}

// is_open: is the site (row, col) open?
bool is_open(percolation *p, int row, int col) {
    if (row < 1 || row > p->width || col < 1 || col > p->width) {
        fprintf(stderr, "is_open: illegal argument\n");
        exit(1);
    }
    return p->open[pos(p, row, col)];
}

// is_full: is the site (row, col) full?
bool is_full(percolation *p, int row, int col) {
    if (row < 1 || row > p->width || col < 1 || col > p->width) {
        fprintf(stderr, "is_full: illegal argument\n");
        exit(1);
    }
    return !p->open[pos(p, row, col)];
}

// number_of_open_sites: returns the number of open sites
int number_of_open_sites(percolation *p) {
    return p->open_sites;
}

// percolates: does the system percolate?
bool percolates(percolation *p) {
    return connected(p->uf, p->width * p->width, p->width + p->width + 1);
}

quick_union.h

#include <stdbool.h>

typedef struct UF UF;

UF *createuf(int);
void connect(UF *, int, int);
bool connected(UF *, int, int);
int find(UF *, int);

quick_union.c

#include <stdlib.h>
#include <stdbool.h>

#define max(a, b) (((a) > (b)) ? (a) : (b))

typedef struct UF {
    int count;
    int *id;
    int *sz;
    int *largest;
} UF;

// createuf: return pointer to UF with n elements
UF *createuf(int n) {
    UF *uf;

    uf = malloc(sizeof *uf);
    uf->count = n;
    uf->id = malloc(n * sizeof *uf->id);
    uf->sz = malloc(n * sizeof *uf->sz);
    uf->largest = malloc(n * sizeof *uf->largest);
    for (int i = 0; i < n; ++i) {
        uf->id[i] = uf->largest[i] = i;
        uf->sz[i] = 1;
    }
    return uf;
}

// connect: connect elements p and q
void connect(UF *uf, int p, int q) {
    int root(UF *, int);
    int i = root(uf, p);
    int j = root(uf, q);

    if (i == j) {
        return;
    }
    if (uf->sz[i] <= uf->sz[j]) {
        uf->id[i] = j;
        uf->sz[j] += uf->sz[i];
        uf->largest[j] = max(uf->largest[i], uf->largest[j]);
    } else {
        uf->id[j] = i;
        uf->sz[i] += uf->sz[j];
        uf->largest[i] = max(uf->largest[j], uf->largest[i]);
    }
}

// connected: return true if elements p and q are connected
bool connected(UF *uf, int p, int q) {
    int root(UF *, int);

    return root(uf, p) == root(uf, q);
}

// find: return largest element in i's connected component
int find(UF *uf, int i) {
    int root(UF *, int);

    return uf->largest[root(uf, i)];
}

// root: return root element of i
int root(UF *uf, int i) {
    while (i != uf->id[i]) {
        uf->id[i] = uf->id[uf->id[i]];
        i = uf->id[i];
    }
    return i;
}

Where did I go wrong?