高性能计算：阶乘

1
2
3
4
5
6
7

int factorial(int n) {
    int result = 1;
    for (int i = 2; i <= n; i++) {
        result *= i;
    }
    return result;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

class HugeInteger {
private:
    vector<int> digits;  // 存储每一位数字
    
public:
    HugeInteger operator*(const HugeInteger& other) const {
        // 实现竖式乘法
        HugeInteger result;
        result.digits.resize(digits.size() + other.digits.size(), 0);
        
        for (size_t i = 0; i < digits.size(); i++) {
            for (size_t j = 0; j < other.digits.size(); j++) {
                result.digits[i + j] += digits[i] * other.digits[j];
                result.digits[i + j + 1] += result.digits[i + j] / 10;
                result.digits[i + j] %= 10;
            }
        }
        
        return result;
    }
};

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

HugeInteger operator*(const HugeInteger& other) const {
    HugeInteger result;
    result.digits.resize(digits.size() + other.digits.size(), 0);
    
    for (size_t i = 0; i < digits.size(); i++) {
        long long carry = 0;
        for (size_t j = 0; j < other.digits.size() || carry; j++) {
            long long cur = result.digits[i + j] 
                          + digits[i] * (j < other.digits.size() ? other.digits[j] : 0) 
                          + carry;
            result.digits[i + j] = cur % BASE;
            carry = cur / BASE;
        }
    }
    
    // 去除前导零
    while (result.digits.size() > 1 && result.digits.back() == 0)
        result.digits.pop_back();
    
    return result;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

void fft(vector<Complex>& a, int limit, int type) {
    // 位反转置换
    for (int i = 0; i < limit; i++) {
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    
    // 蝴蝶操作
    for (int mid = 1; mid < limit; mid <<= 1) {
        Complex wn(cos(PI / mid), type * sin(PI / mid));
        for (int left = 0; left < limit; left += (mid << 1)) {
            Complex w(1, 0);
            for (int k = 0; k < mid; k++) {
                Complex x = a[left + k];
                Complex y = w * a[left + k + mid];
                a[left + k] = x + y;
                a[left + k + mid] = x - y;
                w *= wn;
            }
        }
    }
    
    if (type == -1) {
        for (int i = 0; i < limit; i++) {
            a[i] /= limit;
        }
    }
}

1
2
3
4

for (int i = 1; i <= n / 2; i++) {
    result *= i * (n - i + 1); 
}
if (n % 2 == 1) result *= (n + 1) / 2;  // 奇数时补中间的项

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

void factorial_binary_splitting(mpz_t result, unsigned long n, unsigned long m) {
    if (n - m == 0) {
        mpz_set_ui(result, 1);
        return;
    }
    if (n - m == 1) {
        mpz_set_ui(result, n);
        return;
    }
    
    unsigned long mid = m + (n - m) / 2;
    mpz_t high, low;
    mpz_init(high);
    mpz_init(low);
    
    factorial_binary_splitting(high, n, mid);
    factorial_binary_splitting(low, mid, m);
    
    mpz_mul(result, high, low);
    
    mpz_clear(high);
    mpz_clear(low);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

class FactorialPrimeSwingOptimized {
private:
    std::vector<unsigned long> primes;
    
    void swing(mpz_t result, unsigned long n) {
        mpz_set_ui(result, 1);
        for (unsigned long p : primes) {
            if (p > n) break;
            
            long exponent = 0;
            unsigned long temp = n;
            while (temp > 0) {
                temp /= p;
                exponent += temp;
            }
            temp = n / 2;
            while (temp > 0) {
                temp /= p;
                exponent -= 2 * temp;
            }
            
            if (exponent > 0) {
                mpz_t pow_val;
                mpz_init(pow_val);
                mpz_ui_pow_ui(pow_val, p, exponent);
                mpz_mul(result, result, pow_val);
                mpz_clear(pow_val);
            }
        }
    }
    
    void rec_factorial(mpz_t result, unsigned long n) {
        if (n < 2) {
            mpz_set_ui(result, 1);
            return;
        }
        
        mpz_t rec_half;
        mpz_init(rec_half);
        rec_factorial(rec_half, n / 2);
        
        mpz_t swing_val;
        mpz_init(swing_val);
        swing(swing_val, n);
        
        mpz_mul(result, rec_half, rec_half);
        mpz_mul(result, result, swing_val);
        
        mpz_clear(rec_half);
        mpz_clear(swing_val);
    }
    
public:
    FactorialPrimeSwingOptimized(unsigned long max_n) {
        primes = sieve_primes(max_n);
    }
    
    void factorial(mpz_t result, unsigned long n) {
        if (n < 20) {
            mpz_set_ui(result, 1);
            for (unsigned long i = 2; i <= n; i++) {
                mpz_mul_ui(result, result, i);
            }
            return;
        }
        
        rec_factorial(result, n);
        unsigned long shift = n - __builtin_popcountll(n);
        mpz_mul_2exp(result, result, shift);
    }
};

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

mpz_mul(n=10000)
└─ mpn_mul(an=10000, bn=10000)
   ├─ 若 an≈bn ≥ FFT_THRESHOLD → mpn_fft_mul
   │   ├─ 位打包
   │   ├─ 三层 FFT（蝴蝶）
   │   └─ 分量乘 → 递归 mpn_mul（小尺寸）
   └─ 否则 → mpn_toom44_mul
       ├─ 5 段插值
       ├─ 9 次点乘 → 递归 mpn_mul（更小）
       └─ 插值逆变换

1
2
3
4
5
6

if (n ≤ MUL_THRESHOLD)          // 很小就直接连乘
    return basecase_mul(vec, n);
mid = n >> 1;
left  = fac_prod(vec, mid);
right = fac_prod(vec + mid, n - mid);
return left * right;            // 大数乘法，自动选 Karatsuba/Toom/FFT

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

#include <gmp.h>

int main() {
    mpz_t result;
    mpz_init(result);
    
    // 计算10000!
    mpz_fac_ui(result, 10000);
    
    gmp_printf("10000! = %Zd\n", result);
    
    mpz_clear(result);
    return 0;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

#include <omp.h>

void factorize_factorial_multithread(unsigned int n,
                                     const std::vector<unsigned int>& primes,
                                     std::vector<unsigned int>& exponents) {
    exponents.resize(primes.size(), 0);
    
    #pragma omp parallel for schedule(dynamic)
    for (size_t i = 0; i < primes.size(); i++) {
        unsigned int exp = 0;
        unsigned long long power = primes[i];
        while (power <= n) {
            exp += n / power;
            if (power > n / primes[i]) break;
            power *= primes[i];
        }
        exponents[i] = exp;
    }
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

// GPU核函数：计算质数的幂
template <class params>
__global__ void kernel_factorial(cgbn_error_report_t *report,
                                 instance_t<params> *instances,
                                 uint32_t count) {
    int32_t instance = (blockIdx.x * blockDim.x + threadIdx.x) / params::TPI;
    if (instance >= count) return;
    
    factorial_gpu_t<params> factorial(cgbn_report_monitor, report, instance);
    factorial(instances);
}

// 主机端：质因数分解和结果合并
void calculate_factorial_gpu(unsigned int n, mpz_t result) {
    // 1. CPU生成质数表
    std::vector<unsigned int> primes = generate_primes(n);
    
    // 2. CPU质因数分解（多线程）
    std::vector<unsigned int> exponents;
    factorize_factorial_multithread(n, primes, exponents);
    
    // 3. GPU并行计算质数幂
    process_batch_gpu(primes, exponents, result);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

#include <cufft.h>
#include <cuda_runtime.h>
#include <thrust/device_vector.h>

// 基于 cuFFT 的大数乘法核心流程
void big_integer_multiply_cufft(
    const std::vector<double>& A,   // 被乘数的多项式系数（以 B 为基）
    const std::vector<double>& B,   // 乘数的多项式系数
    std::vector<double>& C)         // 乘积的多项式系数
{
    size_t lenA = A.size();
    size_t lenB = B.size();
    size_t conv_len = lenA + lenB - 1;  // 线性卷积长度

    // 1. 零填充至 2 的幂，避免循环卷积混叠
    size_t fft_len = 1;
    while (fft_len < conv_len) fft_len <<= 1;

    // 2. 将输入数据拷贝至 GPU
    thrust::device_vector<cufftDoubleComplex> d_A(fft_len, make_cuDoubleComplex(0.0, 0.0));
    thrust::device_vector<cufftDoubleComplex> d_B(fft_len, make_cuDoubleComplex(0.0, 0.0));
    for (size_t i = 0; i < lenA; i++) d_A[i] = make_cuDoubleComplex(A[i], 0.0);
    for (size_t i = 0; i < lenB; i++) d_B[i] = make_cuDoubleComplex(B[i], 0.0);

    // 3. 创建 cuFFT plan 并执行正向 FFT
    cufftHandle plan;
    cufftPlan1d(&plan, fft_len, CUFFT_Z2Z, 1);  // Z2Z: 双精度复数到复数
    cufftExecZ2Z(plan, 
        thrust::raw_pointer_cast(d_A.data()),
        thrust::raw_pointer_cast(d_A.data()), 
        CUFFT_FORWARD);
    cufftExecZ2Z(plan,
        thrust::raw_pointer_cast(d_B.data()),
        thrust::raw_pointer_cast(d_B.data()), 
        CUFFT_FORWARD);

    // 4. 频域逐点相乘
    thrust::transform(d_A.begin(), d_A.end(), d_B.begin(), d_A.begin(),
        [] __device__ (cufftDoubleComplex x, cufftDoubleComplex y) {
            return make_cuDoubleComplex(
                x.x * y.x - x.y * y.y,
                x.x * y.y + x.y * y.x);
        });

    // 5. 逆 FFT 并缩放
    cufftExecZ2Z(plan,
        thrust::raw_pointer_cast(d_A.data()),
        thrust::raw_pointer_cast(d_A.data()),
        CUFFT_INVERSE);
    cufftDestroy(plan);

    // 6. 拷贝回主机端并执行进位归一化
    C.resize(conv_len);
    for (size_t i = 0; i < conv_len; i++)
        C[i] = d_A[i].x / fft_len;  // IFFT 的缩放因子
}

1
2
3
4
5
6
7

void factorial_gpu_fft(mpz_t result, unsigned long n) {
    // 1. CPU 端完成 PrimeSwing 分解，生成 Swing 数序列
    // 2. 将大数乘法密集的子任务打包发送到 GPU
    // 3. GPU 端调用 cuFFT 完成频域卷积
    // 4. 主机端回收结果并执行进位归一化
    // 5. 按 PrimeSwing 公式组装最终阶乘
}