这题真是过了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;
}