https://www.hackerearth.com/practice/math/number-theory/basic-number-theory-1/practice-problems/algorithm/riyas-birthday-party-1/

Summary of problem statement

A sequence $a$ is generated as follows:

• $a_1 = 1$
• $a_2 = 6$
• $a_i = (a_{i-1} + a_{i+1})/2 – 2, \forall i \geq 2$.

Given $n$ ($1 \leq n \leq 10^{18}$), calculate $a_n \mod (10^9+ 7)$.

O(n) solution

Using properties of equalities, we convert: $a_i=(a_{i-1} + a_{i+1})/2 – 2$ to:
$a_{i+1} = 2a_i – a_{i-1} + 4$.

From this formula we derive a trivial $O(N)$ algorithm to calculate $a_n$:

typedef long long ll;

int Compute(ll n) {
    vector<int> a(n + 1);
    a[1] = 1;
    a[2] = 6;
    for (int i = 3; i < n; ++i) {
        a[i + 1] = ((ll) 2 * a[i] - a[i - 1] + 4) % 1000000007;
    }
    return a[n];
}

O(log n) solution

To speed up, we should think about matrix exponentiation (similar to the way we calculate the k^th Fibonacci number).

First, the matrix representation of $a$ is:

$\begin{aligned}\left[ \begin{array}{c} a_{i+1} \\ a_i \\ 1 \end{array} \right] = \left[ \begin{array}{ccc} 2 & -1 & 4 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array} \right] . \left[ \begin{array}{c} a_i \\ a_{i-1} \\ 1 \end{array} \right] \end{aligned}$

Let’s put in exponentiation:

$\left[ \begin{array}{c} a_{i+1} \\ a_i \\ 1 \end{array} \right] = \left[ \begin{array}{ccc} 2 & -1 & 4 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array} \right]^{i-1} . \left[ \begin{array}{c} a_2 \\ a_1 \\ 1 \end{array} \right]$

Substitute $i$ with $n$, $a_2$ with $6$, and $a_1$ with $1$:

$\left[ \begin{array}{c} a_{n+1} \\ a_n \\ 1 \end{array} \right] = \left[ \begin{array}{ccc} 2 & -1 & 4 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array} \right]^{n-1} . \left[ \begin{array}{c} 6 \\ 1 \\ 1 \end{array} \right]$

Matrix exponentiation $A^n$ can be done using $\log n$ multiplications. Hence it brings down the total time complexity to $O(\log n)$. Please find sample code below:

#define MOD 1000000007
typedef vector<vector<int>> Matrix;

Matrix MatrixMul(const Matrix& a, const Matrix& b) {
    int m = a.size(), n = a[0].size(), k = b[0].size();
    assert(n == b.size());

    Matrix c(m, vector<int>(k, 0));
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < k; ++j) {
            for (int t = 0; t < n; ++t) {
                c[i][j] += ((ll) a[i][t] * b[t][j]) % MOD;
                c[i][j] %= MOD;
            }}}
    return c;
}

Matrix MatrixPow(const Matrix& a, ll n) {
    if (n == 1) return a;
    Matrix tmp = MatrixPow(a, n/2);
    if (n % 2 == 0) return MatrixMul(tmp, tmp);
    else return MatrixMul(MatrixMul(tmp, tmp), a);
}

int Compute(int n) {
    if (n == 1) return 1;
    Matrix A = { {2, MOD - 1, 4}, {1, 0, 0}, {0, 0, 1} };
    Matrix P = { {6}, {1}, {1} };
    Matrix res = MatrixMul(MatrixPow(A, n - 1), P);
    return res[1][0];
}