P1829 [国家集训队]Crash的数字表格 / JZPTAB

Posted 2021-01-25 dreagonm

推式子太快乐啦！虽然我好蠢而且dummy和maomao好巨（划掉）

思路

莫比乌斯反演的题目
首先这题有(O(sqrt n))的做法但是我没写咕咕咕
然后就是爆推一波式子
[ sum_{i=1}^{n}sum_{j=1}^{m}lcm(i,j) ]

[ sum_{i=1}^{n}sum_{j=1}^{m}frac{i imes j}{gcd(i,j)} ]

设$ gcd(i,j)=d$
[ sum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}i imes j imes[gcd(i,j)=1] ]

[ p_1=lfloorfrac{n}{d} floor\ p_2=lfloorfrac{m}{d} floor ]

设(f(x))为满足条件(1 le i le p_1)且$1 le j le p_2 (且)[gcd(i,j)=x](的)i imes j$的和

则可以推出(F(x))为
[ F(x)=sum_{x|d}f(d) ]
所以(F(x))为满足条件(1 le i le p_1)且$1 le j le p_2 (且)[x|gcd(i,j)](的)i imes j$的和

所以
[ F(x)=x^2 imes sum_{i=1}^{lfloorfrac{p_1}{x} floor}i imes sum_{j=1}^{lfloorfrac{p_2}{x} floor}j ]
因为
[ f(x)=sum_{x|d}mu(frac{d}{x})F(d) ]
所以
[ f(x)=sum_{x|d}mu(frac{d}{x}) imes d^2 imes sum_{i=1}^{lfloorfrac{p_1}{d} floor}i imes sum_{j=1}^{lfloorfrac{p_2}{d} floor}j ]
所以
[ f(x)=sum_{t=1}mu(t) imes (tx)^2 imes sum_{i=1}^{lfloorfrac{p_1}{tx} floor}i imes sum_{j=1}^{lfloorfrac{p_2}{tx} floor}j ]

[ ans=sum_{d=1}^{n}dsum_{t=1}^{lfloorfrac{n}{d} floor}mu(t) imes (t)^2 imes sum_{i=1}^{lfloorfrac{n}{td} floor}i imes sum_{j=1}^{lfloorfrac{m}{td} floor}j ]

然后两个整除分块搞定，复杂度(O(n))

然后dummy教了我一种新的计算方式
[ sum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}i imes j imes[gcd(i,j)=1] ]
发现[gcd(i,j)=1]的形式有点眼熟
直接把(mu)带进去算
[ sum_{d=1}^{n}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}i imes j imes sum_{p|gcd(i,j)}mu(p) ]
算得
[ sum_{d=1}^ndsum_{p=1}^{lfloorfrac{n}{d} floor}p^2mu(p)sum_{i=1}^{lfloorfrac{n}{dp} floor}isum_{j=1}^{lfloorfrac{n}{dp} floor}j ]

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;
const int INV = 10050505;
const int MOD = 20101009;
int isprime[10010000],iprime[10010000],cnt,mu[10010000],summu[10010000],n,m;
void prime(int n){
    isprime[n]=true;
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i])
            iprime[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&iprime[j]*i<=n;j++){
            isprime[iprime[j]*i]=true;
            mu[iprime[j]*i]=-mu[i];
            if(i%iprime[j]==0){
                mu[iprime[j]*i]=0;
                break;
            }
        }
    }
    for(int i=1;i<=n;i++)
        summu[i]=(summu[i-1]%MOD+1LL*(mu[i]%MOD+MOD)%MOD*i%MOD*i%MOD)%MOD;
}
long long calc2(int n,int m){
    long long ans=0;
    for(int l=1,r;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=(ans%MOD+1LL*(summu[r]-summu[l-1]%MOD+MOD)%MOD*(1+n/l)%MOD*(n/l)%MOD*INV%MOD*(1+m/l)%MOD*(m/l)%MOD*INV%MOD)%MOD;
    }
    return ans;
}
long long calc1(int n,int m){
    long long ans=0;
    for(int l=1,r;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=(ans%MOD+1LL*(r-l+1)%MOD*(r+l)%MOD*INV%MOD*calc2(n/l,m/l)%MOD)%MOD;
    }
    return ans;
}
int main(){
    prime(10001000);
    scanf("%d %d",&n,&m);
    if(n<m)
        swap(n,m);
    long long ans=calc1(n,m);
    printf("%lld",ans);
}