学步园

Description

Input

Output

Sample Input

Sample Output

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.