如何用C语言多线程计算Gregory-Leibniz级数n/m项之和

需求说明 #

• 从命令行参数获取两个整数m（线程数）和n（级数总项数），创建m个线程并行计算
• 每个线程负责计算Gregory-Leibniz级数的n/m项之和，级数公式：

  pi = 4 * (1 - 1/3 + 1/5 - 1/7 + 1/9 - ...)

• 线程分工：第一个线程计算第1到n/m项的和，第二个线程计算第n/m+1到2n/m项，以此类推
• 线程完成计算后需打印自身的部分和，并将该值原子性添加到共享全局变量中，同时确保所有m个线程都完成原子加法操作

问题现象 #

尝试了两段C语言代码，但PI的输出值不稳定：例如当m=16、n=1024时，有时输出3.125969，有时输出12.503874、15.629843或6.251937。

初始尝试代码 #

#include<stdio.h>
#include<pthread.h>
#include <stdlib.h>
#include<math.h>

pthread_barrier_t barrier;
int count;
long int term;
double total;

void *thread_function(void *vargp)
{
    int thread_rank = *(int *)vargp;
    
    pthread_barrier_wait(&barrier);
    
    double sum = 0.0;
    int n = count * term;
    int start = n - term;
    
    for(int i = start; i < n; i++) 
    {
        sum += pow(-1, i) / (2*i+1);
    }
    total += sum; 
    count++;
    
    printf("thr %d : %lf \n", thread_rank, sum);
    
    return NULL;
}

int main(int argc,char *argv[])
{
    if (argc <= 2) {
        printf("missing arguments. please pass two num. in arguments\n");
        exit(-1);
    }
    
    int m = atoi(argv[1]);
    int n = atoi(argv[2]);
    
    count = 1;
    term = n / m;

    pthread_t thread_id[m];
    int i, ret;
    double pi;
    
    pthread_barrier_init(&barrier, NULL, m);

    for(i = 0; i < m; i++)
    {
        ret = pthread_create(&thread_id[i], NULL , &thread_function, (void *)&i);
        if (ret) {
            printf("unable to create thread! \n");
            exit(-1);
        } 
    }
    
    for(i = 0; i < m; i++)
    {
        if(pthread_join(thread_id[i], NULL) != 0) {
            perror("Failed to join thread");
        }
    }
    
    pi = 4 * total;
    printf("%lf ", pi);
    
    pthread_barrier_destroy(&barrier);

    return 0;
}

修改后的代码 #

#include <inttypes.h>
#include <math.h>
#include <pthread.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

struct args {
    uint64_t thread_id;
    struct {
        uint64_t start;
        uint64_t end;
    } range;
};

pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
pthread_barrier_t barrier;
long double total = 0;
uint64_t total_iterations = 0;

void *partial_sum(void *arg)
{
    struct args *args = arg;
    long double sum = 0;
    
    printf("waiting for barrier in thread -> %" PRId64 "\n", args->thread_id);
    pthread_barrier_wait(&barrier);
    
    for (uint64_t n = args->range.start; n < args->range.end; n++) 
        sum += pow(-1.0, n) / (1 + n * 2);
    
    if (pthread_mutex_lock(&mutex)) {
        perror("pthread_mutex_lock");
        exit(EXIT_FAILURE);
    }
    
    total += sum;
    total_iterations += args->range.end - args->range.start;
    
    if (pthread_mutex_unlock(&mutex)) {
        perror("pthread_mutex_unlock");
        exit(EXIT_FAILURE);
    }
    
    printf("thr %" PRId64 " : %.20Lf\n", args->thread_id, sum);
    
    return NULL;
}

int main(int argc,char *argv[])
{
    if (argc <= 2) {
        fprintf(stderr, "usage: %s THREADS TERMS.\tPlease pass two num. in arguments\n", *argv);
        return EXIT_FAILURE;
    }
    
    int m = atoi(argv[1]);
    int n = atoi(argv[2]);

    if (!m || !n) {
        fprintf(stderr, "Argument is zero.\n");
        return EXIT_FAILURE;
    }
    
    uint64_t threads = m;
    uint64_t terms = n;
    
    uint64_t range = terms / threads;
    uint64_t excess = terms - range * threads;
    
    pthread_t thread_id[threads];
    struct args arguments[threads];
    
    int ret;
    
    ret = pthread_barrier_init(&barrier, NULL, m);
    if (ret) {
        perror("pthread_barrier_init");
        return EXIT_FAILURE;
    }
    
    for (uint64_t i = 0; i < threads; i++) {
        arguments[i].thread_id = i;
        arguments[i].range.start = i * range;
        arguments[i].range.end = arguments[i].range.start + range;
        
        if (threads - 1 == i)
            arguments[i].range.end += excess;
            
        printf("In main: creating thread %ld\n", i);
        ret = pthread_create(thread_id + i, NULL, partial_sum, arguments + i);
        if (ret) {
            perror("pthread_create");
            return EXIT_FAILURE;
        }
    }
    
    for (uint64_t i = 0; i < threads; i++)
        if (pthread_join(thread_id[i], NULL))
            perror("pthread_join");

    pthread_barrier_destroy(&barrier);

    printf("Pi value is : %.10Lf\n", 4 * total);
    printf("COMPLETE? (%s)\n", total_iterations == terms ? "YES" : "NO");
 
    return 0;   
}

问题根源分析 #

初始代码核心问题 #

1. 线程参数传递错误：创建线程时传递&i，循环变量i会被主线程持续修改，线程启动时可能读取到更新后的i值，导致多个线程的thread_rank重复或错误，引发未定义行为。
2. 共享变量无同步保护：total和count都是全局共享变量，total += sum和count++均为非原子操作，多线程同时读写会引发数据竞争，导致total累加错误、count计数混乱。
3. 计算范围依赖非原子变量：线程通过全局count计算任务范围，count的修改无同步，多个线程可能拿到相同的count值，导致重复计算同一区间或区间计算错误。

修改后代码潜在问题 #

1. pow函数精度误差：用pow(-1.0, n)计算符号位，不仅效率低，还会引入浮点精度误差，累积后导致部分和计算不准确。
2. 编译链接问题：若编译时未正确链接pthread库（需加-lpthread参数）和math库（需加-lm参数），会引发线程或数学函数的未定义行为，导致结果不稳定。

修复方案与优化代码 #

修复要点 #

1. 为每个线程单独分配参数，避免传递循环变量地址。
2. 用奇偶性判断替代pow函数计算符号位，消除精度误差。
3. 确保共享变量total的更新在互斥锁保护下完成。
4. 预先计算每个线程的任务范围，直接传递给线程，避免依赖全局变量动态计算。

修复后的代码 #

#include <inttypes.h>
#include <pthread.h>
#include <stdlib.h>
#include <stdio.h>

struct args {
    uint64_t thread_id;
    uint64_t start;
    uint64_t end;
};

pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
long double total = 0;
uint64_t total_iterations = 0;

void *partial_sum(void *arg)
{
    struct args *args = arg;
    long double sum = 0;

    for (uint64_t n = args->start; n < args->end; n++) {
        // 用奇偶性判断符号，替代pow函数避免精度误差
        long double sign = (n % 2 == 0) ? 1.0L : -1.0L;
        sum += sign / (1 + n * 2);
    }

    if (pthread_mutex_lock(&mutex)) {
        perror("pthread_mutex_lock");
        exit(EXIT_FAILURE);
    }

    total += sum;
    total_iterations += args->end - args->start;

    if (pthread_mutex_unlock(&mutex)) {
        perror("pthread_mutex_unlock");
        exit(EXIT_FAILURE);
    }

    printf("thr %" PRId64 " : %.20Lf\n", args->thread_id, sum);

    return NULL;
}

int main(int argc, char *argv[])
{
    if (argc <= 2) {
        fprintf(stderr, "用法: %s 线程数 总项数\n", *argv);
        return EXIT_FAILURE;
    }

    int m = atoi(argv[1]);
    int n = atoi(argv[2]);

    if (m <= 0 || n <= 0) {
        fprintf(stderr, "线程数和总项数必须大于0\n");
        return EXIT_FAILURE;
    }

    uint64_t threads = m;
    uint64_t terms = n;

    uint64_t range = terms / threads;
    uint64_t excess = terms - range * threads;

    pthread_t thread_id[threads];
    struct args arguments[threads];

    int ret;

    for (uint64_t i = 0; i < threads; i++) {
        arguments[i].thread_id = i;
        arguments[i].start = i * range;
        arguments[i].end = arguments[i].start + range;

        // 最后一个线程处理剩余的项
        if (i == threads - 1) {
            arguments[i].end += excess;
        }

        ret = pthread_create(&thread_id[i], NULL, partial_sum, &arguments[i]);
        if (ret) {
            perror("pthread_create");
            return EXIT_FAILURE;
        }
    }

    for (uint64_t i = 0; i < threads; i++) {
        if (pthread_join(thread_id[i], NULL)) {
            perror("pthread_join");
        }
    }

    pthread_mutex_destroy(&mutex);

    printf("PI计算值: %.10Lf\n", 4 * total);
    printf("计算完成? (%s)\n", total_iterations == terms ? "是" : "否");

    return 0;
}

编译与运行命令 #

# 编译（需链接pthread和math库）
gcc -o pi_calc pi_calc.c -lpthread -lm

# 运行示例
./pi_calc 16 1024

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