C语言实现渗流蒙特卡洛模拟结果与预期不符排查

架构老林

2026-8-30

问题背景

我正在学习Coursera平台开设的普林斯顿大学算法课程，课程教学使用Java语言，我选择使用自己更熟悉的C语言跟随练习。其中一项作业要求编写程序，通过*Monte Carlo simulation（蒙特卡洛模拟）*估算渗流阈值，作业规范地址：作业规范页。
我已经完成了全部代码编写，但程序输出和课程提供的参考实现结果不一致（结果并非完全错误，但存在明显偏差），运行结果对比如下：

运行结果对比

课程参考实现运行结果

~/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]

个人实现运行结果

~/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]

完整代码

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;
}

错误定位

代码共有4处核心错误，直接导致结果偏差：

creategrid中pos宏传参错误，引发内存越界+连通逻辑混乱
pos宏的行、列参数要求是1-based（从1开始计数），但初始化连接虚拟节点的循环中，列参数直接传了0-based的循环变量i：

i=0时列参数为0，计算出的索引为-1，connect操作越界写入堆内存，破坏程序数据
所有顶行、底行的连接整体偏移一列，最后一列格子未被连接，底虚拟节点还被错误连到顶行格子，完全打乱连通性判断。
正确写法为列参数传i+1：

for (int i = 0; i < n; ++i) {
    connect(p->uf, n * n, pos(p, 1, i+1));
    connect(p->uf, n * n + 1, pos(p, n, i+1));
}

percolates函数底部虚拟节点索引计算错误
定义的顶部虚拟节点索引为n*n，底部虚拟节点索引为n*n+1，但函数中错误写为p->width + p->width +1，n>2时该值和真实索引完全不匹配（如n=200时真实索引为40001，错误值为401，对应网格内普通格子），导致需要打开远多于正常值的格子才能触发渗流判断，均值偏高。
正确写法：

bool percolates(percolation *p) {
    return connected(p->uf, p->width * p->width, p->width * p->width + 1);
}

初始化阶段提前连接所有顶/底行格子引发回流问题
初始化时不管格子是否打开，就将所有顶行格子连到顶部虚拟节点、底行格子连到底部虚拟节点。Union-Find不区分格子打开状态，初始化连接后同一行格子会通过虚拟节点提前连通，还会出现回流问题：底部连通的格子通过虚拟节点反向连通到未和顶部连通的区域，导致判断失真。
正确逻辑是初始化时不提前连接任何格子到虚拟节点，仅在格子被open时，如果它在第一行则连接顶部虚拟节点，在最后一行则连接底部虚拟节点。
is_full函数实现不符合定义
按照作业要求，full指「打开的格子能够连通到顶部」，当前实现将未打开的格子判定为full，和定义完全相反。正确实现需要先判断格子是否打开，再判断是否和顶部虚拟节点连通。

内容的提问来源于stack exchange，提问作者vim_overlord

C语言实现渗流蒙特卡洛模拟结果与预期不符排查

运行结果对比 #

课程参考实现运行结果 #

个人实现运行结果 #

完整代码 #

percolation_stats.c #

percolation.h #

percolation.c #

quick_union.h #

quick_union.c #