【51NOD 1223】分数等式的数量

Description

有这样一个分数等式：1/X + 1/Y = 1/N，(X,Y,N > 0)。给出L，求有多少满足X < Y <= L的等式。
例如：L = 12，满足条件的等式有3个，分别是：1/3 + 1/6 = 1/2, 1/4 + 1/12 = 1/3, 1/6 + 1/12 = 1/4。

Solution

$$\frac{1}{x}+\frac{1}{y}=\frac{1}{a}$$
$0<x<y$
$$\frac{x+y}{xy}=\frac{1}{a}$$
$$(x+y)|xy$$
设 $k=\gcd(x,y),a=x/k,b=y/k$
$$k(a+b)|k^2ab$$
$$(a+b)|kab$$
又因为ab，与(a+b)互质，
所以： $$(a+b)|k$$
有因为 $ak,bk<=n$，所以对于每一个 $a,b$，k的取值有 $\lfloor\frac{n}{(a+b)*b}\rfloor$这么多个，

于是把试子最终变成了：
设$m=\sqrt n$

$$\sum_{i=2}^{m}\sum_{j=1}^i \lfloor\frac{n}{(i+j)i}\rfloor(\gcd(i,j)=1)$$

看到有gcd，我们就要信心满满的上莫比乌斯反演，

把$i=1$的情况考虑进去，程序末尾减掉即可，

设：$f_d$表示gcd为d的倍数的（那一坨东西）的和，

$$f_d=\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{j=1}^i \lfloor\frac{n}{d^2(i+j)i}\rfloor$$
设： $g_d$表示gcd为d的（那一坨东西）的和，
$$f_d=\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}g_{id}$$
$$g_d=\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}f_{id}*\mu(i)$$
$$g_d=\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\mu(i)*\left(\sum_{j=1}^{\lfloor\frac{m}{id}\rfloor}\sum_{k=1}^j \lfloor\frac{n}{(id)^2(j+k)j}\rfloor \right)$$
于是：
$$Ans=\sum_{i=1}^m\mu(i)*\left(\sum_{j=1}^{\lfloor\frac{m}{i}\rfloor}\sum_{k=1}^j \lfloor\frac{n}{i^2(j+k)j}\rfloor \right)$$

发现，后面的可以分块，
而又因为$\mu$不是每个都有，所以远远没有$m\log(m)$那么大，

复杂度：$O(\sqrt n*\log(\sqrt n)*\sqrt{\sqrt n})$ (远远没有这么大)

Code

#include <cstdio>
#include <cstdlib>
#include <cmath>
#define fo(i,a,b) for(LL i=a;i<=b;i++)
using namespace std;
const int N=4e5;
typedef long long LL;
LL n,ans,m;
LL mu[N];
int b[100];
bool prz[N];
int pr[N/5];
LL min(LL q,LL w){return q<w?q:w;}
LL fk(LL T,LL j)
{
    LL k=j+1,M=2*j;
    LL ans=0;
    while(k<=M)
    {
        if(!(T/k))break;
        LL nx=min(T/(T/k),M);
        ans+=(T/k)*(nx-k+1);
        k=nx+1;
    }
    return ans;
}
int main()
{
    scanf("%lld",&n);m=sqrt(n);
    mu[1]=1;
    fo(i,2,m)
    {
        if(!prz[i])pr[++pr[0]]=i,mu[i]=-1;
        fo(j,1,pr[0])
        {
            int t=pr[j]*i;
            if(t>m)break;
            prz[t]=1;
            if(i%pr[j]==0)break;
            mu[t]=-mu[i];
        }
    }
    fo(i,2,m)mu[i]+=mu[i-1];
    LL i=1;
    while(i<=m)
    {
        LL nx=i,q=n/i/i;
        while(q==n/((nx+1)*(nx+1)))nx++;
        nx=min(nx,m/(m/i));
        LL t=0,T=n/i/i;
        fo(j,1,m/i)t+=fk(T/j,j);
        ans+=(mu[nx]-mu[i-1])*t;
        i=nx+1;
    }
    ans-=n/2;
    printf("%lld\n",ans);
    return 0;
}