目录php
快速幂ios
快速幂取模c++
矩阵快速幂算法
矩阵快速幂取模函数
HDU1005练习this
幂运算:\(x ^ n\)spa
根据其通常定义咱们能够简单实现其非负整数状况下的函数code
定义法:ip
int Pow (int x, int n) { int result = 1; while(n--) { result *= x; } return result; }
不难看出此时算法的时间复杂度是\(O(n)\),一旦n取较大数值,计算时间就会大大增长,极其容易出现超时的状况。ci
快速幂:
首先要在此列举两个前提:
计算机是经过二进制进行存储数据的,对二进制数据能够直接进行操做。
\(2^n+2^n=2*2^n=2^{n + 1}\)
对于\(x ^ n\),其中n能够表示为m位的二进制形式,即\(n=n_12^0+n_22^1+n_32^3+\cdots+n_m2^{m-1}\)
那么\(x ^ n=x ^ {n_12^0+n_22^1+n_32^3+\cdots+n_m2^{m-1}}\)
即\(x^n=x ^ {n_12^0}x^{n_22^1}x^{n_32^3}\cdots x^{n_m2^{m-1}}\)
根据前提1,计算机能够直接对二进制格式存储的数据进行读取,那么咱们就能够对\(n\)一个个读取再对其进行累乘。
当取到第\(a\)位时,其乘数有通项公式:\(x^{n_a2^{a-1}}\)
经过标准库math,用代码实现:
int Pow (int x, int n) { int result = 1; int a = 1; while(n) { if(n & 1) result *= round( pow(x, 1 << (a - 1)) );//round的做用在于double转int时防止丢失精度,对1进行位运算是一种计算2的n次幂的快速方法 n >>= 1; a++; } return result; }
但实际上这个调用了标准库的代码并无实现快速幂,由于仍然是采用pow()进行的运算
此处由 \(2^n+2^n=2\times2^n=2^{n + 1}\)
即\((x ^ {2 ^ {n}} )^2=x ^ {2\times 2 ^ {n}} =x ^ {2 ^ {n + 1}}\)
所以咱们能够经过对前项进行二次幂运算获得后项
先获得首项\(f(1)=x^{2^{1-1}}=x\)
即令int f = x
具体实现:
int Pow (int x, int n) { int result = 1; int f = x; while(n) { if(n & 1) result *= f; f *= f; n >>= 1; } return result; }
不难发现此时算法的时间复杂度由其次数的二进制位数而决定,即\(O(m)\)也就是\(O(log_2n)\)
另外由于此算法经常用于计算大数,因此int类型最好都换成long long类型,防止溢出。
对\(x^n\)取模运算:\(x^n\mod p\),当n极大时,一样其运算量也会增大
这就须要咱们寻求一种快速解决的方法,此时能够运用咱们上文所提到的快速幂算法
引论1:\((np+a)\mod p=a\mod p\quad (n\in\mathbb{Z} )\)
证实以下:设\(f(n,p,a)=(np+a)\mod p\quad (n\in\mathbb{Z} )\)
则由定义有\((np+a)\mod p\\=\frac{np+d}{p}-([\frac{np+d}{p}+1]-1)\\ =d-p([\frac{d}{p}+1]-1)\)
显而易见,\((np+a)\mod p=a\)与\(n\)无关
令\(n=0\)得\((np+a)\mod p=a\mod p\quad (n\in\mathbb{Z} )\\ Q.E.D\)
引论2:\((mn)\mod p=((m\mod p)(n\mod p))\mod p\)
证实以下:令\(\left\{ \begin{array}{c} m=(ip+c)& (i\in\mathbb{Z} )\\ n=(jp+d) & (j\in\mathbb{Z} )\end{array} \right.\)
则有\(\left\{ \begin{array}{c} m\mod p=c\\ n\mod p=d \end{array} \right.\)
原式\((m*n)%p\\=((ip+c)(jp+d))\mod p\\=(jip^2+jpc+ipd+cd)\mod p\\=(jip^2+jpc+ipd+cd)\mod p\\=((jip+jc+id)p+cd)\mod p\\由于(jip+jc+id)\in\mathbb{Z},由引论1得:\\=(cd)\mod p\\=((m\mod p)(n\mod p))\mod p\)
即\((mn)\mod p=((m\mod p)(n\mod p))\mod p\\ Q.E.D\)
所以对于\(x^n\mod p\)
能够写成\((x ^ {n_12^0}x^{n_22^1}x^{n_32^3}\cdots x^{n_m2^{m-1}})\mod p\\ =((x ^ {n_12^0}\mod p)(x^{n_22^1}\mod p)(x^{n_32^3}\mod p)\cdots (x^{n_m2^{m-1}}\mod p))\mod p\)
而且由以前的推定\((x ^ {2 ^ {n}} )^2=x ^ {2\times 2 ^ {n}} =x ^ {2 ^ {n + 1}}\)
有\((x ^ {2 ^ {n}} \mod p)^2\mod p =(x ^ {2 ^ {n}})^2\mod p=x ^ {2 ^ {n + 1}}\mod p\)
代码实现:
int Mod (int x, int n, int p) { int result = 1; int f = x % p; while(n) { if(n & 1) result = (result * f)%p; f = (f * f)%p; n >>= 1; } return result; }
当\(X\)为任意方块矩阵,即\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nn}\\ \end{matrix} \right|\)时
同理\(X^a\),有\(X^a=X ^ {a_12^0}X^{a_22^1}X^{a_32^3}\cdots X^{a_m2^{m-1}}\)
一样的,任意方块矩阵也适用于快速幂
代码实现:
#include <iostream> #include <cstring> using namespace std; #define R 2 struct Matrix{ int data[R][R]; }; Matrix multi(Matrix a,Matrix b,int rec){ Matrix result; memset(&(result.data), 0, sizeof(result.data)); for(int i = 0; i < rec; i++) for(int j = 0; j < rec; j++) for(int k = 0; k < rec; k++) result.data[i][j] += a.data[i][k] * b.data[k][j]; return result; } Matrix poww (Matrix x, int n) { Matrix result; memset(&(result.data), 0,sizeof( result.data)); for(int i = 0; i < R; i++) result.data[i][i] = 1; Matrix f = x; while(n) { if(n & 1) result = multi(result, f, R); f = multi(f, f, R); n >>= 1; } return result; } void MatrixPrint(Matrix target) { for(int i = 0; i < R; i++) { for(int j = 0; j < R; j++) cout << target.data[i][j]<< " "; cout << endl; } } int main() { Matrix result; memset(&(result.data), 0,sizeof( result.data)); result.data[0][0] = 1; result.data[0][1] = 2; result.data[1][0] = 3; result.data[1][1] = 4; MatrixPrint(result); cout << endl; result = poww(result, 3); MatrixPrint(result); return 0; }
输出结果:
1 2
3 437 54
81 118
运算定义:当\(X\)为任意方块矩阵,即\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nn}\\ \end{matrix} \right|\)时
则\(X\mod p=\left| \begin{matrix} x_{11} \mod p & \cdots & x_{1n}\mod p\\ \vdots & \ddots & \vdots \\ x_{n1}\mod p & \cdots & x_{nn}\mod p\\ \end{matrix} \right|\)
引论1:\((m+n)\mod p=((m\mod p)+(n\mod p))\mod p\)
证实以下:令\(\left\{ \begin{array}{c} m=(ip+c)& (i\in\mathbb{Z} )\\ n=(jp+d) & (j\in\mathbb{Z} )\end{array} \right.\)
则有\(\left\{ \begin{array}{c} m \mod p=c\\ n \mod p=d \end{array} \right.\)
原式\((m+n)\mod p\\=((ip+c)+(jp+d))\mod p\\=((i+j)p+c+d)\mod p\\由于(i+j)\in\mathbb{Z},由上文引论得:\\=(c+d)\mod p\\=((m\mod p)+(n\mod p))\mod p\)
即\((m+n)\mod p=((m\mod p)+(n\mod p))\mod p\\ Q.E.D\)
引论2:有方块矩阵\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{nn} & \cdots & x_{nn}\\ \end{matrix} \right|\),\(Y=\left| \begin{matrix} y_{11} & \cdots & y_{1m}\\ \vdots & \ddots & \vdots \\ y_{n1} & \cdots & y_{nm}\\ \end{matrix} \right|\)
则有\(XY\mod p=((X\mod p)·(Y\mod p))\mod p\)
证实以下:
令\(X=\left| \begin{matrix} x_{11} & \cdots & x_{1n}\\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nn}\\ \end{matrix} \right|\),\(Y=\left| \begin{matrix} y_{11} & \cdots & y_{1n}\\ \vdots & \ddots & \vdots \\ y_{n1} & \cdots & y_{nn}\\ \end{matrix} \right|\)
则
\(X·Y\mod p=\left| \begin{matrix} (x_{11}y_{11} + \cdots +x_{1n}y_{n1})\mod p & \cdots & (x_{11}y_{1n} + \cdots +x_{1n}y_{nn})\mod p\\ \vdots & \ddots & \vdots \\ (x_{n1}y_{11} + \cdots +x_{nn}y_{n1})\mod p & \cdots & (x_{n1}y_{1n} + \cdots +x_{nn}y_{nn})\mod p\\ \end{matrix} \right|\)
\(=\left| \begin{matrix} ((x_{11}y_{11})\mod p + \cdots +(x_{1n}y_{n1})\mod p )\mod p & \cdots & ((x_{11}y_{1n})\mod p + \cdots +(x_{1n}y_{nn})\mod p\mod p\\ \vdots & \ddots & \vdots \\ ((x_{n1}y_{11})\mod p + \cdots +(x_{nn}y_{n1})\mod p)\mod p & \cdots & ((x_{n1}y_{1n})\mod p + \cdots +(x_{nn}y_{nn})\mod p)\mod p\\ \end{matrix} \right|\)
\(=\left| \begin{matrix} (((x_{11}\mod p)(y_{11}\mod p))\mod p + \cdots +((x_{1n}\mod p)(y_{n1}\mod p))\mod p )\mod p & \cdots & (((x_{11}\mod p)(y_{1n}\mod p))\mod p + \cdots +(((x_{1n}\mod p)(y_{nn}\mod p))\mod p)\mod p\\ \vdots & \ddots & \vdots \\ (((x_{n1}\mod p)(y_{11}\mod p))\mod p + \cdots +((x_{nn}\mod p)(y_{n1}\mod p))\mod p)\mod p & \cdots & (((x_{n1}\mod p)(y_{1n}\mod p))\mod p + \cdots +((x_{nn}\mod p)(y_{nn}\mod p))\mod p)\mod p\\ \end{matrix} \right|\)
\(=\left| \begin{matrix} ((x_{11}\mod p)(y_{11}\mod p) + \cdots +(x_{1n}\mod p)(y_{n1}\mod p))\mod p & \cdots & ((x_{11}\mod p)(y_{1n}\mod p) + \cdots +((x_{1n}\mod p)(y_{nn}\mod p))\mod p\\ \vdots & \ddots & \vdots \\ ((x_{n1}\mod p)(y_{11}\mod p) + \cdots +(x_{nn}\mod p)(y_{n1}\mod p))\mod p & \cdots & ((x_{n1}\mod p)(y_{1n}\mod p) + \cdots +(x_{nn}\mod p)(y_{nn}\mod p))\mod p\\ \end{matrix} \right|\)
\(=\left( \left| \begin{matrix} x_{11}\mod p & \cdots & x_{1n}\mod p\\ \vdots & \ddots & \vdots \\ x_{n1}\mod p & \cdots & x_{nn}\mod p\\ \end{matrix} \right|\left| \begin{matrix} y_{11}\mod p & \cdots & y_{1n}\mod p\\ \vdots & \ddots & \vdots \\ y_{n1}\mod p & \cdots & y_{nn}\mod p\\ \end{matrix} \right| \right)\mod p\)
\(=((X\mod p)·(Y\mod p))\mod p \\\text{即}XY\mod p=((X\mod p)·(Y\mod p))\mod p\\Q.E.D\)
说明任意矩阵也符合快速幂取模的算法
具体实现:
#include <iostream> #include <cstring> using namespace std; #define R 2 struct Matrix{ int data[R][R]; }; Matrix multi(Matrix a,Matrix b,int rec,int p){ Matrix result; memset(&(result.data), 0,sizeof( result.data)); for(int i = 0; i < rec; i++) for(int j = 0; j < rec; j++) { for(int k = 0; k < rec; k++) result.data[i][j] += a.data[i][k] * b.data[k][j]; result.data[i][j] %= p; } return result; } Matrix poww (Matrix x, int n, int p) { Matrix result; memset(&(result.data), 0,sizeof( result.data)); for(int i = 0; i < R; i++) result.data[i][i] = 1; Matrix f = x; while(n) { if(n & 1) result = multi(result, f, R, p); f = multi(f, f, R, p); n >>= 1; } return result; } void MatrixPrint(Matrix target) { for(int i = 0; i < R; i++) { for(int j = 0; j < R; j++) cout << target.data[i][j]<< " "; cout << endl; } } int main() { Matrix result; memset(&(result.data), 0,sizeof( result.data)); result.data[0][0] = 1; result.data[0][1] = 2; result.data[1][0] = 3; result.data[1][1] = 4; MatrixPrint(result); cout << endl; result = poww(result, 3, 3); MatrixPrint(result); return 0; }
输出结果:
1 2
3 41 0
0 1
练习(HDU1005):
Problem Description
A number sequence is defined as follows:
\(f(1) = 1, f(2) = 1, f(n) = (A * f(n - 1) + B * f(n - 2)) \mod 7.\)
Given A, B, and n, you are to calculate the value of f(n).
Input
The input consists of multiple test cases. Each test case contains 3 integers A, B and n on a single line (1 <= A, B <= 1000, 1 <= n <= 100,000,000). Three zeros signal the end of input and this test case is not to be processed.
Output
For each test case, print the value of f(n) on a single line.
Sample Input
1 1 3
1 2 10
0 0 0
Sample Output
2
5
Author
CHEN, Shunbao
Source
由题意不难发现,题目输入三个数A,B,n,要求咱们求出广义的斐波那契数列(generalized Fibonacci sequence)第n项的取模
即\(f(n) = (A * f(n - 1) + B * f(n - 2)) \mod 7 = (Aa_{n-1}+Ba_{n-2})= a_n \mod 7\)
所以经过先后项关系,能够构建矩阵\(\left| \begin{matrix} a_{n} & 0\\ a_{n-1} & 0\\ \end{matrix} \right|\mod 7=\left( \left| \begin{matrix} A & B\\ 1 & 0\\ \end{matrix} \right|\left| \begin{matrix} a_{n-1} & 0\\ a_{n-2} & 0\\ \end{matrix} \right| \right)\mod 7\)
\(=\left( \left| \begin{matrix} A & B\\ 1 & 0\\ \end{matrix} \right|^{n-2}\left| \begin{matrix} a_2 & 0\\ a_1 & 0\\ \end{matrix} \right| \right)\mod 7\)
\(=\left(\left( \left| \begin{matrix} A & B\\ 1 & 0\\ \end{matrix} \right|^{n-2}\mod 7\right)\left(\left| \begin{matrix} a_2 & 0\\ a_1 & 0\\ \end{matrix} \right|\mod7\right) \right)\mod 7\)
如此,就转化为一个矩阵快速幂的问题
具体实现(谢谢提醒!原来个人solution忽略了\(n = 1\)和\(n = 2\)的状况,致使无限循环):
#include <iostream> #include <cstring> using namespace std; #define R 2 struct M{ int data[R][R]; }; M multi(M a,M b,int rec,int p){ M result; memset(&(result.data), 0,sizeof( result.data)); for(int i = 0; i < rec; i++) for(int j = 0; j < rec; j++) { for(int k = 0; k < rec; k++) result.data[i][j] += a.data[i][k] * b.data[k][j]; result.data[i][j] %= p; } return result; } M poww (M x, int n, int p) { M result; M f = x; memset(&(result.data), 0,sizeof( result.data)); for(int i = 0; i < R; i++) result.data[i][i] = 1; while(n) { if(n & 1) result = multi(result, f, R, p); f = multi(f, f, R, p); n >>= 1; } return result; } M solve(int a, int b, int n, int p) { M result; M trans; memset(&(result.data), 0,sizeof( result.data)); memset(&(trans.data), 0,sizeof( trans.data)); result.data[0][0] = 1; result.data[1][0] = 1; trans.data[0][0] = a; trans.data[0][1] = b; trans.data[1][0] = 1; trans = poww(trans, n, p); result = multi(trans, result, R, p); return result; } int main() { int a, b, n; while(true) { cin >> a >> b >> n; if(!a && !b && !n) break; if(n == 2 || n == 1 ){ cout << 1 << endl; continue; } cout << solve(a, b, n - 2, 7).data[0][0] << endl; } return 0; }