（poj 2480 Longge's problem）<欧拉函数>

题目

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.

Input

Input contain several test case.
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

大致意思是求一个数的∑gcd（i,n），i∈[1,n)

题解

1. 一个很显然的算法是：直接暴力枚举每一个gcd并求和。当然同样显然的是，它会TLE
2. 正确思路：∑gcd是一个积性函数，因此利用唯一分解定理，将n分为pi^ai，求出每一个的∑gcd，乘起来即可得到答案。

以下是一堆证明：

• 证明公式：∑gcd(i,n)=∑x*φ(N/x)
  因为在1···N中，若gcd(i，N) = x，则gcd(i/x , N/x)=1，因此在1···N中，gcd(i，N) = x 的i的个数的等于φ(N / x)。

• 证明∑gcd是积性函数
  由上一个证明知
  ∑gcd(1,n)=∑k*φ(n/k)
  设F(x)=x , G(x)=φ(x) , H(x)=∑gcd(i,x)
  ∴H(x)=∑F(x)*G(x)
  及∑k*φ(n/k)为F(x)和G(x)的狄利克雷卷积
  积性函数的狄利克雷卷积依然是狄利克雷卷积
  ∴H(x)是积性函数

• 由∑gcd是积性函数和唯一分解定理
  ∴H(n)=H(p1^a1)*H(p2^a2) *H(p3^a3)……H(n^an)，Pi是素数
  ∴只需要算出每一个H(pi^ai)，相乘即得答案

• 计算H(pi^a1)的方法：
  已经证明，如果p是n的约数，那么满足gcd(i,n)==p的i的个数是Φ(n/p)
  以及∑gcd(i,n)=∑x*φ(N/x)
  ∴f(pi^ai)=φ(pi^ai)+pi*φ(pi^(ai-1))+pi^2*φ(pi^(ai-2))+…+pi^(ai-1)* φ(pi)+pi^ai*φ(1)

• 计算φ(pi^ai)的方法：
  小于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)
  ∴H(pi^ai)
  =pi^(ai-1)(pi-1) + pi*pi^(ai-2)(pi-1)….+pi^ai
  = pi^ai*(1+ai*(1-1/pi))

• 化简：
  H(n)
  = p1^a1*p2^a2…pr^ar(1+a1*(1-1/p1))(1+a2(1-1/p2))*…
  = n*(1+a1*(1-1/p1))(1+a2(1-1/p2))*…
  = n*∏(ai*pi+pi-ai)/pi

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long LL;

LL n;

int main(){
    while(scanf("%d",&n)!=EOF){
        LL ans=n;
        for(LL pi=2;pi*pi<=n;++pi){
            if(n%pi==0){
                LL ai=0;    
                while(n%pi==0){
                    n/=pi;
                    ai++;
                }
                ans=ans/pi*(pi*ai+pi-ai);
            }
        }
        if(n>1){
            ans=ans/n*(n+n-1);
        }
        printf("%lld\n",ans);
    }
    return 0;
}