矩阵乘法优化递归式

序：

在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

Central: 
    construct a Matrix X[a, 0, 1] and Matrix F[f(1)]. Cal A = X ^ (n-1) * F. The ans is A[0][0] % mod. 
                        [1, 0, 0]             [f(0)]  
                        [0, 0, 1]             [  c ]
                                               [ f(n) ]   [a, 0, 1]   [f(n-1)]
    Recursion f(n) = (p * f(n-1) + c) % Mod == [f(n-1)] = [1, 0, 0] * [f(n-2)] (all % mod); 
                                               [   c  ]   [0, 0, 1]   [   c  ]
*/

#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