【BZOJ 2154】Crash的数字表格 （莫比乌斯+分块）

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。
 
 
　　这题用了两次分块了~~ 好高级...不过还不是多组的呢~
　　好吧我也还是看题解的，还不是很会推~~
　　http://blog.csdn.net/regina8023/article/details/44243911
　　
　　
最后一步的推法类似 求gcd（a,b）=1 的对数，1改成a*b即可，如下：
   
 
所以进行两次分块，两次都是√n，总时间复杂度：O(n)。
 
代码如下：
  1 #include<cstdio>
  2 #include<cstdlib>
  3 #include<cstring>
  4 #include<iostream>
  5 #include<algorithm>
  6 #include<queue>
  7 #include<cmath>
  8 using namespace std;
  9 #define Mod 20101009
 10 #define Maxn 10000010
 11 #define LL long long
 12 
 13 LL mu[Maxn],pri[Maxn],h[Maxn],pl;
 14 bool q[Maxn];
 15 
 16 LL mymin(LL x,LL y) {return x<y?x:y;}
 17 
 18 void get_mu(LL mx)
 19 {
 20     pl=0;
 21     memset(q,1,sizeof(q));
 22     mu[1]=1;
 23     for(LL i=2;i<=mx;i++)
 24     {
 25         if(q[i])
 26         {
 27             pri[++pl]=i;
 28             mu[i]=-1;
 29         }
 30         for(LL j=1;j<=pl;j++)
 31         {
 32             if(i*pri[j]>mx) break;
 33             q[i*pri[j]]=0;
 34             if(i%pri[j]==0) mu[i*pri[j]]=0;
 35             else mu[i*pri[j]]=-mu[i];
 36             if(i%pri[j]==0) break;
 37         }
 38     }
 39     h[1]=(mu[1]*1*1)%Mod;
 40     for(LL i=2;i<=mx;i++) h[i]=(h[i-1]+mu[i]*i*i)%Mod;
 41 }
 42 
 43 LL get_g(LL x,LL y)
 44 {
 45     return ( ( ((x+1)*x/2)%Mod )*( ((y+1)*y/2)%Mod ) )%Mod;
 46 }
 47 
 48 LL get_f(LL n,LL m)
 49 {
 50     LL t;
 51     if(n>m) t=n,n=m,m=t;
 52     
 53     LL sq=(LL)ceil(sqrt((double)m));
 54     
 55     LL ans=0;
 56     for(LL i=1;i<=mymin(n,sq);i++)
 57     {
 58         ans=( ans+((mu[i]*i*i)%Mod)*get_g(n/i,m/i) )%Mod;
 59     }
 60     for(int i=sq+1;i<=n;)
 61     {
 62         int x=n/i,y=m/i;
 63         int r1=n/x+1,r2=m/y+1;
 64         if(r1>n+1) r1=n+1;
 65         if(r2>n+1) r2=n+1;
 66         int r=mymin(r1,r2);
 67         ans=(ans+ ((h[r-1]+Mod-h[i-1])%Mod)*get_g(x,y) )%Mod;
 68         i=r;
 69     }
 70     return ans;
 71 }
 72 
 73 LL get_ans(int n,int m)
 74 {
 75     LL ans=0;
 76         
 77     LL sq=(LL)ceil(sqrt((double)m));
 78     for(LL i=1;i<=mymin(sq,n);i++)
 79     {
 80         ans=(ans+i*get_f(n/i,m/i) )%Mod;
 81     }
 82     
 83     for(LL i=sq+1;i<=n;)
 84     {
 85         LL x=n/i,y=m/i;
 86         LL r1=n/x+1,r2=m/y+1;
 87         LL r=mymin(r1,r2);
 88         if(r>m+1) r=m+1;
 89         ans=( ans+(((r-i)*(i+r-1)/2)%Mod)*get_f(x,y) )%Mod;
 90         i=r;
 91     }
 92     return ans;
 93 }
 94 
 95 int main()
 96 {
 97     int T;
 98     T=1;
 99     
100     while(T--)
101     {
102         LL n,m,t;
103         scanf("%lld%lld",&n,&m);
104         if(n>m) t=n,n=m,m=t;
105         get_mu(m);
106         
107         LL ans=get_ans(n,m);
108         
109         printf("%lld\n",ans);
110     }
111     return 0;
112 }
[BZOJ2154]

 

2016-08-30 16:00:42

posted @ 2016-08-30 15:57  konjak魔芋  阅读(234)  评论(0)  编辑  收藏  举报