【51nod1223】【分数等式的数量】【莫比乌斯反演】

最新推荐文章于 2022-04-08 09:24:32 发布

原创最新推荐文章于 2022-04-08 09:24:32 发布 · 451 阅读

1 ·

CC 4.0 BY-SA版权

数论同时被 2 个专栏收录

56 篇文章

订阅专栏

51nod

12 篇文章

订阅专栏

本文介绍了一种解决特定分数等式的算法。等式形式为1/X+1/Y=1/N，目标是在X<Y≤L的条件下找出所有可能的解。通过数学推导和优化，文章提供了一个高效的解决方案，并附带了C++实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目大意

有这样一个分数等式：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。

解题思路

$Ans=\sum_{x=1}^N\sum_{y=x+1}^N[(x+y)|(xy)]$

令d=(x,y),x’=x/d,y’=y/d。

$[(dx'+dy')|(d^2x'y')]=[(x'+y')|(dx'y')]$

由于(x’,y’)==1,(x’+y’,x’)==1,(x’+y’,y’)==1,所以(x’+y’)|d。

由于dx<=N,d<=N/d,又(x+y)|d,d取值有N/y/(x+y)个。

$Ans=\sum_{d=1}^N\sum_{x=1}^{N/d}\sum_{y=x+1}^{N/d}[(x,y)==1](N/(dy)/(dx+dy))$

$g[d]=\sum_{x=1}^{N/d}\sum_{y=x+1}^{N/d}[(x,y)==1](N/(dy)/(dx+dy))$

$f[d]=\sum_{x=1}^{N/d}\sum_{y=x+1}^{N/d}(N/(dy)/(dx+dy))$

$f[d]=\sum_{i=1}^{N/d}g[di]$

$g[d]=\sum_{i=1}^{N/d}f[di]mu[i]$

令d=bi。

$Ans=g[1]=\sum_{d=1}^Nmu[d]\sum_{x=1}^{N/d}\sum_{y=x+1}^{N/d}(N/(dy)/(dx+dy))$

$Ans=\sum_{d=1}^Nmu[d]\sum_{x=1}^{N/d}\sum_{y=x+1}^{N/d}(N/d/d/y/(x+y))$

因为N/d/d,所以d<=sqrt(n),因为N/d/d/y/(x+y),所以y<=sqrt(n)/d。

$Ans=\sum_{d=1}^{\sqrt n}mu[d]\sum_{y=1}^{\sqrt N/d}\sum_{x=1}^{y-1}(N/d/d/y/(x+y))$

枚举前两个循环，最后一个分块即可，具体请看代码。

code

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LF double
#define LL long long
#define Min(a,b) ((a<b)?a:b)
#define Max(a,b) ((a>b)?a:b)
#define Fo(i,j,k) for(int i=j;i<=k;i++)
#define Fd(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
int const Mxn=1e6;
LL N;int Tag[Mxn+9],Mu[Mxn+9],Pri[Mxn+9];
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w ",stdout);
    scanf("%lld",&N);int NN=sqrt(N);
    Mu[1]=1;
    Fo(i,2,NN){
        if(!Tag[i])Pri[++Pri[0]]=i,Mu[i]=-1;
        Fo(j,1,Pri[0]){
            if(i*Pri[j]>NN)break;
            Tag[i*Pri[j]]=1;
            Mu[i*Pri[j]]=-Mu[i];
            if(i%Pri[j]==0){Mu[i*Pri[j]]=0;break;}
        }
    }
    LL Ans=0;LL Nd,Y;
    Fo(d,1,NN)Fo(y,2,NN/d){
        Nd=N/d/d/y;Y=Min(y-1,Nd-y);
        for(int i=1,j;i<=Y;i=j+1){
            if(!(Nd/(i+y)))
                break;
            j=Nd/(Nd/(i+y))-y;
            j=(j>Y)?Y:j;
            Ans+=Mu[d]*(Nd/(i+y))*(j-i+1);
        }
    }
    printf("%lld",Ans);
    return 0;
}