【BZOJ2154】Crash的数字表格（莫比乌斯反演）

题意：

思路：如上

From http://blog.csdn.net/regina8023/article/details/44243911

最后的F(x,y)的推法和求gcd(x,y)=1的(x,y)对数差不多，只不过在推导过程中把原来1的地方换成x*y。

那么我们预处理出i^2*u[i]的前缀和，F(x,y)可以在O(sqrt(n))时间内求出来，而ans的式子也可以在O(sqrt(n))中求出来，

因此最后的时间复杂度是O(n)的。

 1 const mo=20101009;
 2 var flag,prime,mu:array[0..10000000]of longint;
 3     sum:array[0..10000000]of int64;
 4     n1,m1,i,j,max,m,t,x,y,t1,t2:longint;
 5     ans,pos:int64;
 6 
 7 function ask(n,m:int64):int64;
 8 begin
 9  n:=n*(n+1) div 2 mod mo;
10  m:=m*(m+1) div 2 mod mo;
11  exit(n*m mod mo);
12 end;
13 
14 function clac(n,m:longint):int64;
15 var t1,t2,i,pos:longint;
16     x,y:int64;
17 begin
18  i:=1; clac:=0;
19  while i<=n do
20  begin
21   x:=n div i; y:=m div i;
22   t1:=n div x; t2:=m div y;
23   if t1<t2 then pos:=t1
24    else pos:=t2;
25   clac:=clac+ask(x,y)*(sum[pos]-sum[i-1]) mod mo;
26   clac:=(clac mod mo+mo) mod mo;
27   i:=pos+1;
28  end;
29 end;
30 
31 begin
32  assign(input,'bzoj2154.in'); reset(input);
33  assign(output,'bzoj2154.out'); rewrite(output);
34  readln(n1,m1);
35  if n1>m1 then max:=n1
36   else max:=m1;
37  mu[1]:=1;
38  for i:=2 to max do
39  begin
40   if flag[i]=0 then
41   begin
42    inc(m); prime[m]:=i;
43    mu[i]:=-1;
44   end;
45   j:=1;
46   while (j<=m)and(prime[j]*i<=max) do
47   begin
48    t:=prime[j]*i; flag[t]:=1;
49    if i mod prime[j]=0 then
50    begin
51     mu[t]:=0;
52     break;
53    end;
54    mu[t]:=-mu[i];
55    inc(j);
56   end;
57  end;
58  for i:=1 to max do sum[i]:=(sum[i-1]+int64(mu[i])*i*i) mod mo;
59  if n1>m1 then
60  begin
61   t:=n1; n1:=m1; m1:=t;
62  end;
63  i:=1;
64  while i<=n1 do
65  begin
66   x:=n1 div i; y:=m1 div i;
67   t1:=n1 div x; t2:=m1 div y;
68   if t1<t2 then pos:=t1
69    else pos:=t2;
70   ans:=ans+(pos-i+1)*(pos+i) div 2 mod mo*clac(x,y) mod mo;
71   ans:=ans mod mo;
72   i:=pos+1;
73  end;
74  writeln(ans);
75  close(input);
76  close(output);
77 end.

UPD（2018.10.22）：C++

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<string>
 4 #include<cmath>
 5 #include<iostream>
 6 #include<algorithm>
 7 #include<map>
 8 #include<set>
 9 #include<queue>
10 #include<vector>
11 using namespace std;
12 typedef long long ll;
13 typedef unsigned int uint;
14 typedef unsigned long long ull;
15 typedef pair<int,int> PII;
16 typedef vector<int> VI;
17 #define fi first
18 #define se second
19 #define MP make_pair
20 #define MAXN 10000000 
21 #define eps 1e-8
22 #define pi  acos(-1)
23 #define oo  1e9
24 #define MOD 20101009
25  
26 int flag[MAXN+10],prime[MAXN+10],mu[MAXN+10];
27 ll s[MAXN+10];
28  
29 ll ask(ll n,ll m)
30 {
31     n=n*(n+1)/2%MOD;
32     m=m*(m+1)/2%MOD;
33     return n*m%MOD;
34 }
35  
36 ll calc(int n,int m)
37 {
38     int i=1; 
39     ll ans=0;
40     while(i<=n)
41     {
42         ll x=n/i; ll y=m/i;
43         ll t1=n/x; ll t2=m/y;
44         ll pos=min(t1,t2);
45         ans+=ask(x,y)*(s[pos]-s[i-1])%MOD;
46         ans=(ans%MOD+MOD)%MOD;
47         i=pos+1;
48     }
49     return ans;
50 }
51  
52 int main()
53 {
54     int N,M;
55     scanf("%d%d",&N,&M);
56     int MAX=max(N,M);
57     mu[1]=1;
58     int m=0;
59     for(int i=2;i<=MAX;i++)
60     {
61         if(!flag[i])
62         {
63             prime[++m]=i;
64             mu[i]=-1;
65         }
66         for(int j=1;j<=m;j++)
67         {
68             int t=prime[j]*i;
69             if(t>MAX) break; 
70             flag[t]=1;
71             if(i%prime[j]==0)
72             {
73                 mu[t]=0;
74                 break;
75             }
76             mu[t]=-mu[i];
77         }
78     }
79      
80     for(int i=1;i<=MAX;i++) s[i]=(s[i-1]+1ll*mu[i]*i*i)%MOD;
81     if(N>M) swap(N,M);
82     int i=1;
83     ll ans=0;
84     while(i<=N)
85     {
86         int x=N/i; int y=M/i;
87         int t1=N/x; int t2=M/y;
88         ll pos=min(t1,t2);
89         ans+=(pos-i+1)*(pos+i)/2*calc(x,y)%MOD;
90         ans%=MOD;
91         i=pos+1;
92     }
93     printf("%lld
",ans); 
94     return 0;
95 }
96