pku1845 Sumdiv

题意：求a^b mod 9901 (0 <= A,B <= 50000000)

用到了二分求等比数列和

 

#include <iostream>
using namespace std;

#define MAXN 7100

int prim[MAXN],pcnt;
bool is_prim[MAXN];

void prepare(){
    __int64 i,j;
    pcnt=0;
    memset(is_prim,true,sizeof(is_prim));
    for(i=2;i<MAXN;i++){
        if(is_prim[i]){
            for(j=i;i*j<MAXN;j++)
                is_prim[i*j]=false;
            prim[pcnt++]=i;
        }
    }
}

__int64 exp_mod(__int64 A,__int64 B,__int64 M){
    __int64 d=1,bit[100],k=0,i;
    while(B){
        bit[k++]=B%2;
        B/=2;
    }//(bit[k],bit[k-1],,bit[1],bit[0])是B的二进制表示
    for(i=k-1;i>=0;i--){
        d=d*d%M;
        if(bit[i]==1)
            d=d*A%M;
    }
    return d;
}


__int64 sum(__int64 p,__int64 k){//1+p+p^2++p^(k-1)
    if(k==1)
        return 1;
    __int64 t=sum(p,k/2);
    if((k&1)==1)
        return ((1+exp_mod(p,k/2+1,9901))*t%9901+exp_mod(p,k/2,9901))%9901;
    else
        return (1+exp_mod(p,k/2,9901))*t%9901;
}



int main(){
    prepare();
    __int64 a,b,ans,i,j;
    while(scanf("%I64d%I64d",&a,&b)!=EOF){
        ans=1;
        //b%=9901;
        for(i=0;i<pcnt && prim[i]*prim[i]<a && a>1;i++){
            for(j=0;a%prim[i]==0;j++)
                a/=prim[i];
            if(j)
                ans=ans*sum(prim[i],j*b+1)%9901;
        }
        if(a>1 || i==pcnt)
            ans=ans*sum(a,b+1)%9901;
        printf("%I64d\n",ans);
    }
    return 0;
}