BZOJ2154: Crash的数字表格 & BZOJ2693: jzptab

【传送门：BZOJ2154&BZOJ2693】

---

简要题意：

　　给出n,m，求$sum_{i=1}^{n}sum_{j=1}^{m}LCM(i,j)$

---

题解：

　　莫比乌斯反演（因为BZOJ2693是多组数据，数据强一点，所以代码用BZOJ2693的）

　　设n<m，原式等于$sum_{i=1}^{n}sum_{j=1}^{m}i*j/gcd(i,j)$

　　然后枚举d值作为i和j的gcd，得到$$sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{m}frac{i*j}{d}[gcd(i,j)==d]$$

　　因为gcd(i,j)==d，所以gcd(i/d,j/d)==1，得到$$sum_{d=1}^{n}d*sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{m}{d}}i*j[gcd(i,j)==1]$$

　　因为莫反的性质：$sum_{d|x}mu(d)=[x==1]$，所以转化为$$sum_{d=1}^{n}d*sum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{m}{d}}i*j*sum_{t|gcd(i,j)}mu(t)$$

　　交换和式得到$$sum_{d=1}^{n} d* sum_{t=1}^{frac{n}{d}} mu(t) * sum_{i=1}^{frac{n}{d}} sum_{j=1}^{frac{m}{d}} i*j[gcd(i,j)==t]$$

　　$$sum_{d=1}^{n}d*sum_{t=1}^{frac{n}{d}}t^{2}*mu(t)*sum_{i=1}^{frac{n}{dt}}sum_{j=1}^{frac{m}{dt}}i*j$$

　　设$T=dt$，$S(x)=sum_{i=1}^{x}i$，得到$$sum_{d=1}^{n}sum_{t=1}^{frac{n}{d}}T*t*mu(t)*S(frac{n}{T})*S(frac{m}{T})$$

　　将$S(frac{n}{T})*S(frac{m}{T})$提前，得到$$sum_{T=1}^{n}T*S(frac{n}{T})*S(frac{m}{T})sum_{d|T}d*mu(d)$$

　　因为$S(frac{n}{T})*S(frac{m}{T})$可以前缀和预处理，显然我们只要将$sum_{d|T}d*mu(d)$快速求出就可以了

　　设$F(T)=sum_{d|T}d*mu(d)$，显然是一个积性函数，在线性筛的时候求就行了

　　然后将$T*F(T)$求前缀和，然后整除分块加速就能过了

---

参考代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long LL;
int prime[11000000],v[11000000];
LL f[11000000],sum[11000000];
LL Mod=1e8+9;
void pre(int n)
{
    f[1]=1;
    int tot=0;
    for(int i=2;i<=n;i++)
    {
        if(v[i]==0)
        {
            v[i]=i;
            prime[++tot]=i;
            f[i]=(1-i+Mod)%Mod;
        }
        for(int j=1;j<=tot;j++)
        {
            if(prime[j]>v[i]||prime[j]>n/i) break;
            v[i*prime[j]]=prime[j];
            if(i%prime[j]==0){f[i*prime[j]]=f[i]%Mod;break;}
            else f[i*prime[j]]=f[i]*f[prime[j]]%Mod;
        }
    }
    for(int i=1;i<=n;i++) f[i]=(f[i]*LL(i)%Mod+f[i-1])%Mod;
    for(int i=1;i<=n;i++) sum[i]=(sum[i-1]+LL(i))%Mod;
}
int main()
{
    pre(10000000);
    int T;
    scanf("%d",&T);
    while(T--)
    {
        int n,m;
        scanf("%d%d",&n,&m);
        if(n>m) swap(n,m);
        LL ans=0;
        for(int i=1,j;i<=n;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans=(ans+(f[j]-f[i-1]+Mod)%Mod*sum[n/i]%Mod*sum[m/i]%Mod)%Mod;
        }
        printf("%lld
",ans);
    }
    return 0;
}