C++ algorithm code for Magical sequence that will generate desired output

Question

The Magical Sequence

A Magical Sequence is defined as shown.

Magical[1] = 0
Magical[2] = 1
Magical[n] = Magical[n-1] + 2*Magical[n-2] + 3*Magical[n-3] + ... (n-1)*Magical[1] + n*1., for n > 2

Given n (1 <= n <= 10^9 ), find Magical[n].

Example 1: input: 3 Output: 4

Explanation:

Magical[n] = 1*Magical[n-1] + 2*Magical[n-2] + 3*1
    Magical[3] = 1*Magical[2] + 2*Magical[1] + 3*1
    Magical[3] = 1*1 + 2*0 + 3*1
    Magical[3] = 4

Example 2: input: 4 Output: 10

Magical[4] = 1*Magical[3]+2*Magical[2]+3*Magical[1]+4*1
          = 1*4+2*1+3*0+4 = 10

Example 3: input: 5 Output: 26

Magical[5] = 1*Magical[4]+2*Magical[3]+3*Magical[2]+4*Magical[1]+5*1
          = 1*10+2*4+3*1+4*0+5 = 26

I tried something like below :-

int CuckooNum(int n)
{
    if (1 == n)
    {
        return 0;
    }
    else if (2 == n)
    {
        return 1;
    }

    std::vector<int> vec;
    vec.resize(n);
    vec[0] = 4;
    vec[1] = 0;
    vec[2] = 1;

    int multiplyer = n;
    int result = 0;
    for (int index=3; index <= n; index++)
    {
        result += multiplyer * vec[index-1];
        vec[index] = result;
        multiplyer--;
    }
    return result;
}

Answer 1

As the size n can be very large (10^9), a direct implementation O(n^2) is not possible. A specific algorithm is needed. I will focus here on the algorithm, and propose a O(log n) solution.

To simplify explanation, I rename magical[] as x[]

Moreover, we can define x[0] = 1. Then,

x[n] = x[n-1] + 2*x[n-2] + 3*x[n-3] + ... (n-1)*x[1] + n*x[0]

As

x[n-1] =        1*x[n-2] + 2*x[n-3] + ... (n-2)*x[1] + (n-1)*x[0]

It follows

x[n] - x[n-1] = x[n-1] + x[n-2] + x[n-3] + ... x[1] + x[0] = S[n-1]

When S[n] represents the sum of the terms until n (x[0] included)

Moreover,

S[n] = S[n-1] + x[n] = 2*S[n-1] + x[n-1]

Therefore, the iterative formula can be represented in a simple matrix form:

(x[n]) = (1 1) (x[n-1])
(S[n])   (1 2) (S[n-1])

Or, defining the vector (x[n] S[n])^t as Z[n]:

Z[n] = A * Z[n-1] where A is the matrix (1 1)
                                        (1 2)

Note: this formula is valid for n>= 4 only, as the first x[n] values do no respect the simple recurrence relation.

It follows that

 Z[n] = A^(n-3) Z[3] with Z[3] = (4 6)^t

Classically, this calculation can be performed with O(log n) complexity, iteratively calculating A^2, A^4, A^8 etc.

Pay attention that the values increase rapidly.

Here is an example of C++ implementation. Note that this implementation is not optimized, as for example it doesn't use the fact that all matrices are symmetric.

#include <iostream>
#include <array>

using Matr22 = std::array<std::array<long long int, 2>, 2>;
using Vect2 = std::array<long long int, 2>;

Matr22 Matrsquare (const Matr22 &m) {
    Matr22 m2;
    m2[0][0] = m[0][0]*m[0][0] + m[0][1]*m[1][0];
    m2[0][1] = m[0][0]*m[0][1] + m[0][1]*m[1][1];
    m2[1][0] = m[1][0]*m[0][0] + m[1][1]*m[1][0];
    m2[1][1] = m[1][0]*m[0][1] + m[1][1]*m[1][1];
    
    return m2;
}

Matr22 Mult (const Matr22 &m1, const Matr22 &m2) {
    Matr22 y;
    y[0][0] = m1[0][0]*m2[0][0] + m1[0][1]*m2[1][0];
    y[0][1] = m1[0][0]*m2[0][1] + m1[0][1]*m2[1][1];
    y[1][0] = m1[1][0]*m2[0][0] + m1[1][1]*m2[1][0];
    y[1][1] = m1[1][0]*m2[0][1] + m1[1][1]*m2[1][1];
    
    return y;
}

Vect2 Mult (const Matr22 &m, const Vect2& x) {
    Vect2 y;
    y[0] = m[0][0] * x[0] + m[0][1] * x[1];
    y[1] = m[1][0] * x[0] + m[1][1] * x[1];
    return y;
}

//  Matrix exponentiation
Matr22 Mult_exp (const Matr22 &m, int exp) {
    Matr22 y = {1, 0, 0, 1};
    if (exp == 0) return y;
    Matr22 M2k = m;
    while (exp) {
        if (exp%2) y = Mult (y, M2k);
        M2k = Matrsquare (M2k);
        exp /= 2;
    };
    return y;
}


long long int Magical (int n) {
    if (n == 1) return 0;
    if (n == 2) return 1;
    if (n == 3) return 4;
    Matr22 A = {1, 1, 1, 2};
    Vect2 z = {4, 6};       // corresponds to n=3
    auto Ak = Mult_exp (A, n-3);
    z = Mult (Ak, z);
    return z[0];
}

int main() {
    int n;
    std::cout << "Input n: ";
    std::cin >> n;
    auto ans = Magical (n);
    std::cout << "Magical[" << n << "] = " << ans << '\n';
}

Answer 2

long long func(int n)
{
  if (n==1) return 0;
  else if (n==2) return 1;
  else return 1*func(n-1)+2*func(n-2)+n;
}