[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)|d2ij,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=L−−√

Ans=∑i=2m∑j=1i−1⌊L(i+j)i⌋[gcd(i,j)=1]

因为当

i=j, gcd(i,j)!=1，所以为了方便，将 j的范围设到i

Ans=∑i=2m∑j=1i⌊L(i+j)i⌋[gcd(i,j)=1]

设

f(d)=∑i=2m∑j=1i⌊L(i+j)i⌋[gcd(i,j)=d]

Ans=f(1)

LF同志的名言：看到GCD，就要信心满满上反演！

设

此时我们暂时加上

i=1的情况，最后再特判一下就好。

g(d)=∑k=1⌊md⌋f(kd)=∑i=1⌊md⌋∑j=1i⌊L(i+j)id2⌋

f(d)=∑k=1⌊md⌋μ(k)g(kd)=∑k=1⌊md⌋μ(k)∑i=1⌊mkd⌋∑j=1i⌊L(i+j)i(kd)2⌋

Ans=f(1)=∑k=1mμ(k)∑i=1⌊mk⌋∑j=1i⌊L(i+j)ik2⌋

前后都分块，复杂度是?
式子有点复杂，带了几个根号，不过它开了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;
}