约数之和-POJ1845&AcWing97（矩阵快速幂及矩阵的构建）

                                                           Sumdiv

Time Limit: 1000MS

Memory Limit: 30000K

Description

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).

Input

The only line contains the two natural numbers A and B, (0 <= A,B <= 50000000)separated by blanks.

Output

The only line of the output will contain S modulo 9901.

Sample Input

2 3

Sample Output

15

Hint

2^3 = 8. 
The natural divisors of 8 are: 1,2,4,8. Their sum is 15. 
15 modulo 9901 is 15 (that should be output). 

---

题目大意：

假设现在有两个自然数A和B，S是A^B的所有约数之和。

请你求出S mod 9901的值是多少。

我们这题先可以直接考虑暴力，之后优化一下就行了。

一个数的所有约数一定可以被它的质因子的乘积所表示，因为对于一个数A来讲,其中pi代表质数。

那么他的约数就一定为p1^0, p1^1, p1^2, ……p1^n, p1^0*p2, p1^1*p2^1……，最后它的所有约数之和就为：，我们拆开括号乘一下就知道它的正确性了。

那么对于A^B我们只需要在k1这里乘上一个B就好了：那么我们的约数之和的公式也需要从0次方加到k1*B次方。

接下来就是对每个pi进行计算了，我们可以直接写出前几项：假设f(i)代表p^i, g(i)代表

f(1)=1     g(1)=1

f(2)=p      g(2)=1+p

f(3)=p^2    g(3)=1+p+p^2

f(4)=p^3     g(4)=1+p+p^2+p^3

……

接下来就不用我说大家也知道f(i)与g(i)的结果了，也就是说，g(n)=g(n-1)*p+1

很显然这个是个递推式，我们可以使用矩阵快速幂很快得出答案。

那么也就是：, 我们只需计算一下g(0)和g(1)然后将左边的矩阵做一个n-1次方的快速幂就行了，于是对于每个质因子p，它所提供的贡献就是, 最后取g(n)就好了。初始矩阵我们直接取

以下是AC代码：

#include <bits/stdc++.h>
using namespace std;

#define IOS ios::sync_with_stdio(false)
#define ll long long

const int mod=9901;

struct Mat
{
    ll a[6][6];
};

Mat mul(Mat a,Mat b)
{
    Mat mp;
	for (int i=1; i<=2; i++)
	    for (int j=1; j<=2; j++){
	 	    mp.a[i][j]=0;
	 	    for (int k=1; k<=2; k++){
	 		    mp.a[i][j]=(mp.a[i][j]+a.a[i][k]*b.a[k][j])%mod;
		    }
	}
	return mp;
}

Mat qick(Mat a,ll b)
{
    Mat sum;
    sum.a[1][1]=1;sum.a[1][2]=0;
    sum.a[2][1]=0;sum.a[2][2]=1;
    while (b>0){
        if (b&1) sum=mul(sum,a);
        b>>=1;
        a=mul(a,a);
    }
    return sum;
}

int main()
{
    int a,b;
    IOS;
    cin>>a>>b;
    int ans=1;
    if (b==0) {
        cout<<1<<endl;return 0;
    }
    if (a==0) {
        cout<<0<<endl;return 0;
    }
    for (int i=2; i<=a; i++){
        int s=0;
        while (a%i==0){
            s++;
            a/=i;
        }
        if (!s) continue;
        Mat sum;
        sum.a[1][1]=i;sum.a[1][2]=1;
        sum.a[2][1]=0;sum.a[2][2]=1;
        sum=qick(sum,(ll)s*b-1);
        Mat mp;
        mp.a[1][1]=i+1;mp.a[1][2]=0;
        mp.a[2][1]=1;mp.a[2][2]=0;
        sum=mul(sum,mp);
        ans=ans*sum.a[1][1]%mod;
    }
    cout<<ans<<endl;
    return 0;
}