poj1845Sumdiv+约数和定理

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的所有约数和；
解法：
$a=p_1^{n1}p_2^{n2}...p_n^{nm}$
两个定理：
一。约数的个数=$(n1+1)*(n2+1)*...*(nm+1)$
二。约数和=$(p_1^0+p_1^1+...+p_1^{n1})*(p_2^0+p_2^1+...+p_2^{n1})*...*(p_n^0+p_n^1+...+p_n^{nm})$

$a^b=p_1^{n1*b}p_2^{n2*b}...p_n^{nm*b}$
直接套公式就好。
$p_1^0+p_1^1+...+p_1^{n1}$就是等比数列求和$\frac {a_0(1-q^n)} {1-q}$
但是，计算一下中间结果就算是除法取膜先膜（（q-1)*mod），如果q是一个很大的质数，中间也可能溢出了。（wa5,6次）

所以要用递归二分来求
对于$1+a^1+a^1+....a^n$

当n为奇数时，可以化简为$(1+a^{\frac{n}{2}+1})*(1+a^1+...+a^{\frac{n}{2}})$

当n为偶数时，可以化简为$(1+a^{\frac{n}{2}+1})*(1+a^1+...+a^{\frac{n}{2}-1})+a^{\frac{n}{2}}$

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<map>
#include<cstring>
using namespace std;
#define LL long long
#define MOD 9901

int prime[10005];
bool visit[10005];
int cnt;
void getprime(){
    cnt=0;
    memset(visit,false,sizeof(visit));
    for(int i=2;i<=10005;i++){
        if(visit[i]==false){
            prime[cnt++]=i;
            for(int j=i*i;j<=10005;j+=i) visit[j]=true;
        }
    }
}
map<int,int>Q;
map<int,int>::iterator it;

LL quick_pow(LL a,LL n,LL mod){
    LL ans=1;
    while(n){
        if(n&1) ans=ans*a%mod;
        a=a*a%mod;
        n>>=1;
    }
    return ans;
}
LL sum(LL p,LL n){  //计算1+p+p^2+````+p^n
    if(p==0) return 0;
    if(n==0) return 1;
    if(n&1) return (1+quick_pow(p,n/2+1,MOD))*sum(p,n/2)%MOD;
    else return ((1+quick_pow(p,n/2+1,MOD))*sum(p,n/2-1)%MOD+quick_pow(p,n/2,MOD))%MOD;
}
int main(){
    getprime();
    LL A,B;
    while(cin>>A>>B){
        Q.clear();
        for(int i=0;i<cnt&&prime[i]*prime[i]<=A;i++){
            if(A%prime[i]==0){
                int cntp=0;
                while(A%prime[i]==0){
                    cntp++;
                    A/=prime[i];
                }
                Q[prime[i]]=cntp;
            }
        }
        if(A>1){
            Q[A]=1;
        }
        //for(it=Q.begin();it!=Q.end();it++) cout<<it->first<<" "<<it->second<<endl;
        LL ans=1;
        for(it=Q.begin();it!=Q.end();it++){
            LL temp=(LL)it->first;//p
            LL temp2=((LL)(it->second))*B;//n+1
            ans=ans*(sum(temp,temp2)%MOD)%MOD;
        }
        cout<<ans<<endl;
    }
    return 0;
}