并发程序对 PI 的错误近似

Question

我正在尝试解决这个并发编程问题：

以原点为中心的单位圆上的点由函数 f(x) = sqrt(1-x²) 定义。回想一下，圆的面积是 pi*r²，其中 r 是半径。第 1 讲中描述的自适应求积例程可以近似 pi 值，计算单位圆右上象限的面积，然后将结果乘以 4。开发一个多线程程序（使用 Pthreads 的 C 语言或 Java 语言）来计算使用给定数量的进程（线程） np 的给定 epsilon 的 pi 被假定为命令行参数（即，它是给定的）。

但是，当使用越来越小的梯形来覆盖所描述的区域并用处理器连续划分问题时，我很难获得 pi 的正确近似值（我得到 2.666...）。

注意：代码中，worker指的是处理器。

这是我的代码：

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

我错过了什么？

谢谢！

Answer 1

你的被积数

f

是错误的。对于单位圆，我们有

x*x + y*y == 1

，你的

f

反而计算

y = 1-x*x

，这意味着你实际计算的是

(1-x*x)

从

x=0

到

x=1

的积分的4倍。积分可以解析为

g(1)-g(0)

，其中

g(x) = x - x*x*x/3

。这就是为什么你会得到

4*(1-1/3)=8/3=2.6666

。

您需要通过集成正确的功能来修复它。由于您正在积分第一象限，因此

x*x + y*y == 1

会转换为

y=sqrt(1-x*x)

，因此您应该修改

f

以返回它。

实现它的一种数值稳定的方法是使用

(1-x)*(1+x)

来计算

(1-x*x)

。例如，如果

x

仅比

sqrt(DBL_EPSILON)

大一点，则在计算

x*x

时，

1-x*x

中的大部分有效数字将丢失。相反，在

(1-x)*(1+x)

中，两个因子相乘之前，大部分有效数字都保留在其中。所以你应该将

f

修改为以下内容，

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