[bzoj2154]Crash的数字表格

来自FallDream的博客,未经允许,请勿转载,谢谢。


题意:求$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)   \\\    n,m\leqslant 10^{7}$$

之前有做过一道类似的题目,nm相同,但是最大有$10^{10}$,可以用欧拉函数+杜教筛解决。

但是这道题n与m不同,欧拉函数就没办法排上用场了,还是换换思路吧。由于求的是lcm,所以考虑枚举gcd是p。假装n<m,得到

$$Ans=\sum_{p=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{p}\rfloor}p*i*j\sum_{d|i,d|j}\mu(d)$$

把d提出来,发现计算的次数是一个等差数列两两相乘,令$G(n)=\sum_{i=1}^{n}i$,那么

$$Ans=\sum_{p=1}^{n}p*\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}\mu(d)*d^{2}*S(\lfloor\frac{n}{pd}\rfloor)*S(\lfloor\frac{m}{pd}\rfloor)$$

然后老套路,令T=pd,得到

$$Ans=\sum_{T=1}^{n}T*S(\lfloor\frac{n}{T}\rfloor)*S(\lfloor\frac{m}{T}\rfloor)*\sum_{d|T}mu(d)*d$$

所以我们只需要预处理出$F[T]=\sum_{d|T}mu(d)*d$就行了。

嗯我推到这里感觉很正确就交了一个类似欧拉筛的东西上去,妥妥T掉了....

但其实并不需要,$F[T]$其实是积性函数,这个可以自己脑补,然后可以线筛出来。

假设p是质数,那么p|i的时候,$F[i*p]=F[i]$,因为新的项的μ都是0.

如果不能整除,那么$F[i*p]=F[i]*F[p]$

这样就可以过啦。最后的计算可以每次$\sqrt{n}$做,当然直接暴力也是可以的。

#include<iostream>
#include<cstdio>
#include<ctime>
#define MN 10000000
#define mod 20101009
#define ll long long
using namespace std;
inline int read()
{
    int x = 0 , f = 1; char ch = getchar();
    while(ch < '0' || ch > '9'){ if(ch == '-') f = -1;  ch = getchar();}
    while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0';ch = getchar();}
    return x * f;
}

int n,m,s[MN],num=0,ans=0,f[MN+5];
bool b[MN+5];

inline int calc(int x,int y){return 1LL*(1+(x/y))*(x/y)%mod*10050505%mod;}

int main()
{
    f[1]=1;
    for(register int i=2;i<=MN;++i)
    {
        if(!b[i]) s[++num]=i,f[i]=mod+1-i;
        for(int j=1;s[j]*i<=MN;++j)
        {
            b[s[j]*i]=1;
            if(i%s[j]==0) {f[s[j]*i]=f[i];break;}
            f[s[j]*i]=1LL*f[i]*f[s[j]]%mod;
        }
    }
    for(register int i=2;i<=MN;++i) f[i]=(1LL*f[i]*i+f[i-1])%mod;
    n=read();m=read();if(n>m) swap(n,m);
    for(register int i=1,last;i<=n;i=last+1)
    {
        last=min(m/(m/i),n/(n/i));
        ans=(ans+1LL*(f[last]-f[i-1]+mod)*calc(n,i)%mod*calc(m,i)%mod)%mod;
    }
    cout<<ans;
    return 0;
}

 

posted @ 2017-04-19 20:18  FallDream  阅读(244)  评论(0编辑  收藏  举报