HDU 4059 The Boss on Mars - 容斥原理計算四次方和

問題描述

定義 S = 1⁴ + 2⁴ + ... + n⁴。從中減去與 n 互質的數的四次方和。求減去的總和模 1,000,000,007。n ≤ 10¹²，多測試案例。

數學轉化

求 T = Σ_{1≤i≤n, gcd(i,n)>1} i⁴。
等價於 T = S - Σ_{gcd(i,n)=1} i⁴。
用容斥原理直接計算 T：設 n 的質因子為 p₁, p₂, ..., p_k。不互質 = 至少被一個質因子整除。
T = Σ_{非空子集 J} (-1)^{|J|+1} * Σ_{i≤n, ∏_{j∈J} p_j | i} i⁴

四次方和公式

Σ_{i=1}^m i⁴ = m(m+1)(2m+1)(3m²+3m-1) / 30。

推導：(m+1)⁵ = m⁵ + 5m⁴ + 10m³ + 10m² + 5m + 1，求和望遠鏡消項，代入平方和、立方和公式化簡得出。

模逆元處理除法

模 MOD = 1,000,000,007 (質數)。30 的逆元 inv30 = 30^{MOD-2} mod MOD (費馬小定理) 或擴展歐幾里得。預計算 inv30 = 233333335。

算法流程

1. 質因子分解 n (試除至 √n，n≤10¹² 最大質因子 ≤ 10⁶，預篩質數表)。
2. 枚舉質因子非空子集 (最多 2^k - 1，k ≤ 11 因 2·3·5·...·31 > 10¹²)。
3. 對每個子集，計算乘積 d，項數 m = n/d，四次方和 sum4(m)，貢獻 d⁴ * sum4(m)。
4. 奇數項加，偶數項減 (容斥系數)。
5. 所有運算模 MOD，乘法用 int128 或分步取模防溢出。

參考代碼框架

#include 
using namespace std;
typedef long long ll;
const ll MOD = 1000000007;
const ll INV30 = 233333335; // 30^{-1} mod MOD

ll mod_pow(ll a, ll e) { ll r=1; for(;e;e>>=1,a=a*a%MOD) if(e&1) r=r*a%MOD; return r; }

ll sum4(ll m) {
    m %= MOD;
    ll t1 = m * (m+1) % MOD;
    ll t2 = (2*m+1) % MOD;
    ll t3 = (3*m%MOD*m%MOD + 3*m - 1) % MOD;
    return t1 * t2 % MOD * t3 % MOD * INV30 % MOD;
}

ll mul_mod(ll a, ll b) { return (__int128)a * b % MOD; } // 需支持 int128

int main() {
    // 篩質數至 1e6
    vector primes;
    vector is_composite(1000001);
    for (int i=2;i<=1000000;i++) if(!is_composite[i]) {
        primes.push_back(i);
        if((ll)i*i<=1000000) for(int j=i*i;j<=1000000;j+=i) is_composite[j]=1;
    }

    int T; scanf("%d", &T);
    while (T--) {
        ll n; scanf("%lld", &n);
        // 質因子分解
        vector pf;
        ll tmp = n;
        for (int p : primes) {
            if ((ll)p*p > tmp) break;
            if (tmp % p == 0) {
                pf.push_back(p);
                while (tmp % p == 0) tmp /= p;
            }
        }
        if (tmp > 1) pf.push_back(tmp);

        ll ans = 0;
        int k = pf.size();
        for (int mask = 1; mask < (1<>i & 1) {
                d = mul_mod(d, pf[i]);
                bits++;
            }
            ll m = n / d;
            ll term = mul_mod(mul_mod(d, d), mul_mod(d, d)); // d^4
            term = mul_mod(term, sum4(m));
            if (bits & 1) ans = (ans + term) % MOD;
            else ans = (ans - term + MOD) % MOD;
        }
        printf("%lld\n", ans);
    }
    return 0;
}

關鍵點

• n ≤ 10¹²，質因子數極少 (≤ 11)，子集枚舉極快 (≤ 2047)。
• 四次方和公式除以 30 → 模逆元。
• 大數乘法必須用 __int128 或分步取模，防 64 位溢出。
• 質因子分解用預篩質數表至 10⁶，試除分解 10¹² 極快。