题意：
给你N^M (N,M<=50000000) 让你求出N^M的所有约数之和mod 9901 (prime)
做法：
首先可以肯定的是 给你N 让你求出N的所有约数之和的做法
便是分解质因数并将N表示成 p1^k1*p2^k2.....pm^km 然后将这些数字分为m类
用公式(1+p1+p1^2...+p1^k1)(1+p2+p2^2...+p2^k2)...(1+pm+pm^2+...+pm^km)便可以计算得到
因为从每类中选一种pi^j乘出来 等价于将pi^j乘到这个约数中，所有的乘法可能之和便是约数和
对于N^M 其实本质一样 p1^(M*k1)*p2^(M*k2).....pm^(M*km)
问题转化为了如何求1+q+q^2+...+q^Q这个等比数列前N项和
高中数学告诉我们可以用(q^Q-1)/(q-1)这个公式快速幂解决
离散数学告诉我们可以用构造矩阵用矩阵乘法
用公式 大部分情况都是对的
但是在比如 (q^Q-1)与(q-1)都能被 9901整除的情况下求出来的肯定是0 不是正确解
（我比较愚昧 不知道如何解决 求解决方法）
用矩阵 其实也很简单 构造一个2*2的矩阵即可
我就只讲讲我的大常数sb方法
A是答案矩阵
A11 表示i次幂的时候当前这个数 A12表示i次幂的时候当前这个数加上之前的和 也就是前i项和
B是用来转移的矩阵
B11 = B12 = Num B21=0 B22=1
初始 A11=A12=1 A21=A22=0
要求 N^M次的时候只要把 B 重新构造 把A乘上B的M次 就可以了
这样就解决了此题 虽然常数不咋地。。
1 #include <cstdio>
2 #include <cstring>
3
4 #define P 9901
5 #define n 3
6 int p[P],C[n][n],Mat[n][n],tmp[n][n],N,M,ret;
7 bool mk[P];
8 inline void mkprime()
9 {
10 for (int i=2;i<P;++i)
11 if (!mk[i])
12 for (int j=i<<1;j<P;j+=i)
13 mk[j]=1;
14 for (int i=2;i<P;++i)
15 if (!mk[i]) p[++p[0]]=i;
16 }
17 inline void matmul(int A[][n],int B[][n])
18 {
19 memset(C,0,sizeof(C));
20 for (int i=1;i<3;++i)
21 for (int j=1;j<3;++j)
22 for (int k=1;k<3;++k)
23 C[i][j]=(C[i][j]+A[i][k]*(long long)B[k][j])%P;
24 memcpy(A,C,sizeof(C));
25 }
26 inline int getlog(int &N,int prime)
27 {
28 int ret=0;
29 for (;N%prime==0;N/=prime,++ret);
30 return ret;
31 }
32 inline void Mult(int prime,int log)
33 {
34 Mat[1][1]=Mat[1][2]=1;
35 tmp[1][1]=tmp[1][2]=prime;
36 tmp[2][1]=0,tmp[2][2]=1;
37 for (;log;log>>=1,matmul(tmp,tmp))
38 if (log&1) matmul(Mat,tmp);
39 ret=(ret*Mat[1][2])%P;
40 }
41 int main()
42 {
43 mkprime();
44 scanf("%d%d",&N,&M);
45 ret=1;
46 for (int i=1,j;i<=p[0]&&p[i]*p[i]<=N;++i)
47 if (N%p[i]==0) Mult(p[i],getlog(N,p[i])*M);
48 if (N>1) Mult(N,M);
49 printf("%d\n",ret);
50 return 0;
51 }
52