[poj1845] Sumdiv（数论，质因数分解，约数和，快速幂）

Sumdiv

Time Limit: 1000MS Memory Limit: 30000K

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).

Source

Romania OI 2002

---

解题思路

题意：求$A^B$的所有约数和模9901

首先我们需要知道约数和公式：

由唯一分解定理，对$\forall A \in Z$，可以将$A$写作：$A = p_1^{k_1}p_2^{k_2}\cdots p_n^{k_n}$，其中$p_i$是质数，那么$A$的约数和为$\prod\limits_{i=1}^n\sum\limits_{j=0}^{k_i} p_i^j$

这道题中，可以先把$A$分解质因数，假设$A = p_1^{k_1}p_2^{k_2}\cdots p_n^{k_n}$，那么$A^B= p_1^{Bk_1}p_2^{Bk_2}\cdots p_n^{Bk_n}$，所以答案就是$\prod\sum\limits_{j=0}^{Bk_i} p_i^j\mod9901$
所以现在的问题变成了：
1. 如何分解质因数
将$A$不断除2直至除尽，然后不断除3直至除尽……由于是从小到大不断相
除，可以保证$A$不会除到一个合数（因为$A$一定已经被这个合数的最小质因子除尽了），除的同时记录$A$的因子及其出现次数即可

int fac = 2;
while(fac * fac <= A){
    while(A % fac == 0){
        if(p[p[0]] != fac)  p[++p[0]] = fac;
        k[p[0]]++;
        A /= fac;
    }
    if(fac == 2)    fac++;
    else    fac += 2;
}
if(A != 1)  p[++p[0]] = A, k[p[0]] = 1;

2.如何快速地求出$\sum\limits_{j=0}^{n} p^j \mod 9901$，即$(p^0+p^1+\cdots+p^n)\mod 9901$

• 法一：
  
  • 如果$n$是奇数， $$(p^0+p^1+...+p^{n-1})+(p^n+p^{n+1}+...+p^{2n-1}) = (p^0+p^1+p^2+...+p^{n-1}) * (1+p^n)$$
  • 如果$n$是偶数， $$(p^0+p^1+...+p^{n-1})+p^n+(p^{n+1}+...+p^{2n}) = p^n+(p^0+p^1+p^2+...+p^{n-1}) * (1+p^{n+1})$$
  • 求和部分递归即可，$p^n$和$p^{n+1}$快速幂即可
    这样就可以$O(log_2N)$地计算了
• 法二：
  当然，也可以直接用等比数列求和公式，$(p^0+p^1+\cdots+p^n)=\frac{1-p^{n+1}}{1-p}$，但是由于要取模，所以注意乘逆元，另外，又因为$9901$是一个小质数，用快速幂时还需要特判

---

Code

#include<cstdio>
#include<cstring>

using namespace std;
typedef long long LL;

const int MOD = 9901;
int A, B, p[10001], k[10001], ans;

LL fPow(LL bs, LL idx){
    LL res = 1;
    bs %= MOD;
    while(idx){
        if(idx & 1) (res *= bs) %= MOD;
        (bs *= bs) %= MOD;
        idx >>= 1;
    }
    return res;
}

LL cal(int p, LL k){
    if(k == 0)  return 1ll;
    if(k & 1)   return cal(p, k/2) * (1ll + fPow(p, k/2+1)) % MOD;
    else    return (fPow(p, k/2) + cal(p, k/2-1) * (1ll + fPow(p, k/2+1)) % MOD) % MOD;
}

int main(){
    while(scanf("%d%d", &A, &B) != EOF){
        ans = 1;
        p[0] = 0;
        int fac = 2;
        while(fac * fac <= A){
            while(A % fac == 0){
                if(p[p[0]] != fac)  p[++p[0]] = fac;
                k[p[0]]++;
                A /= fac;
            }
            if(fac == 2)    fac++;
            else    fac += 2;
        }
        if(A != 1)  p[++p[0]] = A, k[p[0]] = 1;
        for(int i = 1; i <= p[0]; i++)
            (ans *= cal(p[i], 1ll * k[i] * B)) %= MOD;
        printf("%d\n", ans);
    }
    return 0;
}