Description
Input
Output
Sample Input
10 10
100 100
Sample Output
2791
HINT
T = 10000
N, M <= 10000000
文章中/都是整除
首先讲一下求gcd(x,y)=1 (1<=x<=a,1<=y<=b)的二元组(x,y)的个数的方法
这个可以容斥求
但是模拟容斥太慢
为了简化过程,有了个mobius函数
mob(i)=1 (i=1),0 (i 不是无平方因子数,就是说i的质因数指数都为1),-1^k (k是i的质因子个数)
这个实际上就是把用到的数的系数提前算出来……就不用在模拟容斥的过程中算了……而且可以利用这个来简化过程
设f[a,b,k]为gcd(x,y)=wk (1<=x<=a,1<=y<=b,w为任意正整数)的对数
发现f[a,b,k]十分好算,f[a,b,k]=(a/k)*(b/k)
这样容斥算好啦
设gcd(x,y)=1 (1<=x<=a,1<=y<=b)的二元组(x,y)的个数为ans
ans=sigma(mob(k)*f[a,b,k])
可是这样需要枚举k,复杂是平方级的,太慢
来观察a/k和b/k
发现由于都是整除
a/k对于许多k取值是相同的
可证明a/k的取值是sqrt(a)级别的
b/k也同理
那么a/k的取值加上b/k的取值难道是sqrt(a)*sqrt(b)级别的?
显然不是……中学毕业的都知道反比例函数是递增的
所a/k*b/k的取值至多有sqrt(a)+sqrt(b)种,且每种都是连续一段
可以对于a/k*b/k的我们合并同类项,把其系数加起来再乘……
这样我们把mob函数做前缀和好了
对于a/k相同的长度有点难算……
不过我们如果知道左界l可以确定右界r=a/(a/l)
至此这个问题就解决了
接下来说这道题
有个变化,求的不是gcd=1的了,求的是质数
要是是给定的gcd=k,还可以转化成gcd(x,y)=1 (1<=x<=a/k,1<=y<=b/k)来算……
枚举素数由于素数非常多,所以慢的要死
怎么解决呢
设本问题答案为ans
ans=sigma(mob(k)*(a/k/d)*(b/k/d)) (d is prime)
首先说下……整除也是满足结合律的……a/k/d=a/(k*d)……我一直都不知道
观察一下
这个a和b由于是读入的,所以不好动手脚
只能处理系数了
转换一下枚举对象
来枚举k*d,设x=k*d
发现ans=simga(sigma(mob(i)){i|x and i is prime}*(a/x)*(b/x))
…………
sigma(mob(i)){i|x and i is prime}
处理出来这个就跟前头那个一样了嘛……
设这个函数为g(x)好了
g(x)不是很好处理
利用下mob函数的性质
若x有某质因数指数大于等于3或者有两个及以上等于2的话,无论d是什么mob(x/d)都=0,所以g(x)=0
若x有且仅有某质因数v指数等于2,g(x)=mob(x/v)
若x为无平方因子数,若有k个质因数,显然任意mob(x/d)都是相等的,g(x)=(-1)^(k+1)*k
具体实现比较蛋疼,详见代码
program bzoj2820;
const
maxn=10000000;
var
o,t,n,m,i,j,k:longint;
f,p:array [0..maxn] of longint;
w:array [0..maxn] of boolean;
function calc (a,b:longint):int64;inline;
var
i,k:longint;
ans:int64;
begin
if a>b then
begin
a:=a xor b;
b:=a xor b;
a:=a xor b;
end;
if a=0 then exit(0);
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+int64(p[k]-p[i-1])*(a div i)*(b div i);
i:=k+1;
end;
exit(ans);
end;
begin
for i:=2 to maxn do
if f[i]=0 then
for j:=i to maxn div i do
f[i*j]:=i;
p[1]:=0;
for i:=2 to maxn do
if f[i]=0 then p[i]:=1
else
begin
o:=i div f[i];
w[i]:=w[o];
if (p[o]=0)or(w[o] and(o mod f[i] = 0)) then
p[i]:=0
else
if o mod f[i] = 0 then
begin
p[i]:=-p[o] div abs(p[o]);
w[i]:=true;
end
else
if w[o] then
p[i]:=-p[o]
else
p[i]:=-p[o] div abs(p[o]) * (abs(p[o])+1);
end;
for i:=2 to maxn do
p[i]:=p[i-1]+p[i];
read(t);
while t>0 do
begin
dec(t);
read(n,m);
writeln(calc(n,m));
end;
end.