bzoj2671: Calc

链接

　　http://www.lydsy.com/JudgeOnline/problem.php?id=2671

题解

　　天啊噜！这道题好生麻烦
　　我尽量长话短说。
　　观察一下条件：

(a+b)|ab

　　我们数论的算法都是和gcd有关的，因此肯定要转到gcd上来，令d=gcd(a,b)，a=xd,b=yd，写成分数
　　

xyd2d(x+y)

　　即

xydx+y

　　因为gcd(x,y)=1，所以gcd(x,x+y)=1(很显然啊因为提不出公因式)，同理gcd(y,x+y)=1，因此gcd(xy,x+y)=1，所以想要让上式为整数，必有(x+y)|d。而这样的d有多少个？先不要匆忙下结论说有⌊Nx+y⌋个，因为d跟x、y是有关联的。
　　根据条件可以列出以下式子
　　

⎧⎩⎨xd≤Nyd≤Nx+y≤d

　　将第三个式子左右同乘以y，得y(y+x)≤N
　　而x的最小值是1，因此y(y+1)≤N。可以发现，y是N−−√级别的，因此可以枚举。而d≤⌊Ny⌋，所以d是有上界的。
　　那么写出答案
　　

ans=∑x=1N√∑y=x+1N√⌊⌊Ny⌋x+y⌋[gcd(x,y)=1]

　　这样的复杂度是稳定的O(NlogN)。那个log是求gcd时候的复杂度。实际测试：TLE。
　　下面进入丧心病狂的部分，首先进行莫比乌斯反演：
　　

∑x=1N√∑y=x+1N√⌊⌊Ny⌋x+y⌋∑d|gcd(x,y)μ(d)

　　

∑d=1N√μ(d)∑x=1⌊N√d⌋∑y=x+1⌊N√d⌋⌊⌊Nyd⌋y(x+y)⌋

　　注意这里的x和y已经不是原来的x和y了，现在的xd和yd分别对应以前的x和y。
　　

∑d=1N√μ(d)∑x=1⌊N√d⌋∑y=x+1⌊N√d⌋⌊⌊Nyd2⌋x+y⌋

　　这样比原来稍稍快一点，复杂度还是O(N1.5)实际测试T掉了但是速度的确比原来快了接近10倍（有时候光看复杂度也不可靠啊），跑最大的数据在我的笔记本上只用5s，看来再加一个小优化就差不多能过了呢。
　　丧病的优化，让我想了一整节晚自习，令s=x+y：
　　

∑d=1N√μ(d)∑y=2⌊N√d⌋∑s=1+y2y−1⌊⌊Nyd2⌋s⌋

　　分块，复杂度变为O(N1.25)，别看这个复杂度，常数可是极小的，真实的复杂度应该比这个小的多，不过我懒得算了….bzoj上1628s过，最大数据在本机只跑0.4511s

代码

//莫比乌斯反演 #include <cstdio>#include <algorithm>#include <cmath>#define maxn 50000#define ll long longusing namespace std;ll prime[maxn+10], mu[maxn+10], N, mark[maxn+10];void init(){    ll i, j;    mu[1]=1;    for(i=2;i<=maxn;i++)    {        if(!mark[i])prime[++prime[0]]=i,mu[i]=-1;        for(j=1;j<=prime[0] and i*prime[j]<=maxn;j++)        {            mark[i*prime[j]]=1;            if(i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            mu[i*prime[j]]=-mu[i];        }    }}void work(){    ll ans=0, d, t, sqrtn=sqrt(N), tmp, y, last, s;    for(d=1;d*(d+1)<=N;d++)    {        tmp=0;        for(y=2;y<=sqrtn/d;y++)        {            t=N/(d*d*y);            for(s=1+y;s<=y+y-1 and s<=t;s=last+1)            {                last=t/(t/s);                tmp+=(min(last,y+y-1)-s+1)*(t/s);            }        }        ans+=tmp*mu[d];    }    printf("%lld",ans);}int main(){    init();    scanf("%lld",&N);    work();    return 0;}