poj 1845 Sumdiv 矩阵法求幂的和

题意：

给A,B,求A^B所有因子模9901的和(0 <= A,B <= 50000000)。

分析：

因素分解后关键是求1+p^1+p^2+....p^k的和，因为等比数列和公式涉及除法，不用乘法逆元的话不能直接用公式。设s(k+1)=p^0+p^1+...p^k,s0=0,则

然后就能2分幂快速求了。

代码：

//poj 1845
//sep9
#include <iostream>
using namespace std;
typedef __int64 INT;
const INT mod=9901;
const int list_len=32;
INT A,B;
INT num[list_len],v[list_len];
INT mat[2][2],tmp[2][2],a[2][2];

void mul(INT z[][2],INT x[][2],INT y[][2])
{
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			tmp[i][j]=0;
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			for(int k=0;k<2;++k)
				tmp[i][j]=(tmp[i][j]+x[i][k]*y[k][j])%mod;
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			z[i][j]=tmp[i][j];				
}

INT f(INT p,INT tmp_k){
	INT k=tmp_k*B+1;
	mat[0][0]=p;mat[0][1]=0;
	mat[1][0]=1;mat[1][1]=1;
	a[0][0]=1;a[0][1]=0;
	a[1][0]=0;a[1][1]=1;
	while(k){
		if(k&1)
			mul(a,a,mat);//a=a*mat;
		mul(mat,mat,mat);	//mat=mat*mat;
		k/=2;
	}
	return a[1][0];
}

int main()
{
	scanf("%I64d%I64d",&A,&B);
	int l=0;
	memset(num,0,sizeof(num));
	for(INT i=2;i*i<=A;++i)
		if(A%i==0){
			v[l]=i;
			while(A%i==0){
				++num[l];
				A/=i;
			}
			++l;
		}
	if(A!=1){
		v[l]=A;
		num[l++]=1;
	}
	INT ans=1;
	for(int i=0;i<l;++i)
		ans=(ans*f(v[i],num[i]))%mod;
	printf("%I64d",ans);
	return 0;	
}