POJ - 1845 Sumdiv【母函数 + 约数和 + 经典数论】

Sumdi

Time Limit: 1000MS Memory Limit: 30000K
Total Submissions: 25345 Accepted: 6280

Description

Consider two natural numbers A and B. Let S be the sum of all natural divisors of A^B. Determine S modulo 9901 (the rest of the division of S by 9901).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8.
The natural divisors of 8 are: 1,2,4,8. Their sum is 15.
15 modulo 9901 is 15 (that should be output).

题意： 让你求 $a ^ b$的约数和 对 $9901取模$

分析： 首先对a进行质因数分解，然后利用数学公式可以得到

$a = p_1 ^{k_1}·p_2 ^{k_2}·p_3 ^{k_3}·····p_n ^{k_n}$
$a ^b = p_1 ^{k_1·b}·p_2 ^{k_2·b}·p_3 ^{k_3·b}·····p_n ^{k_n·b}$

根据母函数的性质，我们可以得到$a^b$所有的约数为

$s = (1 + p_1^{1} + p_1^{2}+···+ p_1^{k1}) * (1 + p_2^{1} + p_2^{2}+···+ p_2^{k2}) *·····*(1 + p_n^{1} + p_n^{2}+···+ p_n^{kn})$

因为a的因子不会很多，一个个枚举计算即可，如第一个因子，主要是计算$\sum_{i = 0}^{k1} p_1^i$

我一开始wa这了，主要是因为看到等比数列直接用了等比数列的前n项和，这里是不对的因为存在$gcd(a,mod)!=1$,我们只能朴素的想了
这里可以观察数列的性质，二分进行计算即可
比如n为奇数的时候，比如5吧，假设为第一个因子，
——那么数列为：$1 + p1 +p1^2 + p1^3 +p1^4 + p1^5$
——我们可以将数列的前三项和后三项结合，
——在后三项中提取个 $p1^3,后三项变为p1^3 * (1 + p1 +p1^2)$
当n为偶数的时候同理，只不过先吧中间的数额外算即可

通过上面的我们就可以递归的来计算$\sum_{i = 0}^{k1} p_1^i$，从而得到求解

这里分解质因数默认都会啦

参考代码

#include<cstring>
#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cmath>


#define mod(x) ((x) % MOD)

using namespace std;

const int maxn = 1000000 + 100, MOD = 9901;

typedef long long ll;

bool isp[maxn];
ll p[maxn / 7];
int len;

void init() {
    isp[0] = isp[1] = true;
    for (int i = 2; i < maxn; i++) {
        if (!isp[i]) {
            p[len++] = i;
            for (int j = i + i; j < maxn; j += i) {
                isp[j] = true;
            }
        }
    }
}

ll a,b;

ll f[10010], s[10010];
int cnt;

ll qpow(ll a, ll b) {
    ll res = 1;
    while (b) {
        if (b & 1) res = mod(res * a);
        a = mod(a * a);
        b >>= 1;
    }
    return res;
}

ll sum (ll r, ll n) {
    if (n == 0) return 1;
    if (n & 1) return mod((qpow(r, (n + 1) / 2) + 1) * sum(r, n / 2));
    else return mod(qpow(r, n / 2) + mod((qpow(r, n / 2 + 1) + 1) * sum(r, n / 2 - 1)));
}

int main() {
    init();
    while (~scanf("%lld%lld", &a, &b)) {
        if (a == 0) {
            puts("0"); continue;
        }
        memset(f, 0, sizeof f);
        memset(s, 0, sizeof s);
        ll t = a;
        cnt = 0;
        for (int i = 0; i < len; i++) {
            if (t == 1) break;
            if (t % p[i] == 0) {
                f[cnt] = p[i];
                int tt = 0;
                while (t % p[i] == 0) tt++,t /= p[i];
                s[cnt++] = tt;
            }
        }
        if (t != 1) {
            f[cnt] = t;
            s[cnt++] = 1;
        }
        for (int i = 0; i < cnt; i++) s[i] *= b;
        ll res = 1;
        for (int i = 0; i < cnt; i++) {
            res = mod(res * 1ll * sum(f[i], s[i]));
        }
        printf("%lld\n", res);
    }
    return 0;
}