bzoj2671: Calc

链接

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

题解

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

$$(a+b)|ab$$
　　我们数论的算法都是和 $gcd$有关的，因此肯定要转到 $gcd$上来，令 $d=gcd(a,b)$， $a=xd,b=yd$，写成分数
　　 $$\frac{xyd^2}{d(x+y)}$$
　　即 $$\frac{xyd}{x+y}$$
　　因为 $gcd(x,y)=1$，所以 $gcd(x,x+y)=1$(很显然啊因为提不出公因式)，同理 $gcd(y,x+y)=1$，因此 $gcd(xy,x+y)=1$，所以想要让上式为整数，必有 $(x+y)|d$。而这样的 $d$有多少个？先不要匆忙下结论说有$\lfloor\frac{N}{x+y}\rfloor$个，因为 $d$跟$x、y$是有关联的。
　　根据条件可以列出以下式子
　　 $$\begin{cases}xd\le N \\ yd\le N\\ x+y\le d\end{cases}$$
　　将第三个式子左右同乘以 $y$，得$y(y+x)\le N$
　　而 $x$的最小值是1，因此$y(y+1)\le N$。可以发现， $y$是$\sqrt N$级别的，因此可以枚举。而 $d\le \lfloor\frac{N}{y}\rfloor$，所以d是有上界的。
　　那么写出答案
　　 $$ans=\sum_{x=1}^{\sqrt N}\sum_{y=x+1}^{\sqrt N}\lfloor\frac{\lfloor\frac{N}{y}\rfloor}{x+y}\rfloor[gcd(x,y)=1]$$
　　这样的复杂度是稳定的 $O(NlogN)$。那个 $log$是求 $gcd$时候的复杂度。实际测试：TLE。
　　下面进入丧心病狂的部分，首先进行莫比乌斯反演：
　　 $$\sum_{x=1}^{\sqrt N}\sum_{y=x+1}^{\sqrt N}\lfloor\frac{\lfloor\frac{N}{y}\rfloor}{x+y}\rfloor\sum_{d|gcd(x,y)}\mu(d)$$
　　 $$\sum_{d=1}^{\sqrt N}\mu(d)\sum_{x=1}^{\lfloor\frac{\sqrt N}{d}\rfloor}\sum_{y=x+1}^{\lfloor\frac{\sqrt N}{d}\rfloor}\lfloor\frac{\lfloor\frac{ N}{yd}\rfloor}{(xd+yd)}\rfloor$$
　　注意这里的 $x$和$y$已经不是原来的 $x$和$y$了，现在的 $xd$和 $yd$分别对应以前的 $x$和$y$。
　　 $$\sum_{d=1}^{\sqrt N}\mu(d)\sum_{x=1}^{\lfloor\frac{\sqrt N}{d}\rfloor}\sum_{y=x+1}^{\lfloor\frac{\sqrt N}{d}\rfloor}\lfloor\frac{\lfloor\frac{ N}{yd^2}\rfloor}{x+y}\rfloor$$
　　这样比原来稍稍快一点，复杂度还是 $O(N^{1.5})$实际测试T掉了但是速度的确比原来快了接近10倍（有时候光看复杂度也不可靠啊），跑最大的数据在我的笔记本上只用5s，看来再加一个小优化就差不多能过了呢。
　　丧病的优化，让我想了一整节晚自习，令 $s=x+y$：
　　 $$\sum_{d=1}^{\sqrt N}\mu(d)\sum_{y=2}^{\lfloor\frac{\sqrt N}{d}\rfloor}\sum_{s=1+y}^{2y-1}\lfloor\frac{\lfloor\frac{N}{yd^2}\rfloor}{s}\rfloor$$
　　分块，复杂度变为 $O(N^{1.25})$，别看这个复杂度，常数可是极小的，真实的复杂度应该比这个小的多，不过我懒得算了….bzoj上1628s过，最大数据在本机只跑0.4511s

代码

//莫比乌斯反演 
#include <cstdio>
#include <algorithm>
#include <cmath>
#define maxn 50000
#define ll long long
using 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;
}