链接:https://ac.nowcoder.com/acm/contest/549/J
来源:牛客网

小A最近开始研究数论题了,这一次他随手写出来一个式子,
∑i=1n∑j=1mgcd(i,j)2,但是他发现他并不太会计算这个式子,你可以告诉他这个结果吗,答案可能会比较大,请模上1000000007。


∑i=1n∑j=1mgcd(i,j)2 
∑d=1min(n,m)d2(∑i=1n∑j=1mgcd(i,j)=d)(d)=(∑i=1n∑j=1mgcd(i,j)=d)F(t)=∑d|tf(d)=[nd][md]f(d)=∑d|(td)F(t)=∑d|(td)[nd][md]
处理因子直接计算复杂度
O(nlogn)就可以通过本题
莫比乌斯反演
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e6;
ll prime[maxn + 5];
ll isprime[maxn + 5];
ll mod = 1e9+7;
ll mbs[maxn + 5];
int cnt;
ll qpow(ll a, ll n, ll mod) {
    ll ans = 1;
    a = a % mod;
    while(n) {
        if(n&1) {
            ans = (ans * a) % mod;
        }
        a = (a * a) % mod;
        n >>= 1;
    }
    return ans;
}
 
void get_mbs(int n) {
    mbs[1] = 1;
    isprime[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!isprime[i]) {
            prime[cnt++] = i;
            mbs[i] = -1;
        }
        for(int j = 0; j < cnt && i * prime[j] <= n; j++) {
            isprime[i * prime[j]] = 1;
            if(i % prime[j] == 0) {
                break;
            }
            mbs[i * prime[j]] = -mbs[i];
        }
        mbs[i]+= mbs[i - 1];
    }  
}
int main() {
    ll n, m, k;
    scanf("%lld%lld", &n, &m);
    k = 2;
    get_mbs(maxn);
    if(n > m) {
        swap(n, m);
    }
    ll ans = 0;
    for(ll d = 1; d <= n; d++) {
        ll a = n / d;
        ll b = m / d;
        ll res = 0;
        for(ll i = 1, j = 1; i <= a; i = j + 1) {
            j = min(a / (a / i), b / (b / i));
            ll y = ((a / i) * (b / i)) % mod;
            res = (res + y * (mbs[j] - mbs[i - 1])) % mod;
        }
        ll x = qpow(d, k, mod);
        ans = (ans + x * res) % mod;
    }
    printf("%lld\n", ans);
    return 0;
}