POJ1845 Sumdiv（求所有因数和+矩阵快速幂）

题目问$A^B$的所有因数和。

根据唯一分解定理将A进行因式分解可得:A = p1^a1 * p2^a2 * p3^a3 * pn^an.
A^B=p1^(a1*B)*p2^(a2*B)*...*pn^(an*B);
A^B的所有约数之和sum=[1+p1+p1^2+...+p1^(a1*B)]*[1+p2+p2^2+...+p2^(a2*B)]*[1+pn+pn^2+...+pn^(an*B)]

知道这个，问题就变成求出A的所有质因数pi以及个数n，然后$prod(1+p_i+p_i^2+cdots+p_i^{n-1}+p_i^n)$就行了。可以构造矩阵来求：

记$S_n=p_i+p_i^2+cdots+p_i^{n-1}+p_i^n$

$$ egin{bmatrix} p_i & 1 \ 0 & 1 end{bmatrix} imes egin{bmatrix} S_n \ p_i end{bmatrix} = egin{bmatrix} S_{n+1} \ p_i end{bmatrix} $$

$$ egin{bmatrix} S_n \ p_i end{bmatrix} = egin{bmatrix} p_i & 1 \ 0 & 1 end{bmatrix} ^n imes egin{bmatrix} S_0 \ p_i end{bmatrix} $$

A忘了$pmod {9901}$，爆intWA到头疼= =

#include<cstdio>
#include<cstring>
using namespace std;
struct Mat{
    int m[2][2];
};
Mat operator*(const Mat &m1,const Mat &m2){
    Mat m={0};
    for(int i=0; i<2; ++i){
        for(int j=0; j<2; ++j){
            for(int k=0; k<2; ++k){
                m.m[i][j]+=m1.m[i][k]*m2.m[k][j];
                m.m[i][j]%=9901;
            }
        }
    }
    return m;
}
int calu(int a,int n){
    a%=9901;
    Mat e={1,0,0,1},x={a,1,0,1};
    while(n){
        if(n&1) e=e*x;
        x=x*x;
        n>>=1;
    }
    return (e.m[0][1]*a+1)%9901;
}
bool isPrime(int n){
    if(n<2) return 0;
    for(int i=2; i*i<=n; ++i){
        if(n%i==0) return 0;
    }
    return 1;
}
int main(){
    int a,b;
    scanf("%d%d",&a,&b);
    if(isPrime(a)){
        printf("%d",calu(a,b));
        return 0;
    }
    int res=1;
    for(int i=2; i*i<=a; ++i){
        if(a%i) continue;
        if(isPrime(i)){
            int cnt=0,tmp=a;
            while(tmp%i==0){
                ++cnt;
                tmp/=i;
            }
            res*=calu(i,cnt*b);
            res%=9901;
        }
        if(i!=a/i && isPrime(a/i)){
            int cnt=0,tmp=a;
            while(tmp%i==0){
                ++cnt;
                tmp/=i;
            }
            res*=calu(a/i,cnt*b);
            res%=9901;
        }
    }
    printf("%d",res);
    return 0;
}