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
Sample Output
122
【数据规模和约定】
100%的数据满足N, M ≤ 107。
jzp神题
似乎跟gcd问题有点像，实际上区别较大
求sigma(lcm(x,y))(1<=x<=n,1<=y<=m)
明显的，枚举最大公约数d之后就变为sigma(x*y)(1<=x<=n/d,1<=y<=m/d,gcd(x,y)=1)=f[d]
那么gcd=d的lcm和为f[d]*d
虽然d的取值在[1,n]但是n/d和 m/d的取值是sqrt(n)级别的
就可以用解决gcd问题的类似的方法来解决这个假设l到r的n/d和m/d是相等的，我们可以将其一起统计
对答案的贡献为sigma(x*y)(1<=x<=n/d,1<=y<=m/d,gcd(x,y)=1)*(r-l+1)*(l+r)/2
接下来问题是如何快速求sigma(x*y)(1<=x<=n/d,1<=y<=m/d,gcd(x,y)=1)
设这个为g[a,b]好了
那么类似的我们可以求f[a,b,k]=sigma(xy)(1<=x<=a,1<=y<=b,gcd[x,y]=wk)
明显f[a,b,k]可以直接算
f[a,b,k]=(1+..+a/k)*k*(1+..+b/k)*k=(a/k+1)*(a/k)/2*k*(b/k+1)*(b/k)/2*k
显然g[a,b]=sigma(mob(k)*f[a,b,k])
问题来了
每个f[a,b,k]都含有k^2这一项,这个变换不连续……
但是显然的f[a,b,k]只会乘mob(k)
那我们定义o(i)=mob(i)*i*i,w(a,b,k)=(a/k+1)*(a/k)/2*(b/k+1)*(b/k)/2
g[a,b]=sigma(o(k)*w(a,b,k))
发现w(a,b,k)是分sqrt段的
o可以预处理出来
就可以sqrt级别时间求g函数了
那么复杂度为O(sqrt(n)*sqrt(n))=O(n)
program bzoj2154;
const
maxn=20101009;
var
n,m,i,j,k:longint;
ans:int64;
f:array [0..10000001] of longint;
p:array [0..10000001] of int64;
function calc (a,b:longint):int64;inline;
var
i,k:longint;
ans:int64;
begin
ans:=0;
i:=1;
while i<=a do
begin
if a div (a div i)<b div (b div i) then
k:=a div (a div i)
else
k:=b div (b div i);
ans:=(ans+(p[k]-p[i-1]) mod maxn*
(int64(a div i+1)*(a div i) div 2 mod maxn) mod maxn*
(int64(b div i+1)*(b div i) div 2 mod maxn) mod maxn) mod maxn;
i:=k+1;
end;
exit(ans);
end;
begin
read(n,m);
if n>m then
begin
n:=n xor m;
m:=n xor m;
n:=n xor m;
end;
for i:=2 to n do
if f[i]=0 then
for j:=i to n div i do
f[i*j]:=i;
p[1]:=1;
for i:=2 to n do
if f[i]=0 then p[i]:=-1
else
begin
k:=i div f[i];
if k mod f[i] = 0 then p[i]:=0
else p[i]:=-p[k];
end;
for i:=1 to n do
p[i]:=p[i]*i*i mod maxn;
for i:=2 to n do
p[i]:=(p[i]+p[i-1]) mod maxn;
i:=1;
while i<=n do
begin
if n div (n div i)<m div (m div i) then
k:=n div (n div i)
else
k:=m div (m div i);
ans:=(ans+int64(k-i+1)*(k+i) div 2 mod maxn * calc(n div i,m div i)) mod maxn;
i:=k+1;
end;
ans:=(ans + maxn) mod maxn;
writeln(ans);
end.