poj2480

题意:给出n(2^31),求∑gcd(i, n) 1<=i <=n。

分析:首先所有与n互质的数字x都满足gcd(x,n)=1。我们先计算这种等于1的情况,恰好是n的欧拉函数phi(n)。我们可以将上述情况视为跟n最大公约数为1的情况,现在我们将其推广至最大公约数为p的情况。对于对于所有满足gcd(x,n)=p(p为常数)的x,他们与n拥有相同的gcd,那么他们同时除以p之后,就会变得和n/p互质,这种数字x有phi(n/p)个,这些x的和为p*phi(n/p)个。所以我们要计算∑gcd(i, n) 1<=i <=n,只需要根据gcd的值不同,分类进行计算即可,总结成公式:∑p*phi(n/p)(p是n的约数)

对于这个公式我们有一个快速计算的方法,我们先想,若要计算∑p,我们可以用公式(x1^0+x1^1+...+x1^c1)*(x2^0+x2^1+...+x2^c2)*...*(xm^0+xm^1+...+xm^cm)。xi是n的质因数,ci是n中含有xi的个数。同理,根据欧拉函数的公式,每个欧拉函数值也和各个质因子有关,求∑phi(n/p)=∑phi(p)=(1 + x1^0 * (x1 - 1) + x1^1 * (x1- 1) + ...+x1^(c1-1) * (x1- 1))*...*(...)。然后由于二者都是质因子幂相乘的形式,所以可以把两者综合起来,对于每个括号里的每一项,都是约数因式和欧拉函数因式的乘积的形式,即p的因式要乘以phi(n/p)的因式,两因式的质因子指数是互补的,最终公式变为(x1^(c1-1)*(c1-1) + x1^c1)*...。

View Code
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;

int n;

bool is[(1 << 16) + 5];
int prm[(1 << 16) + 5];

int getprm(int n)
{
    int i, j, k = 0;
    int s, e = (int) (sqrt(0.0 + n) + 1);
    memset(is, 1, sizeof(is));
    prm[k++] = 2;
    is[0] = is[1] = 0;
    for (i = 4; i < n; i += 2)
        is[i] = 0;
    for (i = 3; i < e; i += 2)
        if (is[i])
        {
            prm[k++] = i;
            for (s = i * 2, j = i * i; j < n; j += s)
                is[j] = 0;
        }
    for (; i < n; i += 2)
        if (is[i])
            prm[k++] = i;
    return k;
}

long long power(int p, int cnt)
{
    if (cnt < 0)
        return 1;
    long long ret = 1;
    long long temp = p;
    while (cnt)
    {
        if (cnt & 1)
            ret *= temp;
        cnt >>= 1;
        temp = temp * temp;
    }
    return ret;
}

long long cal(int p, int cnt)
{
    return power(p, cnt - 1) * cnt * (p - 1) + power(p, cnt);
}

int main()
{
    //freopen("t.txt", "r", stdin);
    int m = getprm(1 << 16);
    long long ans;
    while (~scanf("%d", &n))
    {
        ans = 1;
        for (int i = 0; i < m && prm[i] * prm[i] <= n; i++)
        {
            int cnt = 0;
            while (n % prm[i] == 0)
            {
                n /= prm[i];
                cnt++;
            }
            ans *= cal(prm[i], cnt);
        }
        if (n != 1)
            ans *= cal(n, 1);
        printf("%lld\n", ans);
    }
    return 0;
}

 

原文地址:https://www.cnblogs.com/rainydays/p/2745672.html