bjfu1109 最小公倍数和

这题真是过了n年才a。最早是在2010年北大培训比赛上看到的这题，当时我不会，竹教主也不会，但他记下来了，研究一段时间后就会了，还把这题加到我校oj上。过了这么多年，我上网搜，关于这个问题的解题报告还是没有，于是我花了几天时间做了出来，发布此解题报告。

题目是要求从1到n的所有数与n的最小公倍数的和，再除了n。第一眼看到数据范围，就知道是不能硬做(即从1到n依次算一遍)的。那么，为了找规律，得把公式化一化。

f(n) = [lcm(1,n)+lcm(2,n)+……+lcm(n,n)]/n
　　 = 1/gcd(1,n)+2/gcd(2,n)+……+n/gcd(n,n)

于是f(n)化成了一堆数的和。

显然，每一个gcd(k,n)都是能整除n的，按gcd(k,n)的值进行分类合并，就可以把f(n)整理成以下形式

f(n) = (……)/1 + (……)/2 +……+ (……)/d

这里d就是所有能整除n的整数。而(……)/1的分子为所有与n互质的数的和；(……)/2的分子为所有与n的最大公约数为2，也就是与n/2互质的数的和；依次类推，得到公式

有了这个公式，虽然是前进了一大步，但依然没法解题。因为这还是必须枚举n的所有约数，而这个复杂度依然为O(n)。

接下来就是比较核心的一步，构造。

f(n) = (g(n) + 1) / 2，求得g(n)就立即可得f(n)了。

这里构造g(n)的原因，是因为g(n)有积性性质。

即

对于任意gcd(m, n) = 1，有g(m * n) = g(m) * g(n)。证明略掉，读者自证。

有了积性性质就好办了，将n质因数分解成pi^ci相乘的形式，n = ∏(pi^ci)，则g(n) = g(∏(pi^ci)) = ∏(g(pi^ci))

g(pi^ci)很好计算。这里的d就是1, pi, pi^2, pi^3, ..., pi^ci，而phi(pi^ci)=(pi-1)*pi^(ci-1)，合起来就是个等比数列，可以推出来一个公式。也就是说可以O(1)时间算出g(pi^ci)。

如此一来，整个问题就可解了。

不过这题还是很变态的，因为数据就有50000组，仅分解质因数写得不好，都可能挂掉。最后我还是过了，代码如下：

/*
 * bjfu1109
 * Author    : ben
 */
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include <iostream>
#include <algorithm>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <functional>
#include <numeric>
#include <cctype>
using namespace std;
#ifdef ON_LOCAL_DEBUG
#else
#endif
typedef long long LL;
typedef int typec;

LL getPow(int a, int b) {
    LL res, temp;
    res = 1, temp = (LL) a;
    while (b) {
        if (b & 1) {
            res = res * temp;
        }
        b >>= 1;
        temp = temp * temp;
    }
    return res;
}

int get_int() {
    int res = 0, ch;
    while (!((ch = getchar()) >= '0' && ch <= '9')) {
        if (ch == EOF)
            return -1;
    }
    res = ch - '0';
    while ((ch = getchar()) >= '0' && ch <= '9')
        res = res * 10 + (ch - '0');
    return res;
}

const int N = 100000;
bool isPrime[N + 3];//多用两个元素以免判断边界
vector<int> pt;
void init_prime_table() {
    memset(isPrime, true, sizeof(isPrime));
    int p = 2, q, del;
    double temp;
    while (p <= N) {
        while (!isPrime[p]) {        p++;        }
        if (p > N) {//已经结束
            break;        }
        temp = (double) p;
        temp *= p;
        if (temp > N)
            break;
        while (temp <= N) {
            del = (int) temp;            isPrime[del] = false;
            temp *= p;        }
        q = p + 1;
        while (q < N) {
            while (!isPrime[q]) {    q++;    }
            if (q >= N) { break;}
            temp = (double) p;
            temp *= q;
            if (temp > N)    break;
            while (temp <= N) {
                del = (int) temp;
                isPrime[del] = false;
                temp *= p;
            }
            q++;
        }
        p++;
    }
    for (int i = 2; i <= N; i++) {
        if (isPrime[i]) {
            pt.push_back(i);
        }
    }
}


/**
 * p为素数表，至少应该包括sqrt(N)以内的所有素数，f保存结果
 * f[i].first表示N的一个素因数,f[i].second为这个素因子的个数
 * typec可以是int、long、long long等
 */
typedef vector<pair<typec, int> > FactorList;
void get_prime_factor(const typec &N_, FactorList &f, const vector<int> &p) {
    int i, t, n, pl = p.size();
    typec N = N_;
    f.clear();
    for(i = 0; i < pl; i++) {
        t = p[i];
        if (t * t > N_) {
            break;
        }
        if(N % t == 0) {
            n = 0;
            while(N % t == 0) {
                n++;    N /= t;
            }
            f.push_back(make_pair(t, n));
        }
        if(N == 1) {    break;    }
    }
    if(N > 1) {
        f.push_back(make_pair(N, 1));
    }
}

inline LL g(LL p, int c) {
    LL ret = getPow(p, 2 * c) - 1;
    ret = ret / (p + 1) * p + 1;
    return ret;
}

FactorList fl;
LL g(int n) {
    get_prime_factor(n, fl, pt);
    int len = fl.size();
    LL ret = 1;
    for (int i = 0; i < len; i++) {
        ret *= g(fl[i].first, fl[i].second);
    }
    return ret;
}

int main() {
#ifdef ON_LOCAL_DEBUG
    freopen("data.in", "r", stdin);
#endif
    LL ans, N;
    init_prime_table();
    while ((N = get_int()) > 0) {
        ans = g(N) + 1;
        ans /= 2;
        printf("%I64d\n", ans);
    }
    return 0;
}