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。

Input

输入1个数L(1 <= L <= 10^11)

Output

输出符合条件的等式的数量。

Input示例

12

Output示例

3

Solution

求${1\over x}+{1\over y}={1\over n}$的数量
化简

$${x+y\over xy}={1\over n}$$
即
$$x+y | xy$$
设 $d=gcd(x,y),a=x/d,b=y/d$
于是
$$(a+b)d|abd^2$$
$$(a+b)|abd$$
因为 $gcd(a,b)=1$所以 $(a+b)|d$
在我知道a和b的情况 $（a<b）$下，d的取值范围有 $\lfloor {n\over (a+b)b} \rfloor$个
设 $m=\sqrt n$
$$ans=\sum _{i=2} ^m \sum _{j=1}^{i-1} \lfloor {n\over (a+b)b} \rfloor [gcd(i,j)=1]$$
这里可以用欧拉函数解决，而我用莫比乌斯反演
先忽略i=1的情况，后面减掉即可
设f[d]表示gcd=d的倍数， $\lfloor {n\over (a+b)b} \rfloor$的和
$$f[d]=\sum_{i=1}^{\lfloor {m\over d} \rfloor}\sum_{j=1}^i\lfloor {n\over (a+b)bd^2} \rfloor$$
为什么这样？自己想，与上面一样的思路
设g[d]表示gcd=d时， $\lfloor {n\over (a+b)b} \rfloor$的和
$$g[d]=\sum _{i=1}^{\lfloor {m\over d} \rfloor}f[id]*\mu(i)$$
g[1]就是答案
$$ans=g[1]=\sum _{i=1}^m\mu(i)(\sum _{j=1}^{\lfloor {m\over i} \rfloor}\sum_{k=1}^j\lfloor {n\over (j+k)ji^2} \rfloor)$$

右边括号里的东西可以分块，也就是在程序中枚举i可以分块，然后枚举j，然后枚举(j+k)注意不是k，这时(j+k)可以分块
时间复杂度$O(\sqrt n*log_2(\sqrt n)*\sqrt {\sqrt n})$

Code

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define ll long long
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define M 320000
using namespace std;
ll n,m,ans=0,mu[M];
int bz[M],p[M];
int main()
{
    scanf("%lld",&n);m=(int)sqrt(n);
    mu[1]=1;
    fo(i,2,316227)
    {
        if(!bz[i]) p[++p[0]]=i,mu[i]=-1;
        fo(j,1,p[0])
        {
            int k=p[j]*i;
            if(k>316227) break;
            if(i%p[j]==0)
            {
                mu[k]=0;bz[k]=1;
                break;
            }
            bz[k]=1;mu[k]=-mu[i];
        }
    }
    fo(i,1,316227) mu[i]+=mu[i-1];
    for(ll i=1,i1=1;i<=m;i++,i1=min(m,min((ll)sqrt(n/(n/i/i)),(ll)m/(m/i))))
    {
        ll t=0,nn=n/i/i;
        fo(j,1,m/i)
        {
            ll nm=nn/j;
            for(ll k=j+1,k1=j+1;k<=2*j;k++,k1=min(2*j,nm/(nm/k)))
            {
                t+=(nm/k)*(k1-k+1);
                if(k1==2*j) break;
                k=k1;if(nm/(k+1)==0) break;
            }   
        }
        ans+=(mu[i1]-mu[i-1])*(t);
        if(i1==m) break;
        i=i1;
    }
    ans-=n/2;
    printf("%lld",ans);
}