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; }
