POJ2480 Longge's problem 欧拉函数应用

题目链接：POJ2480

题目大意：∑gcd(i, N) 1<=i <=N.，Given an integer N(1 < N < 2^31),

代码，思路借鉴：大佬的思路

在数论中的积性函数：对于正整数n的一个函数 f(n)，当中f(1)=1且当a,b互质，f(ab)=f(a)f(b)，在数论上就称它为积性函数。若某函数f(n)符合f(1)=1，且就算a,b不互质，f(ab)=f(a)f(b)，则称它为完全积性函数。

欧拉函数，gcd(n,k)(当k固定时)都是积性函数

且当i,j互素时，gcd(i*j,m)=gcd（i,m)*gcd(j,m),所以gcd(n,k)是积性函数

同时，积性函数的和也是积性函数

下文来源：http://lydws.blog.163.com/blog/static/22621105120152265175340/

大概只有我这种人不会做。。。

看了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；

以上这个公式建议大家手推一下，本人一开始没怎么搞懂，后来耐下心来推了一下公式就明白多了。

代码对原博主的代码做了一些修改，速度快了一倍。

总体来说应该比较好懂，打了一个素数表，对n进行分解成质因子的乘积形式，然后根据公式推就出来了。

AC代码

/*
POJ2480
2017年8月7日15:33:36 
手敲 
打了1e5范围内的素数表
速度快了一倍
AC 
*/ 
#include<stdio.h>
#include<math.h>
#include<string.h>
typedef long long ll;
const ll maxn=1e5+10;
ll flag[maxn];//标记 
ll prime[maxn] ;//存储素数 
ll cnt;//计数 
void isprime(){
	memset(flag,0,sizeof(flag));
	cnt=0;
	flag[1]=1;
	for(ll i=2;i<maxn;i++){
		/*如果当前i没有访问过，那么将其放进素数表*/ 
		if(!flag[i]){
			prime[cnt++]=i;
			/*从 i*i-maxn 范围内，所有该素数的倍数全部打上标记 */ 
			for(ll j=i*i;j<maxn;j+=i) flag[j]=1;
		} 
	}
} 

void solve(ll n){
	ll ans=n;
	for(ll i=0;prime[i]<=sqrt(n);i++){
		if(n%prime[i]==0){
			ll j=0;
			while(n%prime[i]==0) n/=prime[i],j++;
			ans/=prime[i];//考虑到j可能除不尽这个数，造成精度误差，所以对公式进行了变形。
			ans*=prime[i]+prime[i]*j-j; 
		}
	}
	if(n>1){
		ans/=n;
		ans*=n+n-1;
	}
	printf("%I64d\n",ans);
}


int main(){
	ll n;
	isprime();
	while(~scanf("%I64d",&n)){
		solve(n);
	}
	return 0;
}