快速幂算法(SOJ 4454)

SOJ 4454: 写作业

问题：给出递推式$f_{1}=2$, $f_{2}=2$, $f_{n}=f_{n-1}+3f_{n-2}+1$, 求解$f_{n}\%1000000007$.

分析：类似于求解斐波那契数(参考这里)，利用矩阵的快速幂算法，关键是写出矩阵形式。下面是一种最自然的方式，

$\left[\begin{array}{c}f_{n}\\f_{n-1}\end{array}\right]=\left[\begin{array}{cc}1 &3\\ 1 &0 \end{array}\right]\left[\begin{array}{c}f_{n-1}\\f_{n-2}\end{array}\right]+\left[\begin{array}{c}1\\0\end{array}\right]$

因为最后有个常数矩阵，我们无法使用快速幂算法，修改形式如下：

$\left[\begin{array}{c}f_{n}\\f_{n-1}\\1\end{array}\right]=\left[\begin{array}{ccc}1 &3 &1\\ 1 &0 &0\\ 0 &0 &1\end{array}\right]\left[\begin{array}{c}f_{n-1}\\f_{n-2}\\1\end{array}\right]$.

代码：

#include<iostream>
#include<cmath>
using namespace std;
long long I[3][3];
long long con = 1000000007;
void matMul(long long flag)
{
    long long a, b, c, d, e, f, g, h, i;
    a = (((I[0][0] % con)*(I[0][0] % con)) % con + ((I[0][1] % con)*(I[1][0] % con)) % con + ((I[0][2] % con)*(I[2][0] % con)) % con) % con;
    b = (((I[0][0] % con)*(I[0][1] % con)) % con + ((I[0][1] % con)*(I[1][1] % con)) % con + ((I[0][2] % con)*(I[2][1] % con)) % con) % con;
    c = (((I[0][0] % con)*(I[0][2] % con)) % con + ((I[0][1] % con)*(I[1][2] % con)) % con + ((I[0][2] % con)*(I[2][2] % con)) % con) % con;
    d = (((I[1][0] % con)*(I[0][0] % con)) % con + ((I[1][1] % con)*(I[1][0] % con)) % con + ((I[1][2] % con)*(I[2][0] % con)) % con) % con;
    e = (((I[1][0] % con)*(I[0][1] % con)) % con + ((I[1][1] % con)*(I[1][1] % con)) % con + ((I[1][2] % con)*(I[2][1] % con)) % con) % con;
    f = (((I[1][0] % con)*(I[0][2] % con)) % con + ((I[1][1] % con)*(I[1][2] % con)) % con + ((I[1][2] % con)*(I[2][2] % con)) % con) % con;
    g = (((I[2][0] % con)*(I[0][0] % con)) % con + ((I[2][1] % con)*(I[1][0] % con)) % con + ((I[2][2] % con)*(I[2][0] % con)) % con) % con;
    h = (((I[2][0] % con)*(I[0][1] % con)) % con + ((I[2][1] % con)*(I[1][1] % con)) % con + ((I[2][2] % con)*(I[2][1] % con)) % con) % con;
    i = (((I[2][0] % con)*(I[0][2] % con)) % con + ((I[2][1] % con)*(I[1][2] % con)) % con + ((I[2][2] % con)*(I[2][2] % con)) % con) % con;
    I[0][0] = a;
    I[0][1] = b;
    I[0][2] = c;
    I[1][0] = d;
    I[1][1] = e;
    I[1][2] = f;
    I[2][0] = g;
    I[2][1] = h;
    I[2][2] = i;
    if (flag)
    {
        a = (I[0][0] % con + I[0][1] % con) % con;
        b = ((I[0][0] % con)*3) % con;
        c = (I[0][0] % con + I[0][2] % con) % con;
        d = (I[1][0] % con + I[1][1] % con) % con;
        e = ((I[1][0] % con) * 3) % con;
        f = (I[1][0] % con + I[1][2] % con) % con;
        g = (I[2][0] % con + I[2][1] % con) % con;
        h = ((I[2][0] % con) * 3) % con;
        i = (I[2][0] % con + I[2][2] % con) % con;
        I[0][0] = a;
        I[0][1] = b;
        I[0][2] = c;
        I[1][0] = d;
        I[1][1] = e;
        I[1][2] = f;
        I[2][0] = g;
        I[2][1] = h;
        I[2][2] = i;
    }
}
int main()
{
    long long N, n;
    int i;
    long long digts;
    long long flag;
    long long temp;
    long long a = 1;
    while (~scanf("%lld", &N))
    {
        if (N <= 2)
        {
            printf("2\n");
            continue;
        }
        I[0][0] = 1;
        I[0][1] = 0;
        I[0][2] = 0;
        I[1][0] = 0;
        I[1][1] = 1;
        I[1][2] = 0;
        I[2][0] = 0;
        I[2][1] = 0;
        I[2][2] = 1;
        n = N - 2;
        digts = floor(log(n*1.0) / log(2.0)) + 1;
        for (i = 1; i <= digts; i++)
        {
            temp = a << (digts - i);
            flag = n / temp;
            n -= flag*temp;
            matMul(flag);
        }
        printf("%lld\n", (2*I[0][0] + 2*I[0][1]+I[0][2]) % con);
    }
    return 0;
}