数论--莫比乌斯反演&分块--bzoj2154 Crash的数字表格

大佬博客

http://www.cnblogs.com/cjyyb/p/8253033.html

//1.mu[1] = 1
//2.long long 注意
//3.2次分块好神奇，而且min(n/ (n/ i),m/ (m/ i))可以解决2个变量一起分块
//4.直接预处理1e7超时，max(n,m) 7s+
#include <cstdio>
#include <iostream>
#include <set>
using namespace std;
const int mod = 20101009;
const int maxn = 1e7 + 5;

typedef long long ll;
ll n,m;
int mu[maxn];
bool isprime[maxn];
int prime[maxn / 10],sz = 0;
ll sum[maxn];

void init(int maxn)
{
    for(int i = 2;i < maxn;i ++){
        if(!isprime[i]) prime[sz ++] = i,mu[i] = -1;
        for(int j = 0;j < sz && 1ll * i * prime[j] < maxn;j ++){
            isprime[1ll * i * prime[j]] = 1;
            if(i % prime[j] == 0){
                mu[1ll * i * prime[j]] = 0;
                break;
            }
            mu[1ll * i * prime[j]] = -mu[i];
        }
    }
    mu[1] = 1;//！！！
    for (int i = 1; i < maxn; i ++) {
        sum[i] = sum[i - 1] + 1ll * mu[i] * i * i % mod;
        sum[i] %= mod;
    }
}

//sum{i=1 to x,sum{j = 1 to y,i * j} }
ll f(ll x,ll y)
{
    x = (1 + x) * x / 2 % mod;
    y = (1 + y) * y / 2 % mod;
    return x * y % mod;
}
ll g(ll x,ll y)
{
    ll t = 0;
    int l,r;
    for (l = 1; l <= x; l = r + 1) {
        r = min(x / (x / l),y / (y / l));
        t += (sum[r] - sum[l - 1]) * f(x / l,y / l) % mod;
        t %= mod;
    }
    return t;
}
int main()
{
    scanf("%lld%lld",&n,&m);
    init(max(n,m) + 3);
    if(n > m) swap(n, m);
    ll ans = 0;
    int l,r;
    for (l = 1; l <= n; l = r + 1) {
        r = min(n / (n / l),m / (m / l));
        ans += 1ll * (l + r) * (r - l + 1) / 2 % mod * g(n/ l,m / l) % mod;
        ans %= mod;
    }
    ans = (ans % mod + mod) % mod;
    printf("%lld\n",ans);
    return 0;
}