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

 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.