Description
求有多少组正整数对 ((a, b)) 满足
- (a + b ≤ n)
- (a + b | ab)
(n ≤ 10^14)
Solution
这题有点绕啊
设 (gcd(a,b)=d),(a=a'd),(b=b'd)
对于第二个式子:((a'+b')d | a'b'd^2),所以 (a'+b' | a'b'd)
然后因为已经提出来了一个 (d) , 所以(gcd(a',b')=1) 即 (a'ot b')
所以 (gcd(a'+b',a')=gcd(a'+b',b')=gcd(a',b')=1)(模拟exgcd的过程)
意思就是 (a') 和 (b') 均不包含 (a'+b') 的任意因数,所以 (a'b') 也肯定不整除 (a'+b')
所以 (a'+b' | d)
那么有了这个时候,对于第一个式子:((a'+b')d leq n),而 (a'+b' | d) ,所以 (a'+b' leq sqrt n)
然后我们设 (t=a'+b'),从 (2) 到 (sqrt n) 枚举,每次算 (t) 固定后每个 (t) 的贡献
- 固定了 (t) 后,就要找 (d) 有多少种取值
设 (d=mt),(d) 有多少种取值就是 (m) 有多少种取值
看第一个式子,代进去,(t^2mleq n),所以只要 (mleq frac{n}{t^2}),就都可以
所以 (m) 共有 (lfloor frac{n}{t^2} floor) 种取值,所以 (d) 有 (lfloor frac{n}{t^2} floor) 种取值 - 再看固定了 (t) 后,把 (t) 分解成互质的 (a') 和 (b') 有多少种方案
我们要求有多少 (a') 和 (b') ,(a'+b'=t),(a'ot b')
因为(gcd(a',b')=1),所以(gcd(a',t-a')=1),根据更相减损术,(gcd(a',t)=1),(a')有多少取值,分解 (a') 和 (b') 就有多少方案
而 (a') 的取值方案就是 (phi (t))了
所以最后的答案就是(ans=sum_{i=2}^{sqrt n}lfloor frac{n}{i^2}
floorphi (i))
预处理后枚举求和
#include<bits/stdc++.h>
#define ll long long
#define db double
#define ld long double
const int MAXN=10000000+10;
int cnt,vis[MAXN],prime[MAXN];
ll phi[MAXN],n,res;
template<typename T> inline void read(T &x)
{
T data=0,w=1;
char ch=0;
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
x=data*w;
}
template<typename T> inline void write(T x,char c=' ')
{
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
if(c!=' ')putchar(c);
}
template<typename T> inline void chkmin(T &x,T y){x=(y<x?y:x);}
template<typename T> inline void chkmax(T &x,T y){x=(y>x?y:x);}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline void init()
{
memset(vis,1,sizeof(vis));
vis[0]=vis[1]=0;
phi[1]=1;
for(register int i=2;i<MAXN;++i)
{
if(vis[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for(register int j=1;j<=cnt&&i*prime[j]<MAXN;++j)
{
vis[i*prime[j]]=0;
if(i%prime[j])phi[i*prime[j]]=phi[i]*phi[prime[j]];
else
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
}
}
}
int main()
{
init();
read(n);
for(register ll i=2,limit=sqrt(n);i<=limit;++i)res+=(n/(i*i))*phi[i];
write(res,'
');
return 0;
}