快速幂与矩阵快速幂（代码分析）

快速幂

  我们都知道3¹⁵⁶，如果普通求解，就是一个for循环，一直循环156次，如果我们把156分开，即：

     156(10)=10011100(2)   <把156十进制转化为二进制>

  这里156=2²+2³+2⁴+2⁷=b
  所以3¹⁵⁶=3^b=3^ (2²) * 3⁽²3^) * 3^ (2⁴) * 3^ (2⁷)
  那我们如何用代码表示呢？

执行代码：

ll Mul(int a,int b,int mod)
{
	ll ans=1;
	while(b)
	{
		if(b&1)//判断是否为奇数，是则进入
			ans=ans*a%mod;//-------------①
		a=a*a;//--------------②
		b>>=1;//b减半
	}
	return ans;
}

分析：

  a表示被积数，即例子中的3；b表示次数，即156；mod表示取模，可有可无，依据题目来定。
   b每次除于2，则相当于

所以每次a=aa；
如果b为奇数，则表明每一次多了一个a，即，该a不能成双结对，因此ans=ansa，a需要额外乘上。
上面的方法我们想一下，如果2^15654564足够大，那乘的过程中a是不是可能会超，因此这里提出了一种优化方法，把①和②变化一下；

优化代码：

ll _mul(int a,int b,int mod)
{
	ll ans=0;
	while(b)
	{
		if(b&1)
			ans=(ans+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return 0;
}

ll MUL(int a,int b,int mod)
{
	ll ans=1;
	while(b)
	{
		if(b&1)
			ans=_mul(ans,a,mod);
		a=_mul(a,a,mod);
		b>>=1;
	}
	return ans;
}

这里mul(a,b,mod)就相当于b*a -> a个b的数量。

---

矩阵快速幂

快速矩阵幂，首先我们要了解矩阵是如何乘积的：

这里求矩阵幂的运算和快速幂非常相像，由上图我们知道，就是把第一个的横行（a b c）与第二个的竖行（g i k） 依次相乘相加就是第一个的

因为这里需要传递的是数组，所以我们构造一个结构体来进行传递：

const int mmax=1005;
struct Mat{
	ll m[mmax][mmax];
}ans,a;

求幂的代码为：

Mat _mul(Mat a,Mat b,int mod)
{
	Mat res;
	memset(res.m,0,sizeof(res,m));
	for(int i=0;i<N;i++)//这里的N要依据之前你定义的temp数组大小,后面会提到这个N
		for(int j=0;j<N;j++)
			for(int k=0;k<N;k++)
				res.m[i][j]=(res.m[i][j]+a.m[i][k]*b.m[k][j])%mod;
	return res;

运用矩阵快速幂的重点在于如何去构造：
这里我们通过斐波那契求解来了解如何去构造
  已知这个条件为：
递推公式： F[n]=F[n−1]+F[n−2]， 由F[0]=0，F[1]=1；
所以构造为：

所以F[n]=ans.m[0][0]*F[2]+ans.m[0][1]*F[1];

代码为：

Mat _power(int p,int b)
{
	Mat temp;
	temp.m[0][0] = 1, temp.m[0][1] = 1;
	temp.m[1][0] = 1, temp.m[1][1] = 0;
	
	memset(ans.m, 0, sizeof(ans.m));
	for (int i = 0; i < 2; i++) //上面的N就是这里的2，因为temp是2*2矩阵
		ans.m[i][i] = 1;
	while (b)
	{
		if (b & 1)
			ans = _mul(ans, temp, p);
		temp = _mul(temp, temp, p);
		b >>= 1;
	}
	return ans;
}

题记：
传送门

这个题就很迷，一个简单的快速幂怎么都过不了
当时代码如下：

//#include<bits/stdc++.h>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<map>
#include<vector>

using namespace std;
#define ll long long
const ll inf = 10e9+7;

struct mat {
	ll a[103][103];
}ab;

mat _mul(mat a, mat b, int p)
{
	mat ans_1;
	memset(ans_1.a, 0, sizeof(ans_1.a));
	for (int i = 0; i < p; i++)
		for (int j = 0; j < p; j++)
			for (int k = 0; k < p; k++)
				ans_1.a[i][j] = (ans_1.a[i][j] + a.a[i][k] * b.a[k][j]%inf) % inf;
	return ans_1;
}

mat _power(int p, ll b)
{
	mat temp;
	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < p; j++)
		{
			cin >> temp.a[i][j];
		}
	}
	mat ans;
	memset(ans.a, 0, sizeof(ans.a));
	for (int i = 0; i < p; i++)
		ans.a[i][i] = 1;
	while (b)
	{
		if (b & 1)
			ans = _mul(ans, temp, p);
		temp = _mul(temp, temp, p);
		b >>= 1;
	}
	return ans;
}

int main()
{
	int p;
	ll k;
	cin >> p >> k;
	mat res = _power(p, k);
	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < p; j++)
		{
			cout << res.a[i][j] << " ";
		}
		cout << endl;
	}
	return 0;
}

结果全wa了，呵，我孝了啊啊啊啊啊

最后··········
神tm的ac了，有没有人可以帮我找找啊~~~~~~~~

这是AC代码：

//#include<bits/stdc++.h>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<map>
#include<vector>

using namespace std;
#define ll long long
#define inf 1000000007

struct mat {
	ll a[103][103];
};

mat mul(mat a, mat b, int p)
{
	mat ans_1;
	memset(ans_1.a, 0, sizeof(ans_1.a));
	for (int i = 1; i <= p; i++)
		for (int j = 1; j <= p; j++)
			for (int k = 1; k <= p; k++)
				ans_1.a[i][j] = (ans_1.a[i][j] + a.a[i][k] * b.a[k][j])%inf;
	return ans_1;
}

mat power(int p, ll b)
{
	mat temp;
	for (int i = 1; i <= p; i++)
	{
		for (int j = 1; j <= p; j++)
		{
			cin >> temp.a[i][j];
		}
	}
	mat ans = temp;
	b = b - 1;
	while (b)
	{
		if (b & 1)
			ans = mul(ans, temp, p);
		temp = mul(temp, temp, p);
		b >>= 1;
	}
	return ans;
}

int main()
{
	int p;
	ll k;
	cin >> p >> k;
	mat res = power(p, k);
	for (int i = 1; i <= p; i++)
	{
		for (int j = 1; j <= p; j++)
		{
			cout << res.a[i][j] << " ";
		}
		cout << endl;
	}
	return 0;
}