nju1017 Heresy

原题网址：http://acm.nju.edu.cn:8000/problem.php?id=1017
原题大意：求$\sum_{i=1}^{n}\sum_{j=1}^{m}i^{2}j^{2}gcd(i,j)$
和BZOJ2154这类题目相似，都是莫比乌斯反演。
记$F(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}i^{2}j^{2}[gcd(i,j)=1]$
则由$\sum_{d|n}\mu(d)=[n=1]$可知
$F(n,m)=\sum_{i=1}^{n}i^{2}\sum_{j=1}^{m}j^{2}\sum_{d|gcd(i,j)}\mu(d)\\=\sum\mu(d)\sum_{1\leq i\leq n\ and\ d|i}i^{2}\sum_{1\leq j\leq m\ and\ d|j}j^{2}\\=\sum\mu(d)d^{4}\sum_{i=1}^{[\frac{n}{d}]}i^{2}\sum_{j=1}^{[\frac{m}{d}]}j^{2}\\=\sum\mu(d)d^{4}Sum([\frac{n}{d}],[\frac{m}{d}])$
其中$Sum(i,j)=\frac{1}{6}i(i+1)(2i+1)*\frac{1}{6}j(j+1)(2j+1)$
我们枚举$gcd(i,j)=d$，则有
$ans=\sum_{d=1}^{min(n,m)}d^{5}F([\frac{n}{d}],[\frac{m}{d}])\\=\sum_{d=1}^{min(n,m)}d^{5}\sum_{i=1}^{min([\frac{n}{d}],[\frac{m}{d}])}\mu(i)i^{4}Sum([\frac{n}{di}],[\frac{m}{di}])$

5次方的数据做的时候总是出问题，不能直接用公式求和，最后只能设数组记录数值才过的

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>

using namespace std;
const int N = 1000010;
#define ll long long
const int MOD = 1000000007;
bool check[N + 2];
int prime[N + 2];
int mu[N + 2];
ll s[N + 2];
ll ssum[N + 2];

void Moblus(){
    memset(check, false, sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for (int i = 2; i <= N; i++){
        if (!check[i]){
            prime[tot++] = i;
            mu[i] = -1;
        }
        for (int j = 0; j<tot; j++){
            if (i*prime[j]>N) break;
            check[i*prime[j]] = true;
            if (i%prime[j] == 0){
                mu[i*prime[j]] = 0;
                break;
            }
            else{
                mu[i*prime[j]] = -mu[i];
            }
        }
    }
    s[0] = 0;
    ssum[0] = 0;
    for (ll i = 1; i <= N; i++){
        s[i] = ((s[i - 1] + (1ll * mu[i] * i%MOD*i%MOD*i%MOD*i%MOD)) % MOD + MOD) % MOD;
        ssum[i] = (ssum[i-1] + i*i%MOD*i%MOD*i%MOD*i%MOD)%MOD;
    }
}

ll sum(ll a, ll b){
    ll x = (a*(a + 1) / 2 * (2 * a + 1) / 3) % MOD;
    ll y = (b*(b + 1) / 2 * (2 * b + 1) / 3) % MOD;
    return (x*y%MOD);

}

ll calf(ll x, ll y){
    ll temp = 0;
    ll pos = 0;
    for (ll i = 1; i <= x; i = pos + 1){
        pos = min(x / (x / i), y / (y / i));
        temp = (temp + ((s[pos] - s[i - 1]) % MOD + MOD) % MOD*sum((ll) x / i, (ll)y / i) % MOD) % MOD;
    }
    return temp;
}

int main(){
    int T;
    scanf("%d", &T);
    int cas = 1;
    Moblus();
    while (T--){
        int n, m;
        scanf("%d%d", &m, &n);
        if (m>n) swap(m, n);
        ll ans = 0;
        ll pos = 0;
        for (ll d = 1; d <= m; d=pos+1){
            pos = min(m / (m / d), n / (n / d));
            ans = (ans + 1ll * ((ssum[pos] - ssum[d-1]) % MOD*calf(m / d, n / d) % MOD + MOD) % MOD) % MOD;
        }
        printf("Case #%d: %lld\n", cas++, ans);
    }
    //system("pause");
    return 0;
}