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=pn11pn22...pnmn
两个定理：
一。约数的个数=(n1+1)∗(n2+1)∗...∗(nm+1)
二。约数和=(p01+p11+...+pn11)∗(p02+p12+...+pn12)∗...∗(p0n+p1n+...+pnmn)

ab=pn1∗b1pn2∗b2...pnm∗bn
直接套公式就好。
p01+p11+...+pn11就是等比数列求和a0(1−qn)1−q
但是，计算一下中间结果就算是除法取膜先膜（（q-1)*mod），如果q是一个很大的质数，中间也可能溢出了。（wa5,6次）

所以要用递归二分来求
对于1+a1+a1+....an

当n为奇数时，可以化简为(1+an2+1)∗(1+a1+...+an2)

当n为偶数时，可以化简为(1+an2+1)∗(1+a1+...+an2−1)+an2

#include<cstdio>#include<algorithm>#include<iostream>#include<map>#include<cstring>using namespace std;#define LL long long#define MOD 9901int 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;}