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互质的数的和;依次类推,得到公式

bjfu1109 最小公倍数和

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

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

bjfu1109 最小公倍数和

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 = , temp = (LL) a;
while (b) {
if (b & ) {
res = res * temp;
}
b >>= ;
temp = temp * temp;
}
return res;
} int get_int() {
int res = , ch;
while (!((ch = getchar()) >= '' && ch <= '')) {
if (ch == EOF)
return -;
}
res = ch - '';
while ((ch = getchar()) >= '' && ch <= '')
res = res * + (ch - '');
return res;
} const int N = ;
bool isPrime[N + ];//多用两个元素以免判断边界
vector<int> pt;
void init_prime_table() {
memset(isPrime, true, sizeof(isPrime));
int p = , 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 + ;
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 = ; 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 = ; i < pl; i++) {
t = p[i];
if (t * t > N_) {
break;
}
if(N % t == ) {
n = ;
while(N % t == ) {
n++; N /= t;
}
f.push_back(make_pair(t, n));
}
if(N == ) { break; }
}
if(N > ) {
f.push_back(make_pair(N, ));
}
} inline LL g(LL p, int c) {
LL ret = getPow(p, * c) - ;
ret = ret / (p + ) * p + ;
return ret;
} FactorList fl;
LL g(int n) {
get_prime_factor(n, fl, pt);
int len = fl.size();
LL ret = ;
for (int i = ; 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()) > ) {
ans = g(N) + ;
ans /= ;
printf("%I64d\n", ans);
}
return ;
}
上一篇:azkaban使用--传入动态参数


下一篇:[转]Asp.net MVC 中Ajax的使用