bzoj2154 Crash的数字表格

题目链接：bzoj2154
题目大意：
给定n、m，求$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i*j)$。
1≤n，m≤10000000

题解：
莫比乌斯反演
把原式化开

$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i*j)=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i*j}{gcd(i,j)}$$
枚举约数(省略前面的了)
$$\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d|i,d|j}\frac{i*j}{d}[gcd(i,j)=d]$$
交换枚举倍数和约数
$$\sum_{d=1}^{min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\frac{id*jd}{d}[gcd(i,j)=1]$$
$[gcd(i,j)=1]就是\epsilon((i,j))$，套 $1*\mu=\epsilon$再约下分什么的就变成了
$$\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}{i*j}\sum_{t|i,t|j}\mu(t)$$
再次交换枚举倍数和约数
$$\sum_{d=1}^{min(n,m)}d\sum_{t=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\sum_{i=1}^{\lfloor\frac{n}{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dt}\rfloor}{it*jt}*\mu(t)$$
$$\sum_{d=1}^{min(n,m)}d\sum_{t=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}t^2*\mu(t)\sum_{i=1}^{\lfloor\frac{n}{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dt}\rfloor}{i*j}$$
设 $sum(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}{i*j}=\sum_{i=1}^{n}i \sum_{j=1}^{m}j$这个直接上求和公式
则上式变成 $$\sum_{d=1}^{min(n,m)}d\sum_{t=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}t^2*\mu(t)*sum({\lfloor\frac{n}{dt}\rfloor},{\lfloor\frac{m}{dt}\rfloor})$$
设 $f(n,m)=\sum_{t=1}^{min(n,m)}t^2*\mu(t)*sum({\lfloor\frac{n}{t}\rfloor},{\lfloor\frac{m}{t}\rfloor})$这个分块求。
于是上式又变成了 $$\sum_{d=1}^{min(n,m)}d\times f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)$$
这个也可以分块了。
所以时间复杂度大概是 $O(n)$..吧
LL、mod什么的真的要注意啊不然WA到死。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
#define N 10000100

const LL ny=10050505;
const LL mod=20101009;
int cnt,pri[N/10];bool ispri[N];LL mu[N],pre[N];
int mymin(int x,int y){return (x<y)?x:y;}
void mobius(int lim)
{
    cnt=0;mu[1]=1;
    for (int i=2;i<=lim;i++)
    {
        if (!ispri[i]) {pri[++cnt]=i;mu[i]=-1;}
        for (int j=1;j<=cnt && i*pri[j]<=lim;j++)
        {
            ispri[i*pri[j]]=true;
            if (i%pri[j]==0)
            {
                mu[i*pri[j]]=0;
                break;
            }
            mu[i*pri[j]]=-mu[i];
        }
    }
}
LL sum(int n,int m)
{
    LL s1=(LL)(n)*(LL)(n+1)%mod*ny%mod;
    LL s2=(LL)(m)*(LL)(m+1)%mod*ny%mod;
    return s1*s2%mod;
}
LL f(int x,int y)
{
    int r,lim=mymin(x,y);LL ans=0;
    for (int i=1;i<=lim;i=r+1)
    {
        r=mymin(x/(x/i),y/(y/i));
        ans=(ans+sum(x/i,y/i)*(pre[r]-pre[i-1])%mod)%mod;
    }while (ans<0) ans+=mod;
    return ans;
}
int main()
{
    //freopen("a.in","r",stdin);
    //freopen("a.out","w",stdout);
    int n,m,r,i,lim;LL ans=0;
    scanf("%d%d",&n,&m);
    lim=mymin(n,m);mobius(lim);
    pre[1]=mu[1];for (i=2;i<=lim;i++) pre[i]=(pre[i-1]+(LL)i*i%mod*mu[i]%mod)%mod;
    for (i=1;i<=lim;i=r+1)
    {
        r=mymin(n/(n/i),m/(m/i));
        ans=(ans+(LL)(r-i+1)*(r+i)%mod*ny%mod*f(n/i,m/i)%mod)%mod;
    }while (ans<0) ans+=mod;
    printf("%lld\n",ans);
    return 0;
}