C++并行实现Trapezoidal Rule运行报Floating point exception错误

1 浮点异常直接诱因 #

你创建线程时传入的线程编号从0开始计数，第一个线程拿到的rank值为0，线程函数中执行stepSize = (1 - 0)/subInterval时触发除零错误，直接导致Floating point exception崩溃。

2 其他必须修复的问题 #

• 互斥锁未初始化：仅定义了pthread_mutex_t mutex，未调用pthread_mutex_init(&mutex, NULL)就直接加解锁，存在未定义行为
• 变量类型完全错误：步长stepSize、积分计算结果都是浮点值，你全部定义为long long整数类型，计算过程会丢失全部精度
• 并行逻辑错误：没有按线程拆分积分区间，每个线程独立计算完整积分没有并行意义，且计算结果未累加到全局sum变量，最终打印的sum是未初始化的无效值
• 宏定义无括号保护：#define f(x) 1/(1+pow(x,2))存在运算优先级隐患，调用时可能出现不符合预期的运算结果

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

const int MAX_THREADS = 1024;

long thread_count;
long long n;
double sum;
pthread_mutex_t mutex;

#define f(x) (1/(1+pow(x,2)))

void Usage(char* prog_name) {
   fprintf(stderr, "usage: %s <number of threads> <n> \n", prog_name);
   fprintf(stderr, "   n is the number of terms and should be >= 1\n");
   exit(0);
}  

void Get_args(int argc, char* argv[]) {
   if (argc != 3) Usage(argv[0]);
   thread_count = strtol(argv[1], NULL, 10);  
   if (thread_count <= 0 || thread_count > MAX_THREADS) Usage(argv[0]);
   n = strtoll(argv[2], NULL, 10);
   if (n <= 0) Usage(argv[0]);
}  

void* Trapezoidal_Rule(void* rank) {
    long my_rank = (long) rank;
    double h = 1.0 / n;
    long long local_n = n / thread_count;
    long long local_start = my_rank * local_n;
    long long local_end = (my_rank == thread_count - 1) ? n : (my_rank + 1) * local_n;
    
    double local_sum = f(local_start * h) + f(local_end * h);
    for (long long i = local_start + 1; i < local_end; i++) {
        local_sum += 2 * f(i * h);
    }
    local_sum = local_sum * h / 2.0;
    
    pthread_mutex_lock(&mutex);
    sum += local_sum;
    pthread_mutex_unlock(&mutex);

    return NULL;
}

int main(int argc, char* argv[]) {
    long thread;  
    pthread_t* thread_handles;
    
    Get_args(argc, argv);
    
    // 初始化互斥锁
    pthread_mutex_init(&mutex, NULL);
    thread_handles = (pthread_t *) malloc(thread_count * sizeof(pthread_t));
    
    sum = 0.0;
    for (thread = 0; thread < thread_count; thread++)
        pthread_create(&thread_handles[thread], NULL, Trapezoidal_Rule, (void*)thread);
    
    for (thread = 0; thread < thread_count; thread++)
        pthread_join(thread_handles[thread], NULL);
    
    printf("integration is = %.15f \n", sum);
    
    pthread_mutex_destroy(&mutex);
    free(thread_handles);
    return 0;
} 

运行./output 3 100000会输出接近0.785398163397448的结果，即π/4的近似值，符合积分预期。

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