一种写程序快速计算常系数线性齐次递推关系的指定项的方法

上次学习到了利用Fibonacci恒等式在log(n)时间内计算出F(n)的方法。这种方法也适用于任何常系数线性齐次递推数列的计算。在这里把能想到的东西小结一下。
例一：扩展Euclid算法。
两个整数a,b的最高公因数为gcd(a,b)。扩展Euclid算法在计算出gcd(a,b)的同时计算出gcd(a,b)=sa+tb中的系数s,t。因为：

这样，只要计算出中间那些矩阵的乘积，就很容易计算出系数s,t。其妙处就在于将(a,b)和(b,r1)通过矩阵联系起来。
例二：Fibonacci数。
F(n)=F(n-1)+F(n-2)，这是一个二阶递推式。所以存在一个二阶矩阵。类似上面的方法，将F(n),F(n-1)和F(n-1),F(n-2)通过矩阵联系起来：

计算n次方可以使用快速模幂算法在log(n)内完成。而二阶矩阵的一次乘法所用时间是常量。所以可以在log(n)时间复杂度内计算F(n)。
还可以把上述等号两边的Fibonacci数合成一个2*2矩阵。即：

通过最后那个式子可以看出F(n+1)与F(n)互素。
以上两个算法的共同点是通过矩阵建立起一个递推关系。常系数线性递推关系的解总是指数形式的。都可以使用这样的方法在对数时间复杂度里计算出指定项的值。
下面的代码是用来计算指定项的Fibonacci数和另一个递推关系A(k)=A(k-1)+A(k-2)+A(k-3),A(0)=1,A(1)=A(2)=2的。代码用C++写，比较冗长。先记在这里。

// 利用矩阵恒等式快速计算第N个fibonacci数 #include <iostream> #include <cstring> #include <cstddef> #include <cassert> using namespace std; template <typename T, size_t D> class Matrix { T data[D][D]; static size_t cnt; public: static Matrix<T,D> Identity; public: Matrix() { memset(data, 0, sizeof(data)); if(cnt == 0) { memset(Identity.data, 0, sizeof(Identity.data)); for(int i = 0; i < D; i++) Identity.data[i][i] = 1; } cnt++; } Matrix(const T* arr, size_t n) { memset(data, 0, sizeof(data)); size_t k = 0, r = 0, c = 0; while(k < n && r < D) { data[r][c++] = arr[k++]; if(c == D) { r++; c = 0; } } if(cnt == 0) { memset(Identity.data, 0, sizeof(Identity.data)); for(int i = 0; i < D; i++) Identity.data[i][i] = 1; } cnt++; } Matrix(const Matrix& m) { memcpy(data, m.data, sizeof(data)); } ~Matrix() { cnt--; } Matrix& operator=(const Matrix& m) { if(this != &m) { memcpy(data, m.data, sizeof(data)); } return *this; } Matrix operator*(const Matrix& m) { Matrix result; for(size_t i = 0; i < D; i++) { for(size_t j = 0; j < D; j++) { result.data[i][j] = 0; for(size_t k = 0; k < D; k++) { result.data[i][j] += data[i][k] * m.data[k][j]; } } } return result; } Matrix& operator*=(const Matrix& m) { return (*this = *this * m); } T* const operator[](int index) { assert(index >= 0 && index < D); return data[index]; } Matrix power(size_t n) { Matrix result = Identity, tmp = *this; //while(n) { // if(n & 1) { // result *= tmp; // n--; // } else { // tmp *= tmp; // n /= 2; // } //} for(; n; n >>= 1) { if(n & 1) result *= tmp; tmp *= tmp; } return result; } friend istream& operator>>(istream& in, Matrix<T,D>& m) { for(size_t i = 0; i < D; i++) for(size_t j = 0; j < D; j++) in >> m.data[i][j]; return in; } friend ostream& operator<<(ostream& out, const Matrix<T,D>& m) { out.setf(ios::left, ios::adjustfield); for(size_t i = 0; i < D; i++) { for(size_t j = 0; j < D; j++) { out.width(4); out << m.data[i][j]; } out << endl; } return out; } }; template <typename T, size_t D> Matrix<T,D> Matrix<T,D>::Identity; template <typename T, size_t D> size_t Matrix<T,D>::cnt = 0; int main() { Matrix<int, 2> fib; int n; fib[0][0] = 1; fib[0][1] = 1; fib[1][0] = 1; fib[1][1] = 0; cin >> n; int F[30] = {0, 1}; for(int i = 2; i < 30; i++) F[i] = F[i-1] + F[i-2]; cout << "F[" << n << "] is: " << F[n] << endl; cout << "F[" << n << "] = " << fib.power(n-1)[0][0] << endl; // A(k) = A(k-1) + A(k-2) + A(k-3), A(0) = 1, A(1) = 2, A(2) = 2 int AM[] = {1,1,1,1,0,0,0,1,0}; Matrix<int, 3> am(AM, sizeof(AM)); cin >> n; Matrix<int, 3> ret = am.power(n-2); cout << "A[" << n << "] = " << ret[0][0] * 2 + ret[0][1] * 2 + ret[0][2] << endl; int A[30] = {1, 2, 2}; for(int i = 3; i < 30; i++) A[i] = A[i-1] + A[i-2] + A[i-3]; cout << "A[" << n << "] is : " << A[n] << endl; return 0; }