CUDA统计1-N素数总数异常问题求助

错误根源 #

你的代码在size≥9时结果错误，直接原因是浮点数精度问题：使用floor(pow(number, 0.5))计算平方根时，对于完全平方数（如9），pow可能返回略小于实际整数平方根的值（例如2.999999999而非3.0），导致floor后得到比正确值小1的结果。这会让循环提前终止，漏掉对平方根本身的检查，使得完全平方数被错误地判定为素数。

修复方案 #

1. 修正素数判断逻辑 #

用整数运算替代浮点数计算，直接通过i*i <= number判断是否需要继续循环，彻底避免浮点数精度问题：

for (int i = 2; i * i <= number; i++) {
    if (number % i == 0)
        return;
}

2. 优化线程/块配置 #

原代码每个线程对应一个块的配置效率极低，且无法支持超大size（如1e7）。改用256线程/块的标准配置，计算所需块数为(size + 255) / 256（向上取整），确保覆盖所有数字。

3. 移除不必要的数组 #

无需预先分配并拷贝数组到设备，每个线程可直接通过全局索引计算对应的数字idx + 1，大幅节省内存和数据传输开销。

修正后的完整代码 #

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <cstdio>

static void HandleError(cudaError_t err, const char* file, int line) {
    if (err != cudaSuccess) {
        printf("%s in %s at line %d\n", cudaGetErrorString(err), file, line);
        exit(EXIT_FAILURE);
    }
}
#define HANDLE_ERROR(err) (HandleError(err, __FILE__, __LINE__))

__global__ void countPrimeGPU(int size, int* count) {
    // 计算全局线程索引
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    // 跳过超出size的线程（因块数向上取整）
    if (idx >= size)
        return;
    int number = idx + 1; // 对应数字1~size

    if (number <= 1)
        return;
    // 用整数运算判断素数，避免浮点数误差
    for (int i = 2; i * i <= number; i++) {
        if (number % i == 0)
            return;
    }
    // 原子操作累加素数计数
    atomicAdd(count, 1);
}

int main() {
    const int size = 10000000; // 可修改为任意size，包括1e7
    int host_primeCount = 0;
    int* dev_primeCount;

    HANDLE_ERROR(cudaSetDevice(0));
    HANDLE_ERROR(cudaMalloc((void**)&dev_primeCount, sizeof(int)));
    HANDLE_ERROR(cudaMemcpy(dev_primeCount, &host_primeCount, sizeof(int), cudaMemcpyHostToDevice));

    // 配置线程块：256线程/块，计算所需块数
    const int block_size = 256;
    int grid_size = (size + block_size - 1) / block_size;
    countPrimeGPU<<<grid_size, block_size>>>(size, dev_primeCount);

    HANDLE_ERROR(cudaGetLastError());
    HANDLE_ERROR(cudaDeviceSynchronize());

    HANDLE_ERROR(cudaMemcpy(&host_primeCount, dev_primeCount, sizeof(int), cudaMemcpyDeviceToHost));
    printf("Prime count for the first %d numbers: %d\n", size, host_primeCount);

    cudaFree(dev_primeCount);
    HANDLE_ERROR(cudaDeviceReset());
    return 0;
}

进一步优化建议 #

对于超大size（如1e7），全局原子操作atomicAdd会成为性能瓶颈。可以通过块级共享内存先统计每个块内的素数数量，再将块结果原子累加到全局计数，减少原子操作的次数，大幅提升性能。

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