并发程序PI近似值计算错误问题排查

原点为圆心的单位圆上的点由函数f(x) = sqrt(1-x²)定义。圆的面积公式为pi*r²（r为半径）。利用自适应积分法，通过计算单位圆右上象限的面积再乘以4来近似PI值。开发一个多线程程序（C语言Pthreads），根据给定的epsilon和命令行参数指定的线程数np计算PI。

但在使用梯形覆盖区域并通过线程拆分问题时，得到的PI近似值为2.666...，与正确值不符。

以下是我的代码：

#ifndef _REENTRANT
#define _REENTRANT
#endif

#include <pthread.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>

#define MAXWORKERS 3        // Maximum number of additional allowed workers besides the main thread.

const double EPSILON = 10e-10;
int numWorkers = 0;
double startTime, endTime;
int numCreatedThreads = 1;

// Mutex lock for numCreatedThreads.
pthread_mutex_t createdThreadsLock;

struct Info {
    double a, b;     // Extreme points.
    double fa, fb;   // Evaluation of extreme points.
    double area;     // Area to be approximated.
};

double read_timer() {
    static bool initialized = false;
    static struct timeval start;
    struct timeval end;
    if(!initialized) {
        gettimeofday(&start, NULL);
        initialized = true;
    }
    gettimeofday(&end, NULL);
    return (end.tv_sec - start.tv_sec) + 1.0e-6 * (end.tv_usec - start.tv_usec);
}

void readCommandLine(int argc, char* argv[]) {
    numWorkers = (argc > 1)? atoi(argv[1]) : MAXWORKERS;
    if (numWorkers > MAXWORKERS || numWorkers < 1)
        numWorkers = MAXWORKERS;
}

double f(double x) {
    return 1 - x*x;
}

void displayResults(double piApprox) {
    printf("\n===========================RESULTS===========================\n");
    printf("Pi was approximated up to %.0f decimal places: %f\n", -log10(EPSILON), 4 * piApprox);
    printf("The execution time is %g seconds.\n", endTime - startTime);
    printf("=============================================================\n");
}

void * calculatePI(void *args) {

    struct Info *info = args;

    // Calculate new data.
    double m = (info->a + info->b) / 2;
    double fm = f(m);
    double larea = (info->fa + fm) * (m - info->a) / 2;
    double rarea = (fm + info->fb) * (info->b - m) / 2;

    // Final result that will be returned.
    double *res = malloc(sizeof(double));

    // Check for termination.
    if(fabs(larea + rarea - info->area) < EPSILON) {
        *res = larea + rarea;
        return res;
    }

    // Boolean to control number of threads working.
    bool tooManyThreads = false;
    void* resultLeft, *resultRight;
    struct Info newLeft = {info->a, m, f(info->a), f(m), larea};
    struct Info newRight = {m, info->b, f(m), f(info->b), rarea};

    // Check whether we surpass the number of allowed threads.
    pthread_mutex_lock(&createdThreadsLock);
        tooManyThreads = numCreatedThreads + 1 > MAXWORKERS;
    pthread_mutex_unlock(&createdThreadsLock);

    if(!tooManyThreads) {

        // Update numCreatedThreads atomically.
        pthread_mutex_lock(&createdThreadsLock);
            numCreatedThreads++;
        pthread_mutex_unlock(&createdThreadsLock);

        // Execute left side with new thread, and right side with current thread.
        pthread_t newThread;
        pthread_create(&newThread, NULL, calculatePI, &newLeft);
        resultRight = calculatePI(&newRight);
        pthread_join(newThread, &resultLeft);
    }

    else {
        // Calculate the area of each side recursively on the same thread.
        resultLeft = calculatePI(&newLeft);
        resultRight = calculatePI(&newRight);
    }

    double * resultL = resultLeft;
    double * resultR = resultRight;
    *res = *resultL + *resultR;
    return res;
}

int main(int argc, char* argv[]) {

    // Read command line args if any.
    readCommandLine(argc, argv);

    // Initialize createdThreads mutex;
    pthread_mutex_init(&createdThreadsLock, NULL);

    // Initialize first information.
    struct Info info = {
            .a = 0,
            .b = 1,
            .fa = f(0),
            .fb = f(1),
            .area = (f(0) + f(1)) / 2
    };

    startTime = read_timer();

    // Obtain approximation of pi/4
    double * piApproximation = (double *) calculatePI(&info);

    endTime = read_timer();
    displayResults(*piApproximation);

    return 0;
}

请问我遗漏了什么？
谢谢！

错误原因分析 #

• 核心函数偏离需求：你实现的f(x)是1 - x*x，但题目要求的是单位圆右上象限的函数sqrt(1 - x²)。当前函数计算的是抛物线下方的面积，积分结果为∫₀¹(1-x²)dx = [x - x³/3]₀¹ = 2/3，乘以4后正好是8/3≈2.666...，这就是错误结果的根本原因。
• 线程计数逻辑不合理：MAXWORKERS定义的是除主线程外的最大线程数，但numCreatedThreads初始值为1（代表主线程），判断numCreatedThreads +1 > MAXWORKERS的逻辑会提前限制线程创建，比如MAXWORKERS=3时，最多允许3个额外线程，总线程数应为4，但当前逻辑会在numCreatedThreads=3时就拒绝创建新线程。
• 内存泄漏问题：每个递归调用都用malloc分配了double空间，但没有释放，会导致内存泄漏。

修正建议 #

1. 修复核心函数：将f(x)改为：

double f(double x) {
    return sqrt(1 - x*x);
}

2. 调整线程计数逻辑：
   • 把numCreatedThreads初始值设为0（仅统计额外创建的线程）
   • 判断条件改为numCreatedThreads >= numWorkers（用命令行传入的numWorkers而非MAXWORKERS做判断）
3. 添加内存释放：在合并左右结果后，释放resultL和resultR指向的内存：

*res = *resultL + *resultR;
free(resultL);
free(resultR);
return res;

4. 修正初始面积计算：初始area计算应为(fa + fb) * (b - a) / 2，当前代码少乘了区间长度，建议改为：

.area = (f(0) + f(1)) * (1 - 0) / 2

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