洛谷1829：crash的数字表格

洛谷1829：crash的数字表格

题意：

• 求(sum_{i=1}^Nsum_{j=1}^Mlcm(i,j))。
• 数据范围(:n,mleq10^7)。

思路：

易得：原式

[sum_{i=1}^Nsum_{j=1}^Mfrac{ij}{gcd(i,j)}​ ]

枚举(gcd(i,j)=d)。

[sum_{d=1}^{min(N,M)}sum_{i=1}^Nsum_{j=1}^M[gcd(i,j)=d]frac{ij}{d} ]

把(d)给提出来，其实就是将枚举项(i,j)看成是(di,dj)。

[sum_{d=1}^{min(N,M)}frac{1}{d}sum_{i=1}^{frac{N}{d}}sum_{j=1}^{frac{M}{d}}[gcd(i,j)=1]idjd.\sum_{d=1}^{min(N,M)}dsum_{i=1}^{frac{N}{d}}sum_{j=1}^{frac{M}{d}}[gcd(i,j)=1]ij ]

我们知道莫比乌斯函数的一个重要性质。

[sum_{d|n}mu(d)=[n=1] ]

代入可得：

[sum_{d=1}^{min(N,M)}dsum_{i=1}^{frac{N}{d}}sum_{j=1}^{frac{M}{d}}sum_{x|gcd(i,j)}mu(x)ij ]

改为枚举(x)。

[sum_{d=1}^{min(N,M)}dsum_{x=1}^{min(frac{N}{d},frac{M}{d})}mu(x)sum_{i=1}^{frac{N}{d}}sum_{j=1}^{frac{M}{d}}ij[x|gcd(i,j)] ]

因为后面有(x|gcd(i,j))的条件，所以我们可以将枚举(i,j)改为枚举(ix,jx)。

那么就有:

[sum_{d=1}^{min(N,M)}dsum_{x=1}^{min(frac{N}{d},frac{M}{d})}mu(x)sum_{i=1}^{frac{N}{dx}}sum_{j=1}^{frac{M}{dx}}ijx^2 ]

接着把(x^2)提到前面去：

[sum_{d=1}^{min(N,M)}dsum_{x=1}^{min(frac{N}{d},frac{M}{d})}x^2mu(x)sum_{i=1}^{frac{N}{dx}}isum_{j=1}^{frac{M}{dx}}j ]

后面的两个(sum)其实就是两个等差数列求和。

中间的(sum)直接(O(n))预处理即可。

我们设(sum(frac{N}{d},frac{M}{d})=sum_{x=1}^{min(frac{N}{d},frac{M}{d})}x^2mu(x)sum_{i=1}^{frac{N}{dx}}isum_{j=1}^{frac{M}{dx}}j).

那么原式其实就是(sum d*sum(frac{N}{d},frac{M}{d})).

两次数论分块。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e7+10;
const int mod = 20101009;
const int inv2 = (mod+1)/2;
ll n, m;

int primes[maxn], cnt;
ll mu[maxn];
bool vis[maxn];
void init(int n)
{
    mu[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!vis[i])
        {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1; primes[j] <= n/i; j++)
        {
            vis[primes[j]*i] = 1;
            if(i % primes[j] == 0) break;
            else mu[i*primes[j]] = -mu[i];
        }
    }
    for(ll i = 1; i <= n; i++)
    mu[i] = (mu[i-1]+(mu[i]+mod)*i%mod*i%mod)%mod;
}

ll sum(ll x, ll y){
    return x%mod*(x+1)%mod*inv2%mod*y%mod*(y+1)%mod*inv2%mod;
}

ll work(ll x, ll y)
{
    ll res = 0;
    for(ll l = 1, r; l <= min(x, y); l = r+1)
    {
        r = min(x/(x/l), y/(y/l));
        res = (res+((mu[r]-mu[l-1])%mod+mod)%mod*sum(x/l, y/l)%mod)%mod;
    }
    return res;
}

void solve(ll n, ll m)
{
    ll res = 0;
    for(ll l = 1, r; l <= min(n, m); l=r+1)
    {
        r = min(n/(n/l), m/(m/l));
        res = (res+(l+r)%mod*(r-l+1)%mod*inv2%mod*work(n/l, m/l)%mod)%mod;
    }
    cout << res << endl;
}

int main()
{
    init(maxn-5);
    cin >> n >> m;
    solve(n, m);
    return 0;
}