Description
Longge is good at mathematics and he likes to think about hard mathematical problems which will be solved by some graceful algorithms. Now a problem comes: Given an integer N(1 < N < 2^31),you are to calculate ∑gcd(i, N) 1<=i <=N.
"Oh, I know, I know!" Longge shouts! But do you know? Please solve it.
"Oh, I know, I know!" Longge shouts! But do you know? Please solve it.
Input
Input contain several test case.
A number N per line.
A number N per line.
Output
For each N, output ,∑gcd(i, N) 1<=i <=N, a line
Sample Input
2 6
Sample Output
3 15
大概只有我这种人不会做。。。
看了lyd大神(和我的名字好像。。。)的题解,有些懂了:
大神链接。
主要考察一个叫积性函数的东西:f(x*y)=f(x)*f(y),中学数学经常见到。
积性函数还有个性质就是积性函数的和也是积性函数。
可以得到一个式子f(n)=f(p1^a1)*f(p2^a2)*f(p3^a3)*...*f(pk^ak)。不过目前我们不知道这个式子有什么用,6666666666666666。
再看本题,如果n和m互质,那么gcd(i,n*m)=gcd(i,n)*gcd(i,m)。∴gcd(i,n)是积性函数==>Σgcd(i,n)是积性函数。
设f(n)=Σgcd(i,n),问题转换为求出所有f(pi^ai)。
下面来求f(pi^ai):
首先明确,如果p是n的约数,那么满足gcd(i,n)==p的i的个数是Φ(n/p)。
证明:gcd(i,n)==p。设i=k*p,n=m*p。
则gcd(k,m)=1,也就是k和m互质,要求出满足条件的i的个数,就是求出i所对应的k的个数,即求m的欧拉函数,m=n/p,所以满足条件的i的个数就是Φ(m)=Φ(n/p)。
好了,用这种原则来求f(pi^ai),就是枚举pi^ai的每个约数(其实就是pi^t,t<=ai),然后求使gcd(i,pi^ai)==pi^t满足的i的个数,所以 个数*pi^t就是要求的答案的一部分。
公式:
f(pi^ai)=Φ(pi^ai)+pi*Φ(pi^(ai-1))+pi^2*Φ(pi^(ai-2))+...+pi^(ai-1)* Φ(pi)+pi^ai*Φ(1);
这里可以把求欧拉函数的部分化简,因为f(pi^ai)中只有一个约数pi,所以:
Φ(pi^ai)=pi^ai-pi^(ai-1)。
证明:小于pi^ai的正整数个数为p^ai - 1个;
其中,和pi^ai不互质的正整数有(pi*1,pi*2,...,pi*(pi^(ai-1)-1) )共计 pi^(ai-1)-1个。
∴
Φ(pi^ai)=pi^ai -1 -(pi^(ai-1)-1)=pi^ai-pi^(ai-1)。
然后整理得:f(pi^ai)=
pi^ai*(1+ai*(1-1/pi))。
所以,
f(n)=n*(1+a1*(1-1/p1))*(1+a2*(1-1/p2))*.…… =n* π(ai*pi+pi-ai)/pi;
本题解决!!!
写了两个小时,写完了,发现自己好弱,大神的一句话我都得反复琢磨还要问老师各种证明才能弄懂,差距很大,好大的压力。。。
代码如下:
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
ll N,ans;
void solve()
{
ans=N;
for(ll i=2;i*i<=N;i++)
if(N%i==0)
{
ll j=0;
while(N%i==0) N/=i,j++;
ans/=i;
ans*=i*j+i-j;
//ans*=1+j-j/i;//用上面的求法而不用这个是因为ans/i能整除而j/i不能整除,精度有问题。
}
if(N>1)//此时N的意义就是原来的N的一个大于sqrt(N)的质因子,再算一遍就好了。
{
ans/=N;
ans*=N+N-1;
}
printf("%lld\n",ans);
}
int main()
{
//freopen("a.in","r",stdin);
while(scanf("%lld",&N)!=EOF) solve();
return 0;
}