数论基础——循环节和矩阵快速幂的运用

首先我们来看一道基础题：

    题目链接：HDU1005 Number Sequence

    题目描述：

Number Sequence

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 147421    Accepted Submission(s): 35814

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 31 2 100 0 0

 

Sample Output

25

题意：

给定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;}

当输入2 5时，得到如下的结果：

从图示中我们可以看出，当到第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

题目链接：HDU 4291A Short problem

题目描述：

A Short problem

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 2453    Accepted Submission(s): 858

Problem Description

　　According to a research, VIM users tend to have shorter fingers, compared with Emacs users.
　　Hence they prefer problems short, too. Here is a short one:
　　Given n (1 <= n <= 10¹⁸), You should solve for 
g(g(g(n))) mod 10⁹ + 7
　　where
g(n) = 3g(n - 1) + g(n - 2)
g(1) = 1
g(0) = 0

 

Input

　　There are several test cases. For each test case there is an integer n in a single line.
　　Please process until EOF (End Of File).

 

Output

　　For each test case, please print a single line with a integer, the corresponding answer to this case.

 

Sample Input

012

 

Sample Output

0142837

题意：

给定一个嵌套表达式：

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文档：

ACM福州大学陈鸿数论

而对于这道题呢，给定的表达式是一个嵌套表达式，也就是和我们在学函数时的复合函数类型一样，那么回想一下，当我们

在求解复合函数的时候，往往会把嵌套的内层函数视为一个整体，因此我们这里也可以这样处理。但是首先我们要明白一个问题，

比如说像这样的一个简单的复合函数：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谢谢