poj 1845 Sumdiv(分治+快速幂)

【题目大意】：求A^B的约数和mod 9901.(0 <= A,B <= 50000000)

【解题思路】：

A^B的约数和记为S，则S=(1+p1+p1^2+...+p1^q1)*(1+p2+p2^2+....+p2^q2)*....*(1+pn+pn^2+...pn^qn);

对于上述的式子，pi代表质因子，qi代表个数。

根据二分法分治，我们又可以知道：

1+p+p^2+...+p^n=

当n为奇数，(1+p+p^2+...+p^((n+1)/2))*(1+p^((n+1)/2))

当n为偶数，(1+p+p^2+...+p^((n+1)/2)+p^((n+1)/2)*(1+p*(1+p+p^2+...+p^((n+1)/2)))) 

因此，我们不妨用二分分治加快速幂解决这道题

【代码】：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <queue>
#include <cmath>
#include <string>
#include <cctype>
#include <map>
#include <iomanip>
                   
using namespace std;
                   
#define eps 1e-8
#define pi acos(-1.0)
#define inf 1<<30
#define linf 1LL<<60
#define pb push_back
#define lc(x) (x << 1)
#define rc(x) (x << 1 | 1)
#define lowbit(x) (x & (-x))
#define ll long long
#define c 9901

int prime[7711],total=0;

void isprime(){
    int i,j;
    prime[1]=1;
    for(int i=2; i*i<=7711; i++)
        for(int j=i*i; j<=7711; j+=i)
            prime[j]=1;
    for(int i=2; i<=7711; i++)
        if(!prime[i])
            prime[total++]=i;
}

ll mod_exp(ll a,ll b){
    ll m=1;
    while(b){
        if(b&1)
            m=(m%c)*(a%c)%c;
        b>>=1;
        a=(a%c)*(a%c)%c;
    }
    return m;
}

ll sum(ll p,ll k){
    if(k==1)
        return 1;
    ll t=(sum(p,k/2));
    if(k&1)
        return (t+mod_exp(p,k/2)*(1+p*t))%c;
    else
        return ((1+mod_exp(p,k/2))*t)%c;
}

int main(){
    ll a,b,ans;
    int i,k;
    isprime();
    while(~scanf("%lld%lld",&a,&b)){
        ans=1;
        for(int i=0; i<total && a>1; i++){
            k=0;
            while(a%prime[i]==0){
                k++;
                a/=prime[i];
            }
            ans=(ans*sum(prime[i],k*b+1))%c;
        }
        if(a>1)
            ans=(ans*sum(a,b+1))%c;
        printf("%lld\n",ans);
    }
    return 0;
}