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

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

题目链接：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 3
1 2 10
0 0 0

Sample Output

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

当输入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 <= 1018), You should solve for

g(g(g(n))) mod 109 + 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

0
1
2

Sample Output

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

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谢谢

时间： 2024-08-10 23:29:18

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

hud5451_求循环节加矩阵快速幂

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=5451 题目描述: 对于,给出x和mod,求y向下取整后取余mod的值为多少? 找循环节解析链接:http://blog.csdn.net/ACdreamers/article/details/25616461 裸题链接:http://blog.csdn.net/chenzhenyu123456/article/details/48529039 1 #include <algorithm> 2 #i

HDU1005 找规律 or 循环点 or 矩阵快速幂

http://acm.hdu.edu.cn/showproblem.php?pid=1005 1.一开始就注意到了n的数据范围 <=100 000 000,但是还是用普通的循环做的,自然TLE了,然后朴素打表,= =运行不了,(怎么可能能把数组开到那么大).再然后就想到了寻找下一个1 1 连续在一起的,那就能开始下一轮循环了. 但是,各种WA--(将数组开大一点,寻找到a[ i ] = a[ i -1 ] ==1 即跳出),这个AC代码将102改成100,150,200都可以,但是108,49

【BZOJ 1409】 Password 数论（扩展欧拉+矩阵快速幂+快速幂）

读了一下题就会很愉快的发现,这个数列是关于p的幂次的斐波那契数列,很愉快,然后就很愉快的发现可以矩阵快速幂一波,然后再一看数据范围就......然后由于上帝与集合对我的正确启示,我就发现这个东西可以用欧拉函数降一下幂,因为两个数一定互质因此不用再加一个phi(m),于是放心的乘吧宝贝!! #include <cstdlib> #include <cstring> #include <cstdio> #include <iostream> #include &

HDU 2276 &amp; FZU 1692 （循环同构优化+矩阵快速幂）

HDU 2276 题意: 给定一个01**字符串环**(2<=字符串长度<=100)然后进行m次的变换. 定义变换的规则为:如果当前位置i的左边是1(下标为0的左边是n-1),那么i就要改变状态0->1 , 1->0 比如当前的状态为100101那么一秒过后的状态为010111. 思路: 用公式表示变化状态其实就是: ai=(a(i+n?1)%n+ai)%2 换成位操作就是: ai=a(i+n?1)%n^ ai 于是我们可以建立一个变化矩阵: ??????????1100...00

HDU 1005矩阵快速幂解法 循环节解法

循环节解法: 对于公式 f[n] = A * f[n-1] + B * f[n-2]; 后者只有7 * 7 = 49 种可能,为什么这么说,因为对于f[n-1] 或者 f[n-2] 的取值只有 0,1,2,3,4,5,6 这7个数,A,B又是固定的,所以就只有49种可能值了.由该关系式得知每一项只与前两项发生关系,所以当连续的两项在前面出现过,由于公式不变,那么后面得到的一定是跟前面相重复的.所以这个时候循环节就出现了,注意循环节并不一定会是开始的 1,1:但1,1一定可以作为一个循环节,只不过

Brute-force Algorithm HDU - 3221（指数循环节 矩阵快速幂）

水题一道 推一下就是 f[n] = f[n - 1] + f[n - 2] 发现n很大 所以用矩阵快速幂就好了 还有P很大  那就指数循环节 一定要注意   这个条件 #include <iostream> #include <cstdio> #include <sstream> #include <cstring> #include <map> #include <cctype> #include <set> #incl

hdu 5895 Mathematician QSC 指数循环节+矩阵快速幂

Mathematician QSC Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 131072/131072 K (Java/Others) Problem Description QSC dream of becoming a mathematician, he believes that everything in this world has a mathematical law. Through unremitting e

矩阵快速幂基础知识

一. 先介绍以下矩阵的基础知识 矩阵:有 n 行 m 列组成一个 n*m 的矩阵 1. 矩阵的加减运算满足的条件:两个矩阵的行.列 必须相同 2. 矩阵的乘运算 满足的条件:  A矩阵的列数为 B矩阵的行数 A(ms)*B(sn)=C(mn) 得到的矩阵 C 是 m 行 n 列的 其中 c[i][j] 为A 的第 i 行与B的第j 列对应乘积的和 即:  代码: 1 const int N=100; 2 int c[N][N]; 3 void multi(int a[m][s],int b[s]

LA 3704细胞自动机——循环矩阵&amp;&amp;矩阵快速幂

题目 一个细胞自动机包含 $n$ 个格子,每个格子的取值为 $0 \sim m-1$.给定距离 $d$,则每次操作是将每个格子的值变为到它的距离不超过 $d$ 的所有格子的在操作之前的值的和除以 $m$ 的余数.给出 $n, m, d, k$ 和自动机各个格子的初始值.你的任务是计算 $k$ 次操作以后各格子的值.($1 \leq n\leq 500, 1 \leq m\leq 10^6, 0 \leq d\leq n/2, 1\leq k\leq 10^7$). 分析 如果我们把 $t$ 次操