Longge's problem


Time Limit: 1000MS

 

Memory Limit: 65536K

Total Submissions: 8691

 

Accepted: 2909


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


Source


​​POJ Contest​​,Author:Mathematica@ZSUpu





        正如上一篇文章说的，这题就是求∑gcd(i, N)。


        求这个正是很好很巧妙的利用了函数的积性。首先设g(x)=gcd(x,N)，显然g(x)是具有积性的，那么根据定理，g(x)的和f(N)也是具有积性的。所以如果对于一个数字n，有n=p1^a1*p2^a2*p3^a3……pi^ai，那么f(n)=f(p1^a1)*f(p2^a2)*f(p3^a3)*……f(pi^ai)。现在的关键就是如何求每一个的f(pi^ai)。根据上一篇文章讲的，∑gcd(i, N)=∑dφ(N/d)。那么可以写成f(pi^ai)=1*φ(pi^ai)+pi*φ(pi^(ai-1))+pi^2*φ(pi^(ai-2))+……+pi^ai*φ(1)。然后再根据φ的积性，φ(pi^ai)=(pi-1)*pi^(ai-1)。故有f(pi^ai)=ai*(pi-1)*pi^(ai-1)+pi^ai=pi^ai*(1+ai(1-1/pi))。那么合起来f(n)=p1^a1*(1+a1(1-1/p1))*p2^a2*(1+a2(1-1/p2))*……*pi^ai*(1+ai(1-1/pi))=n*∏(1+ai(1-1/pi))。


        惊喜+意外，最后结果同样也是只要分解质因数即可，而且恰好结果里面合成回了n。具体见代码： 


#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#define ULL long long
#define N 101000
using namespace std;

ULL prime[N],tot,n,m,t;
bool isprime[N];

void getlist(int listsize)
{
    memset(isprime,1,sizeof(isprime));
    isprime[1]=false; tot=0;
    for(int i=2;i<listsize;i++)
    {
        if(isprime[i]) prime[++tot]=i;
         for(int j=1;j<=tot&&i*prime[j]<listsize;j++)
         {
            isprime[i*prime[j]]=false;
            if(i%prime[j]==0)break;
        }
    }
}

int main()
{
    getlist(N);
    while(~scanf("%I64u",&n))
    {
        ULL y=n;
        for(int i=1;prime[i]*prime[i]<=n;i++)
        {
            if (n%prime[i]==0)
            {
                t=0;
                while(n%prime[i]==0)
                    t++,n/=prime[i];
                y=y+y*t*(prime[i]-1)/prime[i];
            }
        }
        if (n>1) y=y+y*(n-1)/n;
        printf("%I64u\n",y);
    }
    return 0;
}