矩阵快速幂1Yet another Number Sequence R - M斐波那契数列 How Many Calls?

本文链接：https://blog.csdn.net/qq_45364953/article/details/99705663

矩阵快速幂和三个题Yet another Number Sequence，R - M斐波那契数列，How Many Calls?

首先什么是快速幂？

矩阵的幂

很经典的问题：斐波那契数列

Yet another Number Sequence

R - M斐波那契数列

How Many Calls?

首先什么是快速幂？

一种能够非常高效地计算幂运算的快速幂运算算法——反复平方法:

当我们想求 $x^{n}$ 时，将n拆分可变成 $n=2^{k1}+2^{k2}+2^{k3}\cdots$ ,就有 $x^{n}=x^{2^{k1}}x^{2^{k2}}x^{2^{k3}}\cdots$ .只要依次求 $x^{2^{ki}}$ 的同时进行计算就行了，最终得到了 $\Theta \left ( logn \right )$ 计算幂运算的算法。

那这里的k1,k2,k3,怎么算呢？我们可以把n用二进制表示，k1,k2,k3是将n表示为二进制时1那一位转换成十进制时的值。

eg： $x^{22}=x^{16}\times x^{4}\times x^{2}$ （22转成二进制数是10110）

例子：https://cn.vjudge.net/problem/UVA-10006

题意：我们把任意的1<x<n都有 $x^{n}\equiv x\left ( mod \right n )$ 成立的合数n称为Carmichael Number。对于给定的整数n，请判断它是不是Carmichael Number。

限制条件：2<n<65000

我们肯定要算x的n次方，从1到n每个数都要算它的n次方，所有此题中有n个待检查的数，如果每个数都按定义 $\Theta \left ( n \right )$ 复杂度来计算幂，则总的复杂度是 $\Theta \left ( n^{2} \right )$ ,会超时不可取，可以考虑加速幂运算的方法。

计算 $x^{n}$ 的模板：

typedef long long ll;
ll mod_pow(ll x, ll n, ll mod)  //x的n次方模mod
{
	ll res = 1;
	while (n > 0) {
		if (n & 1) //如果二进制最低位为1，则乘上x的k次
			res = res * x%mod;
		x = x * x%mod;
		n >>= 1; //位运算
	}
	return res;
}

矩阵的幂

那矩阵快速幂又是什么呢？计算一个矩阵的n次方，就是直接把上述模板中的x和res换成我们我们要计算的矩阵和结果矩阵。所以我们首先要解决一下矩阵相乘的问题才能计算res和x,x和x的相乘问题，这应该在线性代数里学的吧，可本菜鸡没学过线代，现学一下。

矩阵乘法 m行k列×k行l列得m行l列

简单的说矩阵就是二维数组，数存在里面，矩阵乘法的规则:A*B=C

其中c[i][j]为A的第i行与B的第j列对应乘积的和(比如：c11=a11*b11+a12*b21+...+a1n*bn1)，即:

矩阵Q和矩阵A相乘的代码:

matrix 是一个结构体用来定义数组大小和初始化的（也可以不初始化）

ll mod;
matrix mul(matrix Q, matrix A)
{
	matrix C;
	for (int i = 0; i <n; i++) //n是矩阵大小，也就是每一行或者每一列有多少个数
		for (int j = 0; j < n; j++)
			for (int k = 0; k < n; k++)
				C.mat[i][j] = (C.mat[i][j] + (A.mat[i][k] * Q.mat[k][j]) % mod) % mod;
	return C;
}

那么计算矩阵A的n次方的模板就是把上面两个模板加起来：

ll mod,n;
matrix mul(matrix Q, matrix A)
{
	matrix C;
	for (int i = 0; i <m; i++)
		for (int j = 0; j < n; j++)
			for (int k = 0; k < n; k++)
				C.mat[i][j] = (C.mat[i][j] + (A.mat[i][k] * Q.mat[k][j]) % mod) % mod;
	return C;
}
matrix powmul(matrix A, int m)
{
	matrix ans;
        for(int i = 0; i < n; i++)
            ans.mat[i][i] = 1;  //单位矩阵，相当于初始化为1，也就是说任何矩阵乘了它还是它自己：1×A=A
	while (m > 0)
	{
		if (m & 1)
			ans = mul(ans, A);
		A = mul(A, A);
		m >>= 1;
	}
	return ans;
}

很经典的问题：斐波那契数列

问题1：

Yet another Number Sequence

UVA - 10689

Let’s deﬁne another number sequence, given by the following function:
f(0) = a

f(1) = b

f(n) = f(n−1) + f(n−2), n > 1

When a = 0 and b = 1, this sequence gives the Fibonacci Sequence. Changing the values of a and b, you can get many diﬀerent sequences. Given the values of a, b, you have to ﬁnd the last m digits of f(n).
Input
The ﬁrst line gives the number of test cases, which is less than 10001. Each test case consists of a single line containing the integers a b n m. The values of a and b range in [0,100], value of n ranges in [0,1000000000] and value of m ranges in [1,4].
Output
For each test case, print the last m digits of f(n). However, you should NOT print any leading zero.
Sample Input
4

0 1 11 3

0 1 42 4

0 1 22 4

0 1 21 4
Sample Output
89

4296

7711

946

题意：给出F1，F2，n，m，计算结果Fn并且只输出最后m个数

首先要把递推式表达成矩阵

然后怎么求Fn呢？我们要从F0开始推，首先令n=1的时候得F1和F0的值，它们乘左边的转移矩阵1 1 1 0，变大变成了F2，F1。然后继续乘这个转移矩阵1 1 1 0又变大了变成F3，F2，直至变大到了Fn就不继续乘这个转移矩阵了，假设这个转移矩阵乘了n次

表达式就是： F1,F2已经知道了，我们要做的就是用矩阵快速幂求转移矩阵的n次方。

单位矩阵可以直接把它换掉乘了它相当于乘了1，换成，反正最后算出来转移矩阵的n次方还要乘这个矩阵，一开始就乘了它也省事。

代码：

#include<iostream>
#include<cstring>
#include<vector>
#include<cmath>
using namespace std;
typedef long long ll;
struct matrix
{
    ll mat[2][2];
	matrix() {
		memset(mat, 0, sizeof(mat)); //在这初始化全置为零在多次输入时可能有遗留前面数据导致错误
	}
}U,F;
ll mod;
//两个矩阵相乘：
matrix mul(matrix Q, matrix A)
{
	matrix C;
	for (int i = 0; i <2; i++)
		for (int j = 0; j < 2; j++)
			for (int k = 0; k < 2; k++)
				C.mat[i][j] = (C.mat[i][j] + (A.mat[i][k] * Q.mat[k][j]) % mod) % mod;
	return C;
}
//矩阵A的m次方：
matrix powmul(matrix A, int m)
{
	matrix ans = F; //直接等于a,b
	while (m > 0)
	{
		if (m & 1)
			ans = mul(ans, A);
		A = mul(A, A);
		m >>= 1; //=号不要写掉了啊
	}
	return ans;
}
int main() {
	int N; cin >> N;
	while (N--)
	{
		ll n,a,b,m;
		cin >> a >> b >> n >> m;
                //转移矩阵1 1 1 0
		U.mat[0][0] = U.mat[0][1] = U.mat[1][0] = 1;
		U.mat[1][1] = 0;
                //初始的F1,F2的值
		F.mat[0][0] = b;
		F.mat[0][1] = 0;
		F.mat[1][0] = a;
		F.mat[1][1] = 0;
		mod = 1;
		for (int i = 1; i <= m; i++) 
			mod = mod * 10;
		matrix ans = powmul(U, n);
		cout << ans.mat[1][0] << endl;
	}
	return 0;
}

总结：

小细节:最后输出m位数用对10的m次方取模来实现

for (int i = 1; i <= m; i++)
mod = mod * 10;

如果算的是矩阵的n次方的话，Fn存在ans[1][0]里，如果算的是矩阵的n-1次方的话Fn存在ans[0][0]里，可以自己动手算一下。

问题2：

R - M斐波那契数列

HDU - 4549

M斐波那契数列F[n]是一种整数数列，它的定义如下：

F[0] = a
F[1] = b
F[n] = F[n-1] * F[n-2] ( n > 1 )

现在给出a, b, n，你能求出F[n]的值吗？

Input

输入包含多组测试数据；
每组数据占一行，包含3个整数a, b, n（ 0 <= a, b, n <= 10^9 ）

Output

对每组测试数据请输出一个整数F[n]，由于F[n]可能很大，你只需输出F[n]对1000000007取模后的值即可，每组数据输出一行。

Sample Input

0 1 0
6 10 2

Sample Output

0
60

分析：写出几组Fn的值，观察a,b的指数有什么规律。

$F_{0}=ab^{0},F_{1}=a^{0}b,F_{2}=ab,F_{3}=ab^{2},F_{4}=a^{2}b^{3},F_{5}=a^{3}b^{5}$

发现：a的指数变化：1 0 1 1 2 3···

b的指数变化：0 1 1 2 3 5···

所以a,b的指数时分别以F0=1,F1=0和F0=0,F1=1为前两项的斐波那契数列，但是a,b的指数在n很大时会变得非常大，这里我们需要用到

费马小定理：若p是素数， $gcd\left ( a,p \right )=1$ ， $a^{p-1}\equiv 1\left ( mod \right p )$ ，实际上，它是欧拉定理的一个特殊情况。

在此题中：若 $a^{n}mod\: p$ 中n很大，则可以简化为 $a^{n}mod\; p=a^{[n\, mod\, (p-1)]}\, mod\, p$

证明如下：

n=t*(p-1)+r，其中r为n除以(p-1)的余数，即为n mod (p-1)。

a^n=(a^(p-1))^t * a^r 1^t * a^r a^r (mod p)

费马小定理的推广：如果p为质数， $x^{p}-x$ （x是任意正整数）必能被p整除

所以此题只需根据斐波那契数列用矩阵快速幂求出a，b的指数模1000000006的值令它为m和k，再快速幂求a^m和b^k，它俩相乘再模上1000000007.

代码：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;

struct matrix
{
    ll mat[15][15];
	matrix() {
		memset(mat, 0, sizeof(mat));
	}
};

ll mod = 1000000006;
matrix mul(matrix Q, matrix A)
{
	matrix C;
	for (int i = 1; i <= 2; i++)
		for (int j = 1; j <= 2; j++)
			for (int k = 1; k <= 2; k++)
				C.mat[i][j] = (C.mat[i][j] + (A.mat[i][k] * Q.mat[k][j]) % mod) % mod;
	return C;
}
matrix powmul(matrix A, int m)
{
	matrix ans;
	for (int i = 1; i <= 2; i++)
		ans.mat[i][i] = 1;
	while (m > 0)
	{
		if (m & 1)
			ans = mul(ans, A);
		A = mul(A, A);
		m >>= 1;
	}
	return ans;
}

ll modx = 1000000007;
ll pow(ll a, ll b)
{
	ll res = 1;
	while (b)
	{
		if (b & 1)
			res = (res*a) % modx;
		a = (a*a) % modx;
		b >>= 1;
	}
	return res;
}
int main() {
	long long a, b, n;
	while (cin >> a >> b >> n)
	{
		if (n == 0)
			cout << a % 1000000007 << endl;
		else if (n == 1)
			cout << b % 1000000007 << endl;
		else
		{
			matrix A;
			A.mat[1][1] = A.mat[1][2] = A.mat[2][1] = 1;
			A.mat[2][2] = 0;
			A = powmul(A, n - 1);
			long long m, k;
			m = A.mat[1][1]; k = A.mat[1][2];
			cout << (pow(a, k)*pow(b, m)) % 1000000007 << endl;
		}
	}
	return 0;
}

总结：

知识点：

费马小定理：若p是素数， $gcd\left ( a,p \right )=1$ ， $a^{p-1}\equiv 1\left ( mod \right p )$ ，实际上，它是欧拉定理的一个特殊情况。

费马小定理的推广：如果p为质数， $x^{p}-x$ （x是任意正整数）必能被p整除

个人易错：n=0和n=1要单独讨论，得出的结果也要模上mod

问题三：

How Many Calls?

UVA - 10518

The ﬁbonacci number is deﬁned by the following recurrence:

• fib(0) = 0

• fib(1) = 1

• fib(n) = fib(n−1) + fib(n−2)

But we’re not interested in the ﬁbonacci numbers here. We would like to know how many calls does it take to evaluate the n-th ﬁbonacci number if we follow the given recurrence. Since the numbers are going to be quite large, we’d like to make the job a bit easy for you. We’d only need the last digit of the number of calls, when this number is represented in base b.
Input

Input consists of several test cases. For each test you’d be given two integers n (0 ≤ n < $2^{63}$ −1), b (0 < b ≤ 10000). Input is terminated by a test case where n = 0 and b = 0, you must not process this test case.
Output
For each test case, print the test case number ﬁrst. Then print n, b and the last digit (in base b) of the number of calls. There would be a single space in between the two numbers of a line.
Note that the last digit has to be represented in decimal number system.
Sample Input
0 100

1 100

2 100

3 100

10 10

0 0
Sample Output
Case 1: 0 100 1

Case 2: 1 100 1

Case 3: 2 100 3

Case 4: 3 100 5

Case 5: 10 10 7

题意：这道题是让我们求在我们求Fn的时候调用了多少次（假设这个次数是ans)上面那个递推循环，所以现在我们要找找ans有什么规律。

第几个数：n=0 n=1 n=2 n=3 n=4 n=5 n=6 n=7 n=8 ···

ans(n)的值: 1 1 3 5 9 15 25 41 67 ···

F(n)的值： 1 1 2 3 5 8 13 21 34 ···（这里的F0=1，F1=1，和题中的斐波那契数列不一样更方便计算）

观察发现规律：1.ans(n)=2×F(n)-1，2.ans(n)=ans(n-1)+ans(n-2)+1

我用的第一个规律，觉得更简单可以把以前的模板直接复制过来。第二个也可以用不过要重新写矩阵，觉得麻烦有时间再补上代码。

所以只需要矩阵快速幂根据斐波那契数列计算出Fn的值再利用第一个规律算出ans(n)就是答案了。题目中说“the last digit (in base b) of the number of calls”实际是要让结果模上b。

代码：

#include<iostream>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
struct matrix
{
	ll mat[2][2];
	matrix() {
		memset(mat, 0, sizeof(mat));
	}
}U, F;
ll mod;
matrix mul(matrix Q, matrix A)
{
	matrix C;
	for (int i = 0; i < 2; i++)
		for (int j = 0; j < 2; j++)
			for (int k = 0; k < 2; k++)
				C.mat[i][j] = (C.mat[i][j] +(A.mat[i][k] * Q.mat[k][j])%mod) % mod;
	return C;
}
matrix powmul(matrix A, ll m)
{
	matrix ans = F;
	while (m > 0)
	{
		if (m & 1)
			ans = mul(ans, A);
		A = mul(A, A);
		m >>= 1;
	}
	return ans;
}
int main() {
	ll n;
	int cut = 0;
	while(cin>>n>>mod){
		if (n == 0 && mod == 0)
			break;
		cut++;
		U.mat[0][0] = 1;U.mat[0][1] = 1;
		U.mat[1][0] = 1;U.mat[1][1] = 0;
		
		F.mat[0][0] = 1; F.mat[1][0] = 1;
		F.mat[1][1] = 0; F.mat[0][1] = 0;
		if (n == 0 || n == 1)
			cout << "Case " << cut << ": " << n << " " << mod << " " << "1"<< endl;
		else
		{
			matrix ans = powmul(U, n);
			cout << "Case " << cut << ": " << n << " " << mod << " " << (ans.mat[1][0]* 2- 1+mod) % mod << endl; 
                        //加上mod很重要，如果ans.mat[1][0]是0的话减一会变成负数，错了好几次都没发现
		}
	}
	return 0;
}