題意：
給你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