[bzoj 2154]Crash的数字表格
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。
【数据规模和约定】
100%的数据满足N, M ≤ 10^7。
题解:
∑i∑jlcm(i,j) (i<=n,j<=m)
=∑i∑j i*j/gcd(i,j)
=∑d∑i∑j i*j[gcd(i,j)=d]/d
=∑d∑i∑j i*j*d*d[gcd(i,j)=1]/d (i<=n/d,j<=m/d)
=∑dd∑i∑j i*j∑pμ(p) (p|i,p|j)
=∑dd∑pμ(p)∑ii*p∑jj*p (p<=n/d,i<=n/dp,j<=m/dp)
令x=[n/dp],y=[m/dp]
=∑dd∑pμ(p)*p*p*((x+1)*x/2)*((y+1)*y/2)
手写公式好累啊,求推荐
接下来就可以分块了
令F(n)=∑pμ(p)*p*p*((x+1)*x/2)*((y+1)*y/2) (p<=n)
因为n/d的值可以分块,所以在外层分块i~pos,调用F(n/i)
在F()内再分块,因为(n/d)/p也可以分块
还有一个改动就是∑dd∑pμ(p)*p*p*((x+1)*x/2)*((y+1)*y/2)
这也不难,因为每一块内的值都相同,对于前面的d,套用等差数列求和再*F(n/i)即可
对于p*p,同理,在维护莫比乌斯函数μ前缀和时稍作修改
mu[i]+=mu[i-1]-> mu[i]=mu[i-1]+mu[i]*i*i
对于分块方法不懂的话
1 #include<iostream> 2 #include<cstdio> 3 #include<algorithm> 4 #include<cstring> 5 using namespace std; 6 typedef long long lol; 7 int Mod=20101009; 8 lol mu[10000001]; 9 int tot; 10 lol prime[5000001]; 11 bool vis[10000001]; 12 lol n,m,ans; 13 void mobius() 14 {lol i,j; 15 mu[1]=1; 16 for (i=2;i<=m;i++) 17 { 18 if (vis[i]==0) 19 { 20 tot++; 21 prime[tot]=i; 22 mu[i]=-1; 23 } 24 for (j=1;j<=tot,i*prime[j]<=m;j++) 25 { 26 vis[i*prime[j]]=1; 27 if (i%prime[j]==0) 28 { 29 mu[i*prime[j]]=0; 30 break; 31 } 32 mu[i*prime[j]]=-mu[i]; 33 } 34 } 35 for (i=1;i<=m;i++) 36 mu[i]=(mu[i-1]+(mu[i]*(i*i)%Mod)%Mod)%Mod; 37 } 38 lol min(lol x,lol y) 39 { 40 if (x<y) return x; 41 else return y; 42 } 43 lol sum(lol x,lol y) 44 { 45 lol s1=((1+x)*x/2)%Mod,s2=(((1+y)*y/2)%Mod); 46 return s1*s2%Mod; 47 } 48 lol work(lol x,lol y) 49 { 50 lol s=0; 51 lol pos=1,i; 52 for (i=1;i<=x;i=pos+1) 53 { 54 pos=min(x/(x/i),y/(y/i)); 55 s=(s+(((mu[pos]-mu[i-1]+Mod)%Mod)*(sum(x/i,y/i)))%Mod)%Mod; 56 } 57 return s; 58 } 59 int main() 60 { 61 cin>>n>>m; 62 if (n>m) swap(n,m); 63 mobius(); 64 lol pos=1,i; 65 for (i=1;i<=n;i=pos+1) 66 { 67 pos=min(n/(n/i),m/(m/i)); 68 ans=(ans+((((i+pos)*(pos-i+1)/2)%Mod)*work(n/i,m/i))%Mod)%Mod; 69 } 70 cout<<ans; 71 }