标签:
首先我们来看一道基础题:
题目描述:
1 1 3 1 2 10 0 0 0
2 5
题意:
给定f(1),f(2)以及一个f(n)的递推式f(n) = (A * f(n - 1) + B * f(n - 2)) % 7,给定三个数,分别是A,B(1 <= A,B <= 1000),n(1 <= n <= 10^8),输出f(n).
解析:
显然,这道题如果暴力去做的话,在1000ms的时限下必然是超时的,所以我们很容易想到这道题是否可以通过找规律的方法,
看看前n项是否有什么规律,因此一开始的时候,我们可以根据样例打表找规律:
从我们打表的结果可以看出,第一个样例的结果是16个一循环的,然后我就按照打表后的规律提交了:
#include<cstdio> int main(){ int A,B,n; while(scanf("%d %d %d",&A,&B,&n)==3&&(A+B+n)!=0){ int sum,f[16]; f[0] = 1,f[1] = 1; for(int i = 2;i <= 15;++i){ f[i] = (A*f[i-1] + B*f[i-2]) % 7; } printf("%d\n",f[(n-1)%16]); } return 0; }
很快,也得到了这样的答案:
可是,这样的结果是正确的吗?显然不是,当我们写出这样的程序寻找一般规律时:
#include<cstdio> int main(){ int f1 = 1,f2 = 1,sum,A,B; scanf("%d %d",&A,&B); printf("第%d项: %d\n",1,f1); printf("第%d项: %d\n",2,f2); for(int i = 3;;++i){ sum = (A * f1 + B * f2) % 7; printf("第%d项: %d\n",i,sum); f1 = f2; f2 = sum; if(f1==1&&sum==1){ printf("第%d项开始循环\n",i); break; } } return 0; }
从图示中我们可以看出,当到第49项时,结果才开始循环,很明显,这道题是测试数据太弱或者是测试数据太单一所致。
那么,这种取巧的方法行不通,那还有什么更好的办法吗?如果你有了解过斐波那契数列的矩阵递推式的话,那么这道题很快就有
想法了。
斐波那契数列的定义(一般递推式):
矩阵表示时的递推式:
而恰好斐波那契数列的一般递推式为多项式类型,因此对于这道题,我们也可以试着往矩阵方向探索。
一般递推式:
f(n) = (A * f(n - 1) + B * f(n - 2));
经过简单的验证之后,我们可以发现
递推式的前一个矩阵为系数矩阵,设为A矩阵,递推式的后一个矩阵为B矩阵,因此有了以上的递推式,我们可以计算到
A * B^(n - 2)(注意矩阵运算不满足交换律,矩阵位置不可随意调换),最后的值再取矩阵反对角线的任意元素即可。
而要求B^(n - 2),我们利用快速幂的思想用矩阵快速幂求解,然后在求解的过程中不断取模即可。
矩阵快速幂的相关知识:
因此,有了以上思路之后,得出以下完整代码:
#include<cstdio> typedef long long ll; typedef struct node matrix; ll a,b,n; const int mod = 7; struct node{ ll _matrix[2][2]; }; //系数矩阵的初始化 void Init_orgin(matrix &temp){ temp._matrix[0][0] = 1; temp._matrix[0][1] = 1; temp._matrix[1][0] = 1; temp._matrix[1][1] = (a+b) % mod; } //指数矩阵的初始化 void Init(matrix &temp){ temp._matrix[0][0] = 0; temp._matrix[0][1] = b; temp._matrix[1][0] = 1; temp._matrix[1][1] = a; } //用于计算两个矩阵相乘 matrix multi(matrix &t1,matrix &t2){ matrix temp; temp._matrix[0][0] = (t1._matrix[0][0]*t2._matrix[0][0] + t1._matrix[0][1]*t2._matrix[1][0]) % mod; temp._matrix[0][1] = (t1._matrix[0][0]*t2._matrix[0][1] + t1._matrix[0][1]*t2._matrix[1][1]) % mod; temp._matrix[1][0] = (t1._matrix[1][0]*t2._matrix[0][0] + t1._matrix[1][1]*t2._matrix[1][0]) % mod; temp._matrix[1][1] = (t1._matrix[1][0]*t2._matrix[0][1] + t1._matrix[1][1]*t2._matrix[1][1]) % mod; return temp; } void solve(){ if(n - 2 <= 0){ printf("1\n"); } else{ int temp = n - 3; //利用一个初始的ans矩阵 matrix ans; Init(ans); matrix A = ans; while(temp > 0){ if(temp & 1){ ans = multi(ans,A); } A = multi(A,A); temp >>= 1; } Init_orgin(A); //将A矩阵修改为系数矩阵 ans = multi(A,ans); //系数矩阵与指数矩阵相乘 printf("%I64d\n",ans._matrix[0][1]); } } int main(){ while(scanf("%I64d%I64d%I64d",&a,&b,&n)==3&&(a+b+n)!=0){ solve(); } return 0; }
基础拓展:
于是,这道题才算真正的解决了,那么既然用二维矩阵快速幂的方法,能够解决f(n) = A * f(n - 1) + B * f(n - 2)
这一类多项式递推式的问题,那么我们举一反三,如果在多项式递推式后面有常数呢?即此时的多项式递推式是:
f(n) = A * f(n - 1) + B * f(n - 2) + C,那么这时我们要怎么样解决这样的一类问题呢?
同样的,既然不加常数的递推式能用矩阵解决,那么加了常数之后,当然也可以用矩阵解决,
我们可以把常数项看成多项式第三维的部分,那么经过简单推导后,我们同样可以得到这样的矩阵递推式:
同样的,递推式的前一个矩阵为系数矩阵,设为A矩阵,递推式的后一个矩阵为B矩阵,因此有了以上的递推式,
为了方便我们可以计算到A * B^(n - 1)(注意矩阵运算不满足交换律,矩阵位置不可随意调换),最后的值再取矩阵的第一行第一列元素即可。
进阶拓展:
那么,在有了以上的知识基础之后,我们可以对求解的问题进行进一步的拓展,比如有下面的一道题,这道题是来自
2012
ACM/ICPC Asia Regional Chengdu Online
题目描述:
0 1 2
0 1 42837
题意:
给定一个嵌套表达式:
g(g(g(n)))mod(10^9+7) 并且有g(n) = 3g(n - 1) + g(n - 2),g(1) = 1,g(0) = 1。n (1 <= n <= 10^18)
然后输入一个数n,要求输出嵌套表达式求值后的结果。
解析:
对于这道题,一开始的时候,我以为这道题很容易,因此三下五除二按照矩阵快速幂方式,并且计算过程中不断对mod(10^9+7)
取模,所以很快就写成了代码,并且通过了测试样例,然后信心满满的提交了,可是很快评测结果便出来了:
于是我便非常费解,为什么测试样例都通过了,可是还是wrong answer呢?
仔细看这道题给定的嵌套表达式: g(g(g(n)))mod(10^9 + 7)
从上面的分析,我们可以得知,当一个表达式对一个数取模的时候,必然会出现循环节,循环节的意思就是比如说有如下的一串数列:
1,2,3,1,2,3,1,2,3.....
这样的话,这串数列的循环节就是1,2,3了,而这串数列则是三个一循环,那么对于任何给定的数n,要求这串数列中第n个数是多少
那么只需将n对3取模即可。
证明过程我们这里就不深究了,具体可以自行百度,或者参考福州大学陈鸿的数列pdf文档:
而对于这道题呢,给定的表达式是一个嵌套表达式,也就是和我们在学函数时的复合函数类型一样,那么回想一下,当我们
在求解复合函数的时候,往往会把嵌套的内层函数视为一个整体,因此我们这里也可以这样处理。但是首先我们要明白一个问题,
比如说像这样的一个简单的复合函数:f(g(x)),我们知道,首先考虑g(x),g(x)的定义域为x的取值范围,而g(x)的值域又是f(g(x))
函数的定义域,因此g(x)和f(g(x))函数的函数图像或者说变化规律显然是不一致的。
那么有了这样的分析,我们可以再来看看这道题,这道题是三层的嵌套表达式,为了方便,我们从外到内分别称这三层为:
最外层,次外层,内层。那么由复合函数的分析方法我们可以得知:
内层函数的值域是次外层函数的定义域,次外层函数的值域是最外层函数的定义域。
那么要求最外层函数的值对一个指定数取模后的值,显然最外层函数值域的循环节即为给定的取模数,而后我们要先找到次外层
函数值域的循环节L1,而次外层函数的值域是决定了最外层函数的定义域,因此可由最外层函数值域的循环节,找到次外层函数
的循环节L1,然后按照同样的方法由向里层递推。
因此有了以上的想法,我们可以先暴力将各层函数循环节找出:
#include<cstdio> typedef long long ll; const ll mod = (ll)1e9 + 7ll; int main(){ ll l1,l2; ll sum,a0 = 0,a1 = 1; for(int i = 2;;++i){ sum = (a0 + 3 * a1) % mod; a0 = a1; a1 = sum; if(a0 == 0 && a1 == 1){ l1 = i - 1; printf("次外层函数值域的循环节为:%I64d\n",l1); //222222224一循环 break; } } a0 = 0,a1 = 1; for(int i = 2;;++i){ sum = (a0 + 3 * a1) % l1; a0 = a1; a1 = sum; if(a0 == 0 && a1 == 1){ l2 = i - 1; printf("内层函数值域的循环节为:%I64d\n",l2); //183120一循环 break; } } return 0; }
当找到各层函数的循环节之后,我们就可以按照最开始矩阵快速幂的求解方式,对这道题进行求解:
#include<cstdio> #include<algorithm> typedef struct matrix Matrix; typedef long long LL; const LL mod1 = (LL)1e9 + 7; const LL mod2 = 222222224LL; const LL mod3 = 183120LL; struct matrix{ LL a[2][2]; }; void Init_Matrix(Matrix &s){ s.a[0][0] = 0; s.a[0][1] = 1; s.a[1][0] = 1; s.a[1][1] = 3; } Matrix Multi(Matrix &s1,Matrix &s2,LL mod){ Matrix temp; temp.a[0][0] = (s1.a[0][0] * s2.a[0][0] % mod + s1.a[0][1] * s2.a[1][0] % mod) % mod; temp.a[0][1] = (s1.a[0][0] * s2.a[0][1] % mod + s1.a[0][1] * s2.a[1][1] % mod) % mod; temp.a[1][0] = (s1.a[1][0] * s2.a[0][0] % mod + s1.a[1][1] * s2.a[1][0] % mod) % mod; temp.a[1][1] = (s1.a[1][0] * s2.a[0][1] % mod + s1.a[1][1] * s2.a[1][1] % mod) % mod; return temp; } void solve(){ LL n; while(scanf("%I64d",&n)!=EOF){ if(n==0||n==1){ printf("%I64d\n",n); continue; } int T = 3; while(T--){ Matrix A,ans; Init_Matrix(ans); A = ans; if(n==0||n==1){ break; } while(n){ if(n & 1){ if(T==2){ ans = Multi(ans,A,mod3); //printf("text3 = %I64d\n",ans.a[0][0]); } else if(T==1){ ans = Multi(ans,A,mod2); //printf("text2 = %I64d\n",ans.a[0][0]); } else{ ans = Multi(ans,A,mod1); //printf("text1 = %I64d\n",ans.a[0][0]); } } if(T==2){ A = Multi(A,A,mod3); } else if(T==1){ A = Multi(A,A,mod2); } else{ A = Multi(A,A,mod1); } n >>= 1; } n = ans.a[0][0]; } printf("%I64d\n",n); } } int main(){ solve(); return 0; }
如有错误,还请指正,O(∩_∩)O谢谢
标签:
原文地址:http://blog.csdn.net/liujian20150808/article/details/51254528