Project Euler 258

Project Euler 258

题目

A lagged Fibonacci sequence

A sequence is defined as:

  • \(g_k = 1\), for \(0 \leq k \leq 1999\)
  • \(g_k = g_{k-2000} + g_{k-1999}\), for \(k \geq 2000\)

Find \(g_k \bmod 20092010\) for \(k = 10^{18}\).

哈密顿-凯莱(Caylay-Camilton)定理

如果一个\(n\)阶矩阵\(A\)的特征多项式为\(f(\lambda)=|A-\lambda I|=\sum_{i=0}^n b_i\lambda^i\),其中\(b^i\)为系数。那么矩阵多项式\(f(A)=\sum_{i=0}^nb_iA^i\)满足\(f(A)=O\).

线性递推

对于一个线性递推\(a_n=\sum_{j=1}^kc_ja_{n-j}\),可以构造如下矩阵进行\(O(k^3\log n)\)的线性递推计算:

\[ \begin{bmatrix} c_1 & c_2 & c_3 &\cdots &c_{k-1} & c_k\\ 1 & 0 & 0 & \cdots & 0 & 0\\ 0 & 1 & 0 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 &\cdots & 0 & 0\\ 0 & 0 & 0 &\cdots & 1 & 0 \end{bmatrix} \begin{bmatrix} a_{n-1}\\ a_{n-2}\\ a_{n-3}\\ \vdots\\ a_{n-k+1}\\ a_{n-k} \end{bmatrix} = \begin{bmatrix} a_n\\ a_{n-1}\\ a_{n-2}\\ \vdots\\ a_{n-k+2}\\ a_{n-k+1} \end{bmatrix} \]

\(A=\begin{bmatrix} c_1 & c_2 & c_3 &\cdots &c_{k-1} & c_k\\ 1 & 0 & 0 & \cdots & 0 & 0\\ 0 & 1 & 0 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 &\cdots & 0 & 0\\ 0 & 0 & 0 &\cdots & 1 & 0 \end{bmatrix},x= \begin{bmatrix} a_{k-1}\\ a_{k-2}\\ a_{k-3}\\ \vdots\\ a_1\\ a_0 \end{bmatrix}\)

解决方案

可以知道,矩阵\(A\)的特征多项式为\(f(\lambda)=|A-\lambda I|=\lambda^k-(\sum_{j=1}^k c_j \lambda^{k-j})\)

代入Caylay-Camilton定理,可以发现,\(A^k=\sum_{j=1}^k c_j A^{k-j}\)

由题目可知\(A^{2000}=A+I\)

Caylay-Camilton定理说明,任意\(m\geq k\)的矩阵幂\(A^k\)都可以表示成一个为\(\sum_{i=0}^{k-1} v_iA^i\)的矩阵多项式,其中\(v_i\)为系数。

而为了求\(a_n\),可以计算\(A^nx\),然后取向量\(A^nx\)的最后一个值\((A^nx)_{k-1}\)即可。

那么令\(A^n=\sum_{i=0}^{k-1} v_iA^i\)

可以知道,\(A^nx=\sum_{i=0}^{k-1} v_iA^ix\)

系数\(v_i\)可以通过快速幂算法进行多项式卷积进行计算,然后将卷积结果系数大于等于\(k\)的部分通过矩阵多项式回退到小于\(k\)

最终结果为\(a_n=\sum_{i=0}^{k-1} v_i\cdot a_i\)

原因:容易发现,计算\(Ax\)只是将整个向量\(x\)向下平移了一次,新的值将在向量\(x\)上面填充。

因此,不必要将整个向量\(A^ix\)计算出来,只需要将最后一个值\((A^ix)_{k-1}\)计算出来,并与系数相乘即可,因为此时的\(i\)不会发生\(i\geq k\)的情况。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int M=2000;
ll N=1e18;
ll mod = 20092010;
// a代表递推式的前M个值
ll a[M];
// p代表Hamilton-Cayley而来的矩阵多项式的系数,题中为F(A)=A^2000=A+I
ll p[M];
// b代表A^(2^i)时所有的A^(2^i)=sum from j=0 to M-1 A^j*b_j时系数b_j的值。
ll b[M];
// v代表解决方案所指的v。
ll v[M];
void mul_ploy(ll a[M],ll b[M],ll p[M]){
//做的是乘法,对应的就是多项式系数的卷积。
ll c[M<<1];memset(c,0,sizeof(c));
for(int i=0;i<M;i++)
for(int j=0;j<M;j++)
c[i+j]=(c[i+j]+a[i]*b[j])%mod;
//利用递推式F(A)的值,将>=M的所有临时系数转化成<M的系数。
for(int i=M*2-1;i>=M;i--)
for(int j=0;j<M;j++)
c[i-M+j]=(c[i-M+j]+c[i]*p[j])%mod;
for(int i=0;i<M;i++)
a[i]=c[i];
}
int main(){
for(int i=0;i<M;i++)
a[i]=1;
p[0]=p[1]=1;
b[1]=1;
v[0]=1;
for(;N;N>>=1){
if(N&1) mul_ploy(v,b,p);
mul_ploy(b,b,p);
}
ll ans=0;
//需要注意的是,当向量和一个矩阵相乘后,我们只要向量中的第一个值。因为向量乘一次这个矩阵递推式,那么所有值都会前移一次。
for(int i=0;i<M;i++)
ans=(ans+v[i]*a[i])%mod;
printf("%lld\n",ans);
}

如果觉得我的文章对您有用,请随意打赏。您的支持将鼓励我继续创作!
Ujimatsu Chiya 微信 微信
Ujimatsu Chiya 支付宝 支付宝