C++实现Prime-k因子问题结果异常求助

核心问题：阶乘溢出导致计算错误 #

你的代码中直接计算factorial(n-i)，当素数p较大时（比如p≥23），(p-1)!的数值会超过unsigned long long的最大值（约1.8×10^19），发生溢出后数值失真，最终导致c%n的结果完全错误。这是你得到错误输出的主要原因。

另外，你的素数生成函数Prime_lst从i=5开始，漏掉了素数2和3，如果题目要求计算所有素数的结果，这部分也会导致总和偏差。

解决方案：利用数论性质避免直接计算大阶乘 #

根据威尔逊定理，对于素数p，有(p-1)! ≡ -1 mod p。结合模逆元的性质，我们可以递推计算每个阶乘项的模p值，所有运算都在模p范围内进行，不会发生溢出：

1. 对于素数p：

   • (p-1)! ≡ -1 mod p → 记为term1 = p-1
   • (p-2)! = (p-1)! × inv(p-1) mod p，其中inv(x)是x在模p下的逆元。由于p-1 ≡ -1 mod p，inv(-1) ≡ -1 mod p，所以term2 = 1
   • 后续项(p-k)!可以通过前一项乘以inv(p-(k-1))递推得到，逆元可以用费马小定理计算：inv(x) = x^(p-2) mod p（因为p是素数）

2. 修正素数生成逻辑，包含2和3，并优化遍历效率。

修正后的代码 #

#include <iostream>    
#include <vector> 

using namespace std;

bool prime(unsigned long long n) {
    if (n == 2 || n == 3) {
        return true; 
    }
    if (n % 2 == 0 || n % 3 == 0) return false;
    for (unsigned long long i = 5; i * i <= n; i += 6) {
        if (n % i == 0 || n % (i + 2) == 0) return false;
    }
    return true;
}

vector<unsigned long long> Prime_lst(unsigned long long n) {
    vector<unsigned long long> v;
    if (n >= 2) v.push_back(2);
    if (n >= 3) v.push_back(3);
    for (unsigned long long i = 5; i <= n; i += 6) {
        if (prime(i)) v.push_back(i);
        if (i + 2 <= n && prime(i + 2)) v.push_back(i + 2);
    }
    return v;
}

unsigned long long mod_pow(unsigned long long base, unsigned long long exp, unsigned long long mod) {
    unsigned long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp % 2 == 1) {
            result = (result * base) % mod;
        }
        base = (base * base) % mod;
        exp /= 2;
    }
    return result;
}

unsigned long long mod_inv(unsigned long long x, unsigned long long mod) {
    return mod_pow(x, mod - 2, mod);
}

unsigned long long prime_p(unsigned long long p) {
    if (p == 2) {
        return 1 % 2; // 仅(p-1)! = 1!
    }
    if (p == 3) {
        return (2 + 1) % 3; // (3-1)! + (3-2)! = 2! +1!
    }
    if (p == 5) {
        return (24 + 6 + 2 + 1 + 1) % 5; // 4!+3!+2!+1!+0!，0! =1
    }
    // 处理p>5的素数
    unsigned long long term1 = p - 1; // (p-1)! ≡ -1 mod p
    unsigned long long term2 = (term1 * mod_inv(p-1, p)) % p;
    unsigned long long term3 = (term2 * mod_inv(p-2, p)) % p;
    unsigned long long term4 = (term3 * mod_inv(p-3, p)) % p;
    unsigned long long term5 = (term4 * mod_inv(p-4, p)) % p;
    
    return (term1 + term2 + term3 + term4 + term5) % p;
}

unsigned long long sum_prime_p(unsigned long long n) {
    unsigned long long total = 0; 
    vector<unsigned long long> prime_lst = Prime_lst(n); 
    for (unsigned long long p : prime_lst) {
        total += prime_p(p); 
    }
    return total; 
}

int main() {
    cout << sum_prime_p(100) << endl; 
    return 0;
}

额外优化建议（针对1e8范围） #

当前的试除法素数生成在1e8范围内会非常慢，建议替换为埃拉托斯特尼筛法（埃氏筛），可以大幅提升素数生成效率。例如：

vector<bool> sieve(unsigned long long n) {
    vector<bool> is_prime(n+1, true);
    is_prime[0] = is_prime[1] = false;
    for (unsigned long long i = 2; i*i <= n; ++i) {
        if (is_prime[i]) {
            for (unsigned long long j = i*i; j <= n; j += i) {
                is_prime[j] = false;
            }
        }
    }
    return is_prime;
}

vector<unsigned long long> Prime_lst(unsigned long long n) {
    vector<unsigned long long> v;
    vector<bool> is_prime = sieve(n);
    for (unsigned long long i = 2; i <= n; ++i) {
        if (is_prime[i]) {
            v.push_back(i);
        }
    }
    return v;
}

这样处理1e8范围的素数生成会高效很多。

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