百万数据sum太慢_给大家表演一个n方过百万（利用GPU暴力计算多项式乘法）

（学生的小打小闹，请各位炼丹大佬绕行）

在2019年noi冬令营上，我们学会（fei）了如何利用cpu指令加速暴力，并达到n方过百万。今天，我在这里表演一个利用gpu暴力算法过百万级数据的多项式乘法。

感兴趣的可以自行百度cuda。

源文件1

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

#include <stdio.h>
#include <iostream>
#include<windows.h>
/*懂的都懂*/
const int P = 998244353;
/*生成测试数据*/
void createData();
/*得到分块信息*/
void createDim(dim3& gridDim, dim3& blockDim);
/*为GPU分配空间*/
void allocDevice();
/*将数据传输到GPU*/
void copyFHTD();
/*核函数*/
__global__ void kernel(int* da, int* db, int* dc, int n);
/*将数据传输到CPU*/
void copyFDTH();
/*释放GPU空间*/
void freeDevice();
/*随机检查结果，若发现错误则返回错误在结果多项式中的下标。若没有错误返回-1*/
int check();
/*输出结果，仅用于调试时。然后清理申请的空间*/
void output();

/*GPU上储存数据*/
int* da, * db, * dc;
/*原多项式为n-1次，结果为n*2-1次*/
int n;
/*内存上储存数据*/
int* a, * b, * c;
/*计时用*/
long long startTime;
/*输出提示信息和时间*/
void outputTime(const char *str) {
    std::cout << str << ":" << GetTickCount64() - startTime << std::endl;
}
int main(){
    startTime = GetTickCount64();
    createData();
    dim3 gridDim, blockDim;
    createDim(gridDim, blockDim);
    allocDevice();
    copyFHTD();
    outputTime("start");
    kernel << <gridDim, blockDim >> > (da,db,dc,n);
    cudaDeviceSynchronize();
    outputTime("finish");
    copyFDTH();
    freeDevice();
    std::cout << check() << std::endl;
    output();
    return 0;
}

源文件2

//#define DEBUG

#include<cstdio>
#include<cstdlib>
#include<ctime>
#include<iostream>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"


extern int n,*a,*b,*c;
extern int* da, * db, * dc;
extern const int P=998244353;

void createData() {
    srand(time(0));
#ifdef DEBUG
    n = 4;
#else
    n = 500000;
#endif
    a = new int[n];
    b = new int[n];
    c = new int[n * 2ll - 1];
    for (int i = 0; i < n; i++) {
        a[i] = rand() % 10;
        b[i] = rand() % 10;
    }
#ifdef DEBUG
    for (int i = 0; i < n; i++)
        printf("%d ", a[i]);
    printf("n");
    for (int i = 0; i < n; i++)
        printf("%d ", b[i]);
    printf("n");
#endif
}
void createDim(dim3& gridDim, dim3& blockDim) {
    blockDim.x = 32;
    int l = n * 2 - 1;
    gridDim.x = l / 32 + bool(l % 32);
}
void allocDevice() {
    cudaMalloc(&da, sizeof(int) * n);
    cudaMalloc(&db, sizeof(int) * n);
    cudaMalloc(&dc, sizeof(int) * (n * 2ll - 1));
}
void copyFHTD() {
    cudaMemcpy(da, a, sizeof(int) * n, cudaMemcpyHostToDevice);
    cudaMemcpy(db, b, sizeof(int) * n, cudaMemcpyHostToDevice);
    cudaMemcpy(dc, c, sizeof(int) * (n * 2ll - 1), cudaMemcpyHostToDevice);
}
__global__ void kernel(int *da,int *db,int *dc,int n) {
    int id = blockDim.x * blockIdx.x + threadIdx.x;
#ifdef DEBUG
    if (id == 0) {
        printf("n:%dn", n);
        for (int i = 0; i < n; i++)
            printf("%d ", da[i]);
        printf("n");
        for (int i = 0; i < n; i++)
            printf("%d ", db[i]);
        printf("n");
    }
#endif
    int sum = 0;
    /*第一个判断对GPU运行效率无影响。
    因为对于每一个线程来说，要么进去，要么不进去，而且一定会有线程进去。*/
    if (id < n * 2 - 1) {
        /*第二个判断有很大几率会导致双倍用时。
        对于一个warp，如果有的线程进入了第一个判断而其他线程进入了第二个判断，
        则会导致必须先等待第一个判断中的线程执行完毕，再执行第二个判断中的线程。
        在划分block时可以规避这个问题，但编码较为复杂我选择双倍。
        */
        if (id < n) {
            for (int i = 0; i <= id; i++) {
                sum += da[i] * db[id - i];
                sum %= P;
            }
        } else {
            for (int i = id - n + 1; i <= id; i++) {
                sum += da[i] * db[id - i];
                sum %= P;
            }
        }
        dc[id] = sum;
#ifdef DEBUG
        printf("%d:%dn", id, sum);
#endif
    }
}
void copyFDTH() {
    cudaMemcpy(c, dc, sizeof(int) * (n * 2ll - 1), cudaMemcpyDeviceToHost);
}
void freeDevice() {
    cudaFree(da);
    cudaFree(db);
    cudaFree(dc);
}
void output() {
#ifdef DEBUG
    for (int i = 0; i < n * 2ll - 1; i++)
        printf("%d ", c[i]);
    printf("n");
#endif
    delete a;
    delete b;
    delete c;
    a = b = c = 0;
}
int check() {
    int k,sum;
    for (int _ = 0; _ < 100; _++) {
        k = rand() % (n * 2 - 1);
        sum = 0;
        for (int i = 0; i <= k; i++) {
            if (k - i < n && i < n) {
                sum = (sum + a[i] * b[k - i]) % P;
            }
        }
        if (sum != c[k]) {
            std::cout << "need:" << sum << ";get:" << c[k] << std::endl;
            return k;
        }

    }
    return -1;
}

可以看出，再kernel函数中，我选择了最最最最暴力的方法进行计算。而且没有加入shared优化在n=5e5，即结果数组1e6的数据下跑出了1.3s的成绩，足以和NTT媲美。

RTX 2060牛逼！