矩阵乘法优化递归式

序：

在OI比赛中，不少状况下咱们能够能经过打表（找规律）或者某些方式发现一个递归式。
例如：f(n) = f(n - 1)+f(n - 2)，（斐波那契数列）。

一般状况下，咱们计算f(n)的时间复杂度就是O(n)（分别计算f(1), f(2) ... f(n - 1)）.
可是当n很大又或者还有其余处理的复杂度一叠加便会超时。

若是不学习矩阵乘法优化的话，咱们恐怕永远不会想到计算递推式还能够进行优化。
实际上利用矩阵乘法，咱们能够将O(n)的计算递归式的复杂度降至O(logn)。

---

优化递归式的特征

形如f(n) = a1 * f(n - 1) + a2 * f(n - 2) + ... + ak * f(n - k)+c (c为常数）

---

本文讨论的范畴

• 形如 f(n) = f(n-1) + f(n-2) + .. + f(n-k) 的递归式(1)

• 形如f(n) = a1 * f(n-1) + a2 * f(n-2) + .. + ak * f(n-k) 的递归式(2)

• 形如f(n) = a1 * f(n-1) + a2 * f(n-2) + .. + ak * f(n-k) + c 的递归式(3)

实际上理解了最简单的第一个递归式的原理，就很容易理解后面的两种状况。每一个前式都是后式的一种特殊状况。

---

理论基础

首先给出斐波那契数列求第n项的O(logn)的作法，由它引出原理。

已知 :
f(n) = f(n -1) + f(n - 2); -- (1)
f(n - 1) = f(n - 1); -- (2)
由(1)(2)能够获得这样的一个式子：
[ f(n) ] = [1, 1] * [f(n - 1)]
[f(n - 1)] [1, 0] * [f(n - 2)]

这就核及到矩阵乘法的运算规则。矩阵乘法的计算方式以下所示：

A = X * Y。必须知足row(X) = colom(Y)。
A[i][j] = sum{x[i][k] * x[k][j]}。
row(A) = row(X), colom(A) = colom(Y)..

对于上式，f(n) = f(n - 1) + f(n - 2),f(n - 1) = f(n - 1) + 0，知足斐波那契数列的规则。
咱们设右边靠左的式子为A，靠右的为F。

那么计算f(n)，咱们只须要计算A^(n - 1) * F。复杂度O(n)。
可是作幂运算，能够经过快速幂将复杂度从O(n)降到O(logn)，所以总复杂度科技降到O(logn)。

经过这个例子咱们能发现什么？很显然的即是A与Y（左式与右二式，一下皆简称为A,X,Y）本质上是一个东西，所以经过迭代，直接计算出第n项。

斐波那契数列是一个最简单的例子，它也是(1)的典型例子。

广义斐波那契数列

定义f(n) = a * f(n - 1) + b * f(n - 2).
与以前的区别仅仅在于前面的系数不是1，那么构造出等式只须要照葫芦画瓢便可。

[ f(n) ] = [a, b] [f(n - 1)]
[f(n - 1)] [1, 0] [f(n - 2)],本质没有变化。
那么对于f(n) = a1 * f(n - 1) + .... + ak * f(n - k)，咱们只须要构造一个k * k的矩阵。

(*)式即是通常状况的式子，其实只须要使得第一行知足数列公式，其余行知足f(i) = f(i)便可。
所以只须要令f[i][i - 1] = 1(下标从0开始)，其余都是0便可使等式成立。

带常数与系数的递归式

经过以前的经验咱们不难看出，矩阵A,Y中的元素必定是每一次计算的时候必需的元素。此次多了一个常数c，所以c也要在A,Y中出现。
对此，咱们在须要在最后一行加上c，构形成一个(k + 1) * (k + 1)的矩阵，以下。

观察一下(**)式，第1和k+1行的最后都为1，其余新增的空位全是0，如此便构造完成了。
至于原理，再经过矩阵乘法验证一下不就行了吗？

最后给出代码实现（NOI2012 d1t3）
这道题就是一个裸的(3)式状况。

/*
About: Matrix Fast Exponentiation to calculate Recursion with Coefficient and Constant
From: NOI2012 d1t3
Description:
    Input: n, p, c, f(0), Mod, mod;
    Recursion: f(n) = (p * f(n-1) + c) % Mod;
    Output: f(n) % mod;
Auther: kongse_qi   
Date:2017/05/19
*/

#include <cstdio>
#include <cstring>

#define read(x) scanf("%d", &x)

typedef long long ll; 

long long mod;

ll Mul (ll a, ll b)
{
    ll ans = 0;
    for(; b; b >>= 1, a = (a + a) % mod) 
    {
        if(b & 1) ans = (ans + a) % mod;
    }
    return ans;
}

struct Matrix{
    ll n, m, x[5][5];
    
    Matrix(){}
    
    Matrix(int a, int b):
        n(a), m(b){}
    
    Matrix operator * (const Matrix& a) 
    {
        Matrix y(n, a.m);
        memset(y.x, 0, sizeof y.x);
        for (int i = 0; i < n; i++)
        {
            for (int k = 0; k < m; k++)
            {
                if (!x[i][k]) continue;
                for (int j = 0; j < a.m; j++)
                {
                    y.x[i][j] = (Mul(x[i][k], a.x[k][j]) + y.x[i][j]) % mod;
                }
            }
        }
        return y;
    }
    
    Matrix Poww(ll times)
    {
        Matrix y = *this, z = y;
        for(--times; times > 0; times >>= 1, z = z * z)
        {
            if(times & 1) y = y * z;
        }
        return y;
    }
};


void Create(Matrix& x, Matrix& f, ll a1, ll p, ll c)
{
    f.n = x.n = x.m = 3;
    f.m = 1;
    memset(x.x, 0, sizeof x.x);
    x.x[0][0] = p;
    x.x[0][2] = x.x[1][0] = x.x[2][2] = 1;
    f.x[0][0] = (Mul(a1, p) + c) % mod; 
    f.x[1][0] = a1;
    f.x[2][0] = c;
    return ;
}

int main()
{
    Matrix x, f;
    ll a1, n, Mod, c, p;
    scanf("%lld%lld%lld%lld%lld%lld", &mod, &p, &c, &a1, &n, &Mod);
    Create(x, f, a1, p, c);
    Matrix y = x.Poww(n - 1) * f; 
    printf("%lld\n", y.x[0][0] % Mod);
    return 0;
}

---

至于一些更简单的状况，如下有几道练习。
luogu P1962 斐波那契数列
luogu P1349 广义斐波那契数列
Vijos 1067 Warcraft III 守望者的烦恼
NOI2012 随机数生成器

最后一题值得注意的是，最后几个点乘法过程会炸出longlong，所以须要使用类快速幂的方法快速求和（积），边加边取模。

至此结束。
箜瑟_qi ** 7:58 2017.05.25**