POJ1845——等比数列的分治思想

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的因子和
方法
将a进行质因子分解，变为p1,p2,p3…pk,显然题目的答案就是(p1^0+p1^1+p1^2+…+p1^a1*b)(p2^0+p2^1+p2^2+p2^3+…+p2^a2*b)…(pk^0+pk^1+pk^2+pk^3+…+pk^ak*b)。

所以关键就是在于求等比数列pi。运用分治的思想，pi^0+pi^1+pi^2+pi^3+…+pi^ai*b 可以从中部分成两段看成

pi^0+pi^1+pi^2+pi^3+…+pi^(ai*b/2)
pi^(ai*b/2+1)+pi^(ai*b/2+2)+pi^(ai*b/2+3)+…+pi^(ai*b/2)

合并后变成(1+pi^(ai*b/2+1))*？。
这里的问号要分奇偶讨论，详见代码。

#include<cstdio>
#include<cmath>
using namespace std;
const int tt=9901;
int a,b,n,ans;
int qsm(int a,int b){
    int w=a%tt,sum=1;
    while(b>0){
        if ((b&1)==1) sum=(sum*w)%tt;
        w=(w*w)%tt;
        b=b>>1;
    }
    return sum;
}
int work(int x,int y){
    if (y==0) return 1;
    if ((y&1)==1) return ((long long) (1+qsm(x,y/2+1))*work(x,y>>1))%tt;
    else return ((long long) (1+qsm(x,y/2+1))*work(x,(y>>1)-1)+qsm(x,y>>1))%tt;
}
int main(){
    freopen("exam.in","r",stdin);
    freopen("exam.out","w",stdout);
    scanf("%d%d",&a,&b);
    n=a;ans=1;
    for (int i=2;i<=sqrt(n);i++){
        int sum=0;
        while (a%i==0) a/=i,sum++;
        ans=(ans*work(i,sum*b))%tt;
    }
    if (a!=1) ans=(ans*work(a,b))%tt;
    printf("%d\n",ans);
    return 0;
}