[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。
L<=10^11

Solution

化简一下式子。

$$1/x + 1/y = 1/N$$
那么 $$Nx+Ny=xy$$
都是整数
$$x+y|xy$$
然后就是套路。
$$d=gcd(x,y)$$
$x=di,y=dj$
$$d(i+j)|d^2ij,gcd(i,j)=1$$
因为
$$gcd(ij,i+j)=1$$
所以 $$i+j|d$$
设 $d=k(i+j)$
则有 $$x=k(i+j)i,y=k(i+j)j$$
再设 $m=\sqrt L$
$$Ans=\sum\limits_{i=2}^{m}\sum\limits_{j=1}^{i-1}\lfloor{L\over (i+j)i}\rfloor[gcd(i,j)=1]$$
因为当 $i=j$, $gcd(i,j)!=1$，所以为了方便，将 $j$的范围设到$i$
$$Ans=\sum\limits_{i=2}^{m}\sum\limits_{j=1}^{i}\lfloor{L\over (i+j)i}\rfloor[gcd(i,j)=1]$$

设

$$f(d)=\sum\limits_{i=2}^{m}\sum\limits_{j=1}^{i}\lfloor{L\over (i+j)i}\rfloor[gcd(i,j)=d]$$
$Ans=f(1)$
LF同志的名言：看到GCD，就要信心满满上反演！
设
此时我们暂时加上 $i=1$的情况，最后再特判一下就好。
$$g(d)=\sum\limits_{k=1}^{\lfloor{m\over d}\rfloor}f(kd)=\sum\limits_{i=1}^{\lfloor{m\over d}\rfloor}\sum\limits_{j=1}^{i}\lfloor{L\over (i+j)id^2}\rfloor$$
$$f(d)=\sum\limits_{k=1}^{\lfloor{m\over d}\rfloor}\mu(k)g(kd)=\sum\limits_{k=1}^{\lfloor{m\over d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor{m\over kd}\rfloor}\sum\limits_{j=1}^{i}\lfloor{L\over (i+j)i(kd)^2}\rfloor$$
$$Ans=f(1)=\sum\limits_{k=1}^{m}\mu(k)\sum\limits_{i=1}^{\lfloor{m\over k}\rfloor}\sum\limits_{j=1}^{i}\lfloor{L\over (i+j)ik^2}\rfloor$$

前后都分块，复杂度是?
式子有点复杂，带了几个根号，不过它开了6s。

LF同志说过，人要有梦想！

Code

#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cmath>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define LL long long 
#define N 400005
using namespace std;
LL n,s,pr[N],mu[N],m,l,sum[N];
bool bz[N];
void prp()
{
    mu[1]=sum[1]=1;
    fo(i,2,m)
    {
        if(!bz[i]) mu[i]=-1,pr[++l]=i;
        for(int j=1;j<=l&&pr[j]*i<=m;j++)
        {
            bz[i*pr[j]]=1;
            if(i%pr[j]==0)
            {
                mu[i*pr[j]]=0;
                break;
            }
            mu[i*pr[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
}
int main()
{
    cin>>n;
    memset(bz,0,sizeof(bz));
    m=sqrt(n);
    l=0;
    prp();
    LL k=1,s=0;
    while(k<=m)
    {
        LL k1=min(m/(m/k),(LL)sqrt(n/(n/(k*k)))),i=1,s1=0;
        if(n/(k*k)==0) break;
        if(k==1) i=2;
        while(i<=m/k)
        {
            LL j=1;
            while(j<=i)
            {
                LL v=(i+j)*i*k*k;
                if(n/v==0) break;
                LL j1=min(n/(n/v)/k/k/i-i,i);
                s1+=(j1-j+1)*(n/v);
                j=j1+1;
            }
            i++;
        }
        s+=s1*(sum[k1]-sum[k-1]);
        k=k1+1;
    }
    cout<<s;
}