bjfu1109 最小公倍數和

標籤：style   blog   http   ar   io   color   os   sp   for   

這題真是過了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#endiftypedef 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;}

 

bjfu1109 最小公倍數和