POJ1825 Sumdiv

问题描述
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).

翻译
给A，B，求A^B的因子和 mod 9901

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

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

样例输入
2 3

样例输出
15

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

来源
Romania OI 2002

代码

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#define LL long long
#define maxn 9999
#define m 9901
using namespace std;

bool isprime[maxn+100];
int prime[maxn+100],primesize;

void getlist(int listsize = maxn){
    memset(isprime,1,sizeof(isprime));
    isprime[1] = false;
    for(int i = 2;i <= listsize;i++){
        if(isprime[i]) prime[++primesize] = i;
        for(int j = 1;j <= primesize && i*prime[j] <= listsize;j++){
            isprime[i*prime[j]] = false;
            if(i%prime[j] == 0) break;
        }
    }
}

LL ksm(LL a,LL k){
    if(!k) return 1;
    LL b = 1;
    while(k){
        if(k&1) b = b*a %m;
        a = a*a %m;
        k >>= 1;
    }
    return b;
}

LL sum(LL p,LL n){
    if(n == 0) return 1;
    LL half = sum(p,(n-1)/2) %m;
    LL step = ksm(p,(n+1)/2) %m;
    LL ans;
    if(n % 2 == 0){
        ans = half %m + (half %m)*(step %m) + ksm(p,n);
    }else{
        ans = half %m + (half %m)*(step %m);
    }
    return ans %m;
}

LL solve(LL A,LL B){

    getlist();

    LL num,ans = 1;

    for(LL i = 1;prime[i] <= A && A != 1 && i <= primesize;i++){
        if(A % prime[i] == 0)
        {
            num = 0;
            while(A % prime[i] == 0){
                A /= prime[i];
                num++;
            }
            ans *= sum(prime[i],num*B) %m;
            ans %= m;
        }
    }
    if(A > 1)
    {
        ans *= sum(A,B) % m;
        ans %= m;
    }
    return ans;
}

int main(){
    LL A,B;
    scanf("%lld%lld",&A,&B);
    printf("%lld",solve(A,B));
    return 0;
}

待续