<p>This is a very interesting sequence. It is almost but not quite the order-4 Fibonacci (a.k.a. Tetranacci) numbers. Having extracted the <a href="https://stackoverflow.com/a/40964862/6732794">doubling formulas for Tetranacci</a> from its companion matrix, I could not resist doing it again for this very similar recurrence relation. </p>

<p>Before we get into the actual code, some definitions and a short derivation of the formulas used are in order. Define an integer sequence <code>A</code> such that:</p>

<pre><code>A(n) := A(n-1) + A(n-3) + A(n-4)
</code></pre>

<p>with initial values <code>A(0), A(1), A(2), A(3) := 1, 1, 1, 2</code>.</p>

<p>For <code>n &gt;= 0</code>, this is the number of <a href="https://en.wikipedia.org/wiki/Composition_(combinatorics)" rel="nofollow noreferrer">integer compositions</a> of <code>n</code> into parts from the set <code>{1, 3, 4}</code>. This is the sequence that we ultimately wish to compute.</p>

<p>For convenience, define a sequence <code>T</code> such that:</p>

<pre><code>T(n) := T(n-1) + T(n-3) + T(n-4)
</code></pre>

<p>with initial values <code>T(0), T(1), T(2), T(3) := 0, 0, 0, 1</code>.</p>

<p>Note that <code>A(n)</code> and <code>T(n)</code> are simply shifts of each other. More precisely, <code>A(n) = T(n+3)</code> for all integers <code>n</code>. Accordingly, as elaborated by <a href="https://stackoverflow.com/a/41112123/6732794">another answer</a>, the companion matrix for both sequences is:</p>

<pre><code>[0  1  0  0]
[0  0  1  0]
[0  0  0  1]
[1  1  0  1]
</code></pre>

<p>Call this matrix <code>C</code>, and let:</p>

<pre><code>a, b, c, d := T(n), T(n+1), T(n+2), T(n+3)

a', b', c', d' := T(2n), T(2n+1), T(2n+2), T(2n+3)
</code></pre>

<p>By induction, it can easily be shown that:</p>

<pre><code>[0  1  0  0]^n = [d-c-a  c-b  b-a  a]
[0  0  1  0]     [  a    d-c  c-b  b]
[0  0  0  1]     [  b    b+a  d-c  c]
[1  1  0  1]     [  c    c+b  b+a  d]
</code></pre>

<p>As seen above, for any <code>n</code>, <code>C^n</code> can be fully determined from its rightmost column alone. Furthermore, multiplying <code>C^n</code> with its rightmost column produces the rightmost column of <code>C^(2n)</code>:</p>

<pre><code>[d-c-a  c-b  b-a  a][a] = [a'] = [a(2d - 2c - a) + b(2c - b)]
[  a    d-c  c-b  b][b]   [b']   [     a^2 + c^2 + 2b(d - c)]
[  b    b+a  d-c  c][c]   [c']   [     b(2a + b) + c(2d - c)]
[  c    c+b  b+a  d][d]   [d']   [     b^2 + d^2 + 2c(a + b)]
</code></pre>

<p>Thus, if we wish to compute <code>C^n</code> for some <code>n</code> by repeated squaring, we need only perform matrix-vector multiplication per step instead of the full matrix-matrix multiplication. </p>

<hr>

<p>Now, the implementation, in Python:</p>

<pre class="lang-python prettyprint-override"><code># O(n) integer additions or subtractions
def A_linearly(n): 
    a, b, c, d = 0, 0, 0, 1 # T(0), T(1), T(2), T(3)

    if n &gt;= 0:
        for _ in range(+n):
            a, b, c, d = b, c, d, a + b + d
    else: # n &lt; 0
        for _ in range(-n):
            a, b, c, d = d - c - a, a, b, c

    return d # because A(n) = T(n+3)

# O(log n) integer multiplications, additions, subtractions.
def A_by_doubling(n):
    n += 3 # because A(n) = T(n+3)

    if n &gt;= 0:
        a, b, c, d = 0, 0, 0, 1 # T(0), T(1), T(2), T(3)
    else: # n &lt; 0
        a, b, c, d = 1, 0, 0, 0 # T(-1), T(0), T(1), T(2)

    # Unroll the final iteration to avoid computing extraneous values
    for i in reversed(range(1, abs(n).bit_length())):
        w = a*(2*(d - c) - a) + b*(2*c - b)
        x = a*a + c*c + 2*b*(d - c)
        y = b*(2*a + b) + c*(2*d - c)
        z = b*b + d*d + 2*c*(a + b)

        if (n &gt;&gt; i) &amp; 1 == 0:
            a, b, c, d = w, x, y, z
        else: # (n &gt;&gt; i) &amp; 1 == 1
            a, b, c, d = x, y, z, w + x + z

    if n &amp; 1 == 0:
        return a*(2*(d - c) - a) + b*(2*c - b) # w
    else: # n &amp; 1 == 1
        return a*a + c*c + 2*b*(d - c)         # x


print(all(A_linearly(n) == A_by_doubling(n) for n in range(-1000, 1001)))
</code></pre>

<p>Because it was rather trivial to code, the sequence is extended to negative <code>n</code> in the usual way. Also provided is a simple linear implementation to serve as a point of reference. </p>

<p>For <code>n</code> large enough, the logarithmic implementation above is 10-20x faster than directly exponentiating the companion matrix with <code>numpy</code>, by a simple (i.e. not rigorous, and likely flawed) timing comparison. And by my estimate, it would still take ~100 years to compute <code>A(10**12)</code>! Even though the algorithm above has room for improvement, that number is simply too large. On the other hand, computing <code>A(10**12) mod M</code> for some <code>M</code> is much more attainable. </p>

<hr>

<h3>A direct relation to Lucas and Fibonacci numbers</h3>

<p>It turns out that <code>T(n)</code> is even closer to the Fibonacci and <a href="https://en.wikipedia.org/wiki/Lucas_number" rel="nofollow noreferrer">Lucas numbers</a> than it is to Tetranacci. To see this, note that the characteristic polynomial for <code>T(n)</code> is <code>x^4 - x^3 - x - 1 = 0</code> which factors into <code>(x^2 - x - 1)(x^2 + 1) = 0</code>. The first factor is the characteristic polynomial for Fibonacci &amp; Lucas! The 4 roots of <code>(x^2 - x - 1)(x^2 + 1) = 0</code> are the two Fibonacci roots, <code>phi</code> and <code>psi = 1 - phi</code>, and <code>i</code> and <code>-i</code>--the two square roots of <code>-1</code>.</p>

<p>The closed-form expression or "Binet" formula for <code>T(n)</code> will have the general form:</p>

<pre><code>T(n) = U(n) + V(n)
U(n) = p*(phi^n) + q*(psi^n)
V(n) = r*(i^n) + s*(-i)^n
</code></pre>

<p>for some constant coefficients <code>p, q, r, s</code>. </p>

<p>Using the initial values for <code>T(n)</code>, solving for the coefficients, applying some algebra, and noting that the Lucas numbers have the closed-form expression: <code>L(n) = phi^n + psi^n</code>, we can derive the following relations:</p>

<pre><code>       L(n+1) - L(n)    L(n-1)   F(n) + F(n-2)
U(n) = ------------- = -------- = ------------
             5            5           5
</code></pre>

<p>where <code>L(n)</code> is the n'th Lucas number with <code>L(0), L(1) := 2, 1</code> and <code>F(n)</code> is the n'th Fibonacci number with <code>F(0), F(1) := 0, 1</code>. And we also have:</p>

<pre><code>V(n) =  1 / 5   if n = 0 (mod 4)
     | -2 / 5   if n = 1 (mod 4)
     | -1 / 5   if n = 2 (mod 4)
     |  2 / 5   if n = 3 (mod 4)
</code></pre>

<p>Which is ugly, but trivial to code. Note that the numerator of <code>V(n)</code> <a href="https://math.stackexchange.com/q/2057991/368127">can also be succinctly expressed</a> as <code>cos(n*pi/2) - 2sin(n*pi/2)</code> or <code>(3-(-1)^n) / 2 * (-1)^(n(n+1)/2)</code>, but we use the piece-wise definition for clarity.</p>

<p>Here's an even nicer, more direct identity:</p>

<pre><code>T(n) + T(n+2) = F(n)
</code></pre>

<p>Essentially, we can compute <code>T(n)</code> (and therefore <code>A(n)</code>) by using Fibonacci &amp; Lucas numbers. Theoretically, this should be much more efficient than the Tetranacci-like approach.</p>

<p>It is known that the Lucas numbers can computed more efficiently than Fibonacci, therefore we will compute <code>A(n)</code> from the Lucas numbers. The most efficient, simple Lucas number algorithm I know of is one by L.F. Johnson (see his <a href="https://arxiv.org/abs/1012.0284" rel="nofollow noreferrer">2010 paper</a>: <em>Middle and Ripple, fast simple O(lg n) algorithms for Lucas Numbers</em>). Once we have a Lucas algorithm, we use the identity: <code>T(n) = L(n - 1) / 5 + V(n)</code> to compute <code>A(n)</code>. </p>

<pre class="lang-python prettyprint-override"><code># O(log n) integer multiplications, additions, subtractions
def A_by_lucas(n):
    n += 3 # because A(n) = T(n+3)
    offset = (+1, -2, -1, +2)[n % 4]
    L = lf_johnson_2010_middle(n - 1)
    return (L + offset) // 5

def lf_johnson_2010_middle(n):
    "-&gt; n'th Lucas number. See [L.F. Johnson 2010a]."
    #: The following Lucas identities are used:
    #:
    #:      L(2n)   = L(n)^2 - 2*(-1)^n
    #:      L(2n+1) = L(2n+2) - L(2n)
    #:      L(2n+2) = L(n+1)^2 - 2*(-1)^(n+1)
    #:
    #: The first and last identities are equivalent.
    #: For the unrolled iteration, the following is also used:
    #:
    #:      L(2n+1) = L(n)*L(n+1) - (-1)^n
    #:
    #: Since this approach uses only square multiplications per loop,
    #: It turns out to be slightly faster than standard Lucas doubling,
    #: which uses 1 square and 1 regular multiplication.
    if n &gt;= 0:
        a, b, sign = 2, 1, +1  # L(0), L(1), (-1)^0
    else: # n &lt; 0
        a, b, sign = -1, 2, -1 # L(-1), L(0), (-1)^(-1)

    # unroll the last iteration to avoid computing unnecessary values
    for i in reversed(range(1, abs(n).bit_length())):
        a = a*a - 2*sign # L(2k)
        c = b*b + 2*sign # L(2k+2)
        b = c - a        # L(2k+1)
        sign = +1

        if (n &gt;&gt; i) &amp; 1:
            a, b = b, c
            sign = -1

    if n &amp; 1:
        return a*b - sign
    else:
        return a*a - 2*sign
</code></pre>

<p>You may verify that <code>A_by_lucas</code> produces the same results as the previous <code>A_by_doubling</code> function, but is roughly 5x faster. Still not fast enough to compute <code>A(10**12)</code> in any reasonable amount of time! </p>
Using matrices to find the number of different ways to write n as the sum of 1, 3, and 4?
<p>This is a question given in this presentation. <a href="https://web.stanford.edu/class/cs97si/04-dynamic-programming.pdf" rel="nofollow noreferrer">Dynamic Programming</a></p>

<p>now i have implemented the algorithm using recursion and it works fine for small values. But when n is greater than 30 it becomes really slow.The presentation mentions that for  large values of n one should consider something similar to
<a href="https://math.stackexchange.com/questions/784710/how-to-prove-fibonacci-sequence-with-matrices">the matrix form of Fibonacci numbers</a> .I am having trouble undestanding how to use the matrix form of Fibonacci numbers to come up with a solution.Can some one give me some hints or pseudocode</p>

<p>Thanks</p>