2480 Longge's problem 欧拉函数变形 求,∑gcd(i, N)

Longge's problem

Time Limit: 1000MS		Memory Limit: 65536K
Total Submissions: 3822		Accepted: 1128

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

先给出一个证明吧：首先假设n有一个约数d，那么怎样计算出1～n中最大公约数为d的个数呢？

很显然，这个个数实质上是等于fin(n/d)(其中先用fin代表欧拉函数），想到这里的话，基本上就确定了策略，

我们先枚举出n的所有约数，然后求出每一个的欧拉函数，然后d*fin(n/d)相加后的结果即为所求，

但是枚举出n的所有约数，这是一个很难的问题，首先那些因数怎么求呢？不过，题目是求和，并不是一个一个地求，

于是我们把欧拉函数的公式套上，可得d*n/d*(1-1/p1)*(1-1/p2)*...*(1-1/pm),

化简得到n*(1-1/p1)(1-1/p2)*..(1-1/pm),那么所有的项其实都可以提一个n出来，

于是关键是求最后一部分的和，首先最后一部分其实是n/d的因数分解所得出的p1,p2,...pm,

那么我们考虑假设n的因数分解是p1^r1*p2^r2*...*pn^rn,那么n的因子d其实都可以表示成

p1^k1*p2^k2*...pn^kn,其中0<=ki<=ri,那么如果ki不为ri的话，n/d这个数中必然含有p1这个素因子，

否则的话就不含p1这个素因子，到了这里，利用排列组合的知识我们可以写出一个最后一部分的和的公式了:

(1+r1*(1-1/p1))*(1+r2*(1-1/p2))*...(1+rn*(1-1/pn));其实是这样的,当不包含p1这个素因子时，第一项选1，

然后若包含p1这项因子时，那么n/d中的p1的指数可以有r1中情况，所以第一项选最后一个r1*(1-1/p1)。

然后得出了最后的公式n*((1+r1*(1-1/p1))*(1+r2*(1-1/p2))*...(1+rn*(1-1/pn)));现在只需要进行因数分解，

这个问题可以在大概的O(sqrt(n))的时间求出。（这是看别人的解法才知道的，自己只推出前半部分）

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
int is[100000],prime[100000],pl=0;
int main()
{
    for(int i=2;i<100000;i++)
    {
        if(is[i]==0) prime[pl++]=i;//0表示素数 
        for(int j=0;j<pl&&prime[j]*i<100000;j++)
        {
            is[prime[j]*i]=1;
            if(i%prime[j]==0) break;
        }
    }
    long long  n;
    while(cin>>n)
    {
        long long sum=n;
        for(int i=0;i<pl&&prime[i]<=sqrt(n*1.0);i++)
        {
            if(n%prime[i]==0)
            {
                int r=0;
                while(n%prime[i]==0)
                {
                    r++;
                    n/=prime[i];
                }
                sum=sum*(1+r)-sum/prime[i]*r;
            }
        }
        if(n>1)
        {
            sum=sum*2-sum/n*1;
        }
        cout<<sum<<endl;
    }
    return 0;
}