[bzoj2154]Crash的数字表格

2154: Crash的数字表格

Time Limit: 20 Sec Memory Limit: 259 MB
Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122

【数据规模和约定】

100%的数据满足N, M ≤ 10^7。

这里我设n < m,不然就swap
题目要求
∑ni=1∑mj=1lcm(i,j)

=∑ni=1∑mj=1ijgcd(i,j)

发现这题的贡献和i，j都有关，并且和gcd有关系。那么我们发现难点主要在gcd，所以我们把gcd拆开做。
那么我们设
f(n,m,d)=∑ni=1∑mj=1i∗j∗[gcd(i,j)=d]

这么设的优点在哪里呢？
我们可以方便的表示出ans的式子
ans=∑nd=1f(n,m,d)d

=∑nd=1f(nd,md,1)∗d∗dd

=∑nd=1f(nd,md,1)∗d

所以我们现在的问题是怎么求f
根据经验，我们知道如果f中存在[xxx=d]的话，那么F中应当有一个[d|xxx]的状况，然后就能倍数和乱搞。
构造F(n,m,d)=∑ni=1∑mj=1i∗j∗[d|gcd(i,j)]使得我们能够反演

现在我们对F进行操作
F(n,m,d)=∑n/di=1∑m/dj=1i∗j∗d2

=d2∗∑n/di=1∑m/dj=1i∗j
通过这个操作，我们去掉了d|gcd(i,j)的条件

对于∑n/di=1∑m/dj=1i∗j

=∑n/di=1(1+md)∗md2

=(1+nd)∗nd2(1+md)∗md2

那么显然这个数据范围大F是不能出现O(n)求的情况，因为ans外面已经套了个∑了，所以必须分段求和优化
设sum(n/d,m/d)=这个式子

所以F(n,m,d)的式子就出来了

F(n,m,d)=d2∗sum(n/d,m/d)这就能分段求和了

那么可以真正上反演了
F(n,m,d)=∑d|nf(n)
f(n,m,d)=∑d|nμ(n/d)F(d)

f(n,m,1)=∑nd=1μ(d)F(d)

=f(n,m,1)=∑nd=1μ(d)∗d2∗sum(n/d,m/d)

这里ans还是不能直接求，发现ans式子中f有除号。那这里分段求和一下就好了。

#include<iostream>#include<cstring>#include<string>#include<cstdlib>#include<cmath>#include<algorithm>#include<cstdio>#define ll long longusing namespace std;const int N=11000000;const int pyz=20101009;int n,m;int prime[N],mu[N],sqri[N];bool is[N];inline int read(){    char c;    bool flag=0;    while((c=getchar())>'9'||c<'0') if(c=='-')flag=1;    int res=c-'0';    while((c=getchar())>='0'&&c<='9') res=(res<<3)+(res<<1)+c-'0';    return flag?-res:res;}inline int sum(ll x,ll y){    return (((x+1)*x/2)%pyz)*(y*(y+1)/2%pyz)%pyz;}inline ll F(int x,int y){    int j=0;    ll ans=0;    for(int i=1;i<=x;i=j+1)    {        j=min(x/(x/i),y/(y/i));        ans=(ans+(ll)(sqri[j]-sqri[i-1]+pyz)%pyz*sum((ll)x/i,(ll)y/i)%pyz)%pyz;    }    return ans;}int main(){//  freopen("2154.txt","r",stdin);    n=read();m=read();    if(n>m) swap(n,m);    mu[1]=1;    for(int i=2;i<=n;++i)    {        if(!is[i]) prime[++prime[0]]=i,mu[i]=-1;        for(int j=1;j<=prime[0]&&i*prime[j]<=n;++j)        {            is[i*prime[j]]=1;            if(!(i%prime[j]))            {                mu[i*prime[j]]=0;;                break;            }            mu[i*prime[j]]=-mu[i];        }    }    for(int i=1;i<=n;++i) sqri[i]=(sqri[i-1]+(ll)i*i%pyz*mu[i])%pyz;     int j=0;    ll ans=0;    for(int d=1;d<=n;d=j+1)    {        j=min(n/(n/d),m/(m/d));        ans=(ans+(ll)(d+j)*(j-d+1)/2%pyz*F(n/d,m/d)%pyz)%pyz;    }    printf("%lld",ans);}