POJ-1485 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

解题思路

把 $A$ 分解质因数，可以得到 $p_1^{c_1} \times p_2^{c_2} \cdots \times p_n^{c_n} $

$A^B$ 可表示为 $p_1^{B \times c_1} \times p_2^{B \times c_2} \cdots \times p_n^{B \times c_n} $

那么 $A^B$ 的约数可以表示为 $p_1^{k_1} \times p_2^{k_2} \times \cdots \times p_n^{k_n} $，其中 $ 0 \leq k_i \leq B \times c_i $

根据乘法分配律，约数之和为：$(1 + p_1 + p_1^2 + \cdots p_1^{B \times c_1}) \times (1+p_2+p_2^2+ \cdots +p_2^{B \times c_2}) \times \cdots \times (1+p_n+p_n^2+ \cdots +p_n^{B \times c_n}) $

我们令 $sum(p,c) = 1 + p + p^2 + \cdots + p^c $，可以通过分治的方式解决

若 $p$ 是奇数
$ sum(p,c) = (1 + p + \cdots + p ^ {\frac{c-1}{2}}) + (p^{\frac{c+1}{2}} + \cdots + p^c) = (1 + p^{\frac{c+1}{2}}) \times sum(p, \frac{c}{2} – 1) $

若 $p$ 是偶数
$ sum(p,c – (1 + p^{\frac{c}{2}}) \times sum(p, \frac{c}{2} – 1) + p^c $

#include 
#include 
#include 
#include 
#include 

using namespace std;
const int MOD = 9901;
const int MAXN = 1000;

inline long long calc(int a, int b){
    long long tmp = 1;
    long long base = a;

    while(b){
        if(b & 1){
            tmp *= base;
            tmp %= MOD;
        }
        base *= base;
        base %= MOD; 

        b >>= 1;
    }

    return tmp;
}

inline long long sum(int p, int c){
    if(c == 1)
        return p + 1;
    else if(c == 0)
        return 1;

    if(c % 2 == 0){
        return ((1 + calc(p, c/2)) * sum(p, c/2 - 1) + calc(p, c)) % MOD;
    }else{
        return ((1 + calc(p, (c+1)/2)) * sum(p, (c-1)/2)) % MOD;    
    }
} 

long long prime[MAXN][2], cnt;
int a, b;
long long ans = 1;

int main(){
    scanf("%d%d", &a, &b);

    for(register int i=2; a!=1 ; i++){
        if(a % i == 0){
            prime[cnt][0] = i;
            while(a % i == 0){
                prime[cnt][1]++;
                a /= i;
            }
            cnt++;
        }
    }

    for(register int i=0; i

[题解] POJ-1845 Sumdiv

BLOG