poj-2480(数论+积性函数性质)

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，让我们计算1到n的每个数i 与n的gcd(i,n)的和

题目分析:这道题目得用到一种数论函数（积性函数）

积性函数指对于所有互质的整数a和b有性质f(ab)=f(a)f(b)的 数论函数 。

百度一下我们可以知道很多常用的积性函数，就这道题目来说，我们只用知道gcd(i,k)是积性函数

即k=n*m   gcd(i,k)=gcd(i,n)*gcd(i,m)

∑gcd(i, N)这个函数也是积性函数。

我们对N质因子分解

N=p1^a1*p2^a2*p3^a3*****pn^an

这根据积性函数的性质:∑gcd(i, N)=∑gcd(i, p1^a1)*∑gcd(i, p2^a2)*∑gcd(i, p3^a3)********∑gcd(i, pn^an)

现在我们就是求每一个∑gcd(i, pk^ak)。

很明显gcd(i,pk^ak)只可能是pk^at次方(0<=at<=ak)

假设gcd(i,pk^ak)=pk^at 那么i的个数为phi[pk^ak/pk^at]

我们拿gcd(i,a)=b证明 

即i=k*b  a=p*b 使得gcd(k*b,p*b)=b, k与p互质的对数就是phi[a/b]

则我们的到来一个式子的答案：

∑gcd(i, pk^ak)=pk^0*phi[pk^(ak-0)]+pk^1*phi[pk^(ak-1)+.....+pk^ak*phi[pk^(ak-ak)]

我然后就写了代码，然后TLE了！

其实还可以变形:

对于phi[pk^ak],因为pk是质数的原因，那么能与不互质的一定是pk的倍数

则phi[pk^ak]=pk^ak-pk^(ak-1)

然后代入：

我们得到了

pk^ak*(1+ak*(1-1/pk)

再进一步化简：

pk^ak+ak*pk^ak-ak*p^(ak-1)

然后就可以做了.

代码如下:

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#define ll long long
using namespace std;

const int maxn=1e6+10;
int prime[maxn],cnt;
bool vis[maxn];
void get_prime()//打表
{
    memset (vis,true,sizeof (vis));
    vis[1]=false;
    for (int i=2;i<maxn;i++) {
        if (vis[i]) prime[cnt++]=i;
        for (int j=0;j<cnt&&i*prime[j]<maxn;j++) {
            vis[i*prime[j]]=false;
            if (i%prime[j]==0) break;
        }
    }
}
int factor[1005][2],fcnt;
void get_factor(int n)//分解质因子
{
    memset (factor,0,sizeof (factor));
    fcnt=0;
    for (int i=0;i<cnt&&prime[i]*prime[i]<=n;i++) {
        if (n%prime[i]==0) {
            factor[fcnt][0]=prime[i];
            while (n%prime[i]==0) {
                factor[fcnt][1]++;
                n=n/prime[i];
            }
            fcnt++;
        }
    }
    if (n!=1) {
        factor[fcnt][0]=n;
        factor[fcnt++][1]=1;
    }
}

int main()
{
    get_prime();
    int n;
    while (scanf("%d",&n)!=EOF) {
        get_factor(n);
        ll ans=1;
        for (int i=0;i<fcnt;i++) {
            ll res=pow(factor[i][0],factor[i][1]);//套入公式里面
            ans*=(res+factor[i][1]*res-factor[i][1]*res/factor[i][0]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}