C多线程编程- 近似求解π

本程序使用蒙特卡洛方法估算圆周率（π）。它首先创建了指定数量的线程，每个线程生成一个随机点并检查该点是否在单位圆内。基于这些线程的结果，程序计算在单位圆内的点的比例，并乘以4来估算π的值。为了对比，程序还直接在主线程中（没有并发）进行了相同的π估算过程（由于每次都是生成随机数，所以这个基准也没啥意义hh~）。最后，程序打印出两种方法得到的π值。

#include 
#include 
#include 

struct arg_t {
    float x;
    float y;
};

void *start(void *arg) {
    struct arg_t *ll = (struct arg_t *)arg;
    float x1 = ll->x;
    float y1 = ll->y;
    long M = 0;
    if (x1*x1 + y1*y1 <= 1.0) {
        M ++;
    }
    pthread_exit((void *)M);
}

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "Please provide a number as an argument.\n");
        exit(1);
    }
    long N1 = atol(argv[1]);
    // printf("%ld\n", N1);
    pthread_t tids[N1];

    // pai(concurrency)
    for (long i = 0; i < N1; ++ i) {
        struct arg_t *arg = malloc(sizeof(struct arg_t));
        int x = rand();
        int y = rand();
        arg->x = 1.0*x / RAND_MAX;
        arg->y = 1.0*y / RAND_MAX;
        pthread_create(&tids[i], 0, start, arg);
    }

    void *res = 0;
    long M = 0;
    for (long i = 0; i < N1; ++ i) {
        pthread_join(tids[i], &res);
        M += (long)res;
    }
    printf("pai = %f\n", 4.0*M/N1);  // concurrency

    // pai(oracle)
    M = 0;
    for (long i = 0; i < N1; ++ i) {
        int x = rand();
        int y = rand();
        float x1 = 1.0*x / RAND_MAX;
        float y1 = 1.0*y / RAND_MAX;
        if (x1*x1 + y1*y1 <= 1.0) {
            M ++;
        }
    }
    printf("pai = %f\n", 4.0*M/N1);  // oracle
    pthread_exit(0);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

测试一下上述程序：

majn@tiger:~$ ./pai 10000
pai = 3.171200
pai = 3.142000
majn@tiger:~$ ./pai 100000
Segmentation fault (core dumped)
1
2
3
4
5

分析Segmentation fault的原因：

在上面的程序中，为每个线程都动态分配了 arg_t 结构的内存，但在线程执行完毕后，这些内存并没有被释放。虽然这不是立即的问题，但长期这样会导致内存泄露。应该在线程函数中或 pthread_join 之后释放这些内存。

改进方案：

一种方法是在线程函数 start 的结尾释放它（插个眼hh~）。但由于在主函数中可能还需要访问这些结构，更安全的方法是在 pthread_join 之后释放这些动态分配的内存。

修改上述代码：

#include 
#include 
#include 

struct arg_t {
    float x;
    float y;
};

void *start(void *arg) {
    struct arg_t *ll = (struct arg_t *)arg;
    float x1 = ll->x;
    float y1 = ll->y;
    long M = 0;
    if (x1*x1 + y1*y1 <= 1.0) {
        M ++;
    }
    pthread_exit((void *)M);
}

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "Please provide a number as an argument.\n");
        exit(1);
    }
    long N1 = atol(argv[1]);

    struct arg_t *args[N1];
    // printf("%ld\n", N1);
    pthread_t tids[N1];

    // pai(concurrency)
    for (long i = 0; i < N1; ++ i) {
        args[i] = malloc(sizeof(*args[i]));
        // struct arg_t *arg = malloc(sizeof(*arg));
        int x = rand();
        int y = rand();
        args[i]->x = 1.0*x / RAND_MAX;
        args[i]->y = 1.0*y / RAND_MAX;
        pthread_create(&tids[i], 0, start, args[i]);
    }

    void *res = 0;
    long M = 0;
    for (long i = 0; i < N1; ++ i) {
        pthread_join(tids[i], &res);
        M += (long)res;
        free(args[i]);
    }
    printf("pai = %f\n", 4.0*M/N1);  // concurrency

    // pai(oracle)
    M = 0;
    for (long i = 0; i < N1; ++ i) {
        int x = rand();
        int y = rand();
        float x1 = 1.0*x / RAND_MAX;
        float y1 = 1.0*y / RAND_MAX;
        if (x1*x1 + y1*y1 <= 1.0) {
            M ++;
        }
    }
    printf("pai = %f\n", 4.0*M/N1);  // oracle
    pthread_exit(0);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

这样，每次线程执行完毕并被主线程收回后，对应的动态分配的内存都会被释放。

测试一下修改后的程序：

majn@tiger:~$ ./pai 100000
pai = 3.721760
pai = 3.137640
majn@tiger:~$ ./pai 1000000
Segmentation fault (core dumped)
1
2
3
4
5

好好好，这样玩儿是吧。

分析Segmentation fault的原因：

仔细观察上面的程序，我使用了一个固定大小的线程数组：pthread_t tids[N1];。对于大的 N1 值，这可能会导致栈溢出。在大多数系统上，默认的栈大小可能不足以容纳大量的 pthread_t 变量。

改进方案：

考虑动态分配线程ID数组的空间。

#include 
#include 
#include 

struct arg_t {
    float x;
    float y;
};

void *start(void *arg) {
    struct arg_t *ll = (struct arg_t *)arg;
    float x1 = ll->x;
    float y1 = ll->y;
    long M = 0;
    if (x1*x1 + y1*y1 <= 1.0) {
        M ++;
    }
    pthread_exit((void *)M);
}

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "Please provide a number as an argument.\n");
        exit(1);
    }
    long N1 = atol(argv[1]);

    struct arg_t *args[N1];
    // printf("%ld\n", N1);
    pthread_t *tids = malloc(N1 * sizeof(pthread_t));

    // pai(concurrency)
    for (long i = 0; i < N1; ++ i) {
        args[i] = malloc(sizeof(*args[i]));
        // struct arg_t *arg = malloc(sizeof(*arg));
        int x = rand();
        int y = rand();
        args[i]->x = 1.0*x / RAND_MAX;
        args[i]->y = 1.0*y / RAND_MAX;
        pthread_create(&tids[i], 0, start, args[i]);
    }

    void *res = 0;
    long M = 0;
    for (long i = 0; i < N1; ++ i) {
        pthread_join(tids[i], &res);
        M += (long)res;
        free(args[i]);
    }
    printf("pai = %f\n", 4.0*M/N1);  // concurrency

    // pai(oracle)
    M = 0;
    for (long i = 0; i < N1; ++ i) {
        int x = rand();
        int y = rand();
        float x1 = 1.0*x / RAND_MAX;
        float y1 = 1.0*y / RAND_MAX;
        if (x1*x1 + y1*y1 <= 1.0) {
            M ++;
        }
    }
    printf("pai = %f\n", 4.0*M/N1);  // oracle
    pthread_exit(0);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

majn@tiger:~$ ./pai 1000000
pai = 3.972176
pai = 3.142136
majn@tiger:~$ ./pai 10000000
Segmentation fault (core dumped)
1
2
3
4
5

？？？不玩儿了！！！ 到此为止，我们把输入数据的规模N1从一开始最大接受10000扩大到现在最大接受1000000。

彩蛋： 要不然把为每个线程都动态分配的 arg_t 结构在线程函数 start 的结尾释放试一试？

#include 
#include 
#include 

struct arg_t {
    float x;
    float y;
};

void *start(void *arg) {
    struct arg_t *ll = (struct arg_t *)arg;
    float x1 = ll->x;
    float y1 = ll->y;
    long M = 0;
    if (x1*x1 + y1*y1 <= 1.0) {
        M ++;
    }
    free(arg);
    pthread_exit((void *)M);
}

int main(int argc, char **argv) {
    if (argc < 2) {
        fprintf(stderr, "Please provide a number as an argument.\n");
        exit(1);
    }
    long N1 = atol(argv[1]);
    // printf("%ld\n", N1);
    // pthread_t tids[N1];
    pthread_t *tids = malloc(N1 * sizeof(pthread_t));

    // pi(concurrency)
    for (long i = 0; i < N1; ++ i) {
        struct arg_t *arg = malloc(sizeof(struct arg_t));
        int x = rand();
        int y = rand();
        arg->x = 1.0*x / RAND_MAX;
        arg->y = 1.0*y / RAND_MAX;
        pthread_create(&tids[i], 0, start, arg);
    }

    void *res;
    long M = 0;
    for (long i = 0; i < N1; ++ i) {
        pthread_join(tids[i], &res);
        M += (long)res;
    }
    printf("pai = %f\n", 4.0*M/N1);  // concurrency

    // pi(oracle)
    M = 0;
    for (long i = 0; i < N1; ++ i) {
        int x = rand();
        int y = rand();
        float x1 = 1.0*x / RAND_MAX;
        float y1 = 1.0*y / RAND_MAX;
        if (x1*x1 + y1*y1 <= 1.0) {
            M ++;
        }
    }
    printf("pai = %f\n", 4.0*M/N1);  // oracle
    pthread_exit(0);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

majn@tiger:~$ ./pai_init 10000000
pai = 3.997218
pai = 3.141381
majn@tiger:~$ ./pai_init 100000000
pai = 3.999722
pai = 3.141420
majn@tiger:~$ ./pai_init 1000000000
^C   // 时间太长了，不想等了
1
2
3
4
5
6
7
8

Amazing！！！怎么这么神奇？

---

注： 蒙特卡洛方法（Monte Carlo method）是一种通过随机抽样来获得数值解的统计方法。这个方法得名于摩纳哥的蒙特卡洛赌场，因为它大量使用随机性和概率。蒙特卡洛方法在物理学、工程学、经济学和许多其他领域都有广泛的应用。

关键概念和特点：

1. 随机抽样：这是蒙特卡洛方法的核心。为了得到一个问题的数值解，这个方法使用随机数或更通常地说，使用伪随机数。

2. 统计结果：通过对大量的随机样本进行统计分析，得到的是一个近似解，而不是确切的解。

3. 精度与样本数量：通常，随着样本数量的增加，估算的精度也会提高。但是，为了使误差减少到原来的一半，样本数通常需要增加四倍。

4. 应用：蒙特卡洛方法在多种应用中都非常有用，尤其是在问题的解析解很难得到或者不存在时。例如，它被用于估算复杂积分、求解难以解析的统计物理问题、进行金融市场模拟等。

5. 示例 - 估算π：一个经典的应用是使用蒙特卡洛方法估算π的值。方法是这样的：随机投掷点到单位正方形内，统计落在单位圆内的点的数量。落在圆内的点数与总点数的比例，乘以4，就给出了π的近似值。

简而言之，蒙特卡洛方法是一种利用随机性来求解问题的技术，通过对大量样本的统计分析来获得结果。