The standard approach to the recursive solution has time complexity O(n), with the required memoization adding a space complexity of O(n).

Using identities derived through the matrix form, we can produce a faster algorithm that scales much better.

Here are the required identities:

The new time complexity for this algorithm is ~O (log( n )). (Technically, it is O ( M ( n ) log( n )) where M ( n ) is the time to multiply two n-digit numbers together, this can be found in the above linked Wikipedia article)

In addition, this optimized solution avoids running into the recursion depth in Python nearly as quickly. Where the normal solution can't even calculate up to n = 1,000, the optimized variation can compute well past n = 1,000,000 without running into the default recursion depth in Python.