Detailed mathematical approach

Problem Summary

Problem 940 asks for values of a two-dimensional recurrence \(A(m,n)\) at Fibonacci-indexed coordinates \((F_i,F_j)\), summed over all \(2 \le i,j \le 50\), with the final result reduced modulo \(1123581313\).

The implementations show that this surface is generated by two fixed \(2\times 2\) matrices. That matrix viewpoint is enough to reconstruct the actual recurrence, the boundary data, and an efficient evaluation strategy for indices as large as \(F_{50}=12586269025\) without ever filling a gigantic grid.

Mathematical Approach

Introduce the matrices and vectors

$C=\begin{pmatrix}0&1\\3&1\end{pmatrix},\qquad M=\begin{pmatrix}1&1\\3&2\end{pmatrix},\qquad v_0=\begin{pmatrix}0\\1\end{pmatrix},\qquad e_1=\begin{pmatrix}1\\0\end{pmatrix}.$

The scalar quantity being sampled is

$A(m,n)=e_1^{\mathsf T} M^m C^n v_0.$

Everything important in the problem comes from the relations between \(C\) and \(M\).

The boundary sequence in the \(n\)-direction

Start with the row \(m=0\). If we write

$u_n=A(0,n)=e_1^{\mathsf T} C^n v_0,$

then \(u_n\) is just the first component of \(C^n v_0\). Let \(C^n v_0=(x_n,y_n)^{\mathsf T}\). Multiplying once more by \(C\) gives

$x_{n+1}=y_n,\qquad y_{n+1}=3x_n+y_n,$

so the first component satisfies

$u_{n+2}=u_{n+1}+3u_n,\qquad u_0=0,\qquad u_1=1.$

The first values are therefore

$u_0=0,\ u_1=1,\ u_2=1,\ u_3=4,\ u_4=7,\ u_5=19,\dots$

This single one-dimensional recurrence is the boundary from which the whole two-dimensional array can be generated.

Why one step in \(m\) adds the next value in \(n\)

The key identity is

$M=I+C.$

Using that directly inside the definition of \(A(m,n)\),

$A(m+1,n)=e_1^{\mathsf T} M^m (I+C) C^n v_0=A(m,n)+A(m,n+1).$

This is the genuinely two-dimensional recurrence. Moving one step in the \(m\)-direction means taking the current value and adding the value one column further in the \(n\)-direction. Once the boundary row \(A(0,n)\) is known, repeated use of this rule fills every higher row.

Recurrences in both directions

The matrix \(C\) satisfies its characteristic equation

$C^2=C+3I.$

Multiplying on the left by \(e_1^{\mathsf T}M^m\) and on the right by \(C^n v_0\) gives a vertical recurrence valid in every row:

$A(m,n+2)=A(m,n+1)+3A(m,n).$

Likewise \(M\) has characteristic polynomial \(x^2-3x-1\), so

$M^2=3M+I,$

and therefore

$A(m+2,n)=3A(m+1,n)+A(m,n).$

The surface is thus anisotropic but perfectly consistent: one recurrence governs vertical motion and another governs horizontal motion.

A binomial-transform interpretation

Because \(M=I+C\) is a polynomial in \(C\), the two matrices commute. Expanding \((I+C)^m\) shows

$A(m,n)=e_1^{\mathsf T}(I+C)^m C^n v_0=\sum_{t=0}^m \binom{m}{t} e_1^{\mathsf T} C^{n+t} v_0=\sum_{t=0}^m \binom{m}{t} u_{n+t}.$

So every row is a binomial transform of the boundary sequence \(u_n\). The implementations do not explicitly expand this formula, but it explains why the matrix model and the recurrence model describe exactly the same object.

Worked example: the first few entries

Using \(u_0=0\), \(u_1=1\), \(u_2=1\), \(u_3=4\), the first rows become

$ \begin{array}{c|cccc} A(m,n) & n=0 & n=1 & n=2 & n=3 \\ \hline m=0 & 0 & 1 & 1 & 4 \\ m=1 & 1 & 2 & 5 & 11 \\ m=2 & 3 & 7 & 16 & 37 \end{array} $

For instance, \(A(2,2)=16\) can be verified in two independent ways. From the horizontal rule,

$A(2,2)=A(1,2)+A(1,3)=5+11=16.$

From the vertical rule,

$A(2,2)=A(2,1)+3A(2,0)=7+3\cdot 3=16.$

This is a compact check that both recurrences are describing the same table.

Fibonacci sampling

Let \(F_0=0\), \(F_1=1\), and \(F_{r+2}=F_{r+1}+F_r\). The target quantity is

$S(k)=\sum_{i=2}^k \sum_{j=2}^k A(F_i,F_j),$

with the final answer being \(S(50)\pmod{1123581313}\).

These indices grow very quickly, so it is not practical to build the grid cell by cell up to \(m,n=F_{50}\). Matrix exponentiation reaches \(M^{F_i}\) and \(C^{F_j}\) in logarithmic time in the exponent, which is exactly what makes the computation feasible.

How the Code Works

The C++, Python, and Java implementations first generate the Fibonacci numbers up to the requested limit. They then apply binary exponentiation under the modulus \(1123581313\) to the two fixed \(2\times 2\) matrices \(M\) and \(C\).

For each Fibonacci index \(F_i\), the implementation stores two reusable objects: the full matrix \(M^{F_i}\) and the vector \(C^{F_i}v_0\). Once those have been prepared, a term of the double sum is only a first-row dot product:

$A(F_i,F_j)=\bigl(\text{first row of }M^{F_i}\bigr)\cdot\bigl(C^{F_j}v_0\bigr).$

This avoids repeating expensive exponentiations inside the final nested loop. The last pass therefore consists only of modular additions and multiplications on already cached 2-component data.

One implementation also checks small identities such as \(A(1,1)=2\), \(A(2,2)=16\), \(S(3)=30\), and \(S(5)=10396\), which match the recurrence table above.

Complexity Analysis

Generating the Fibonacci list costs \(O(k)\). The code then computes \(k-1\) powers of \(M\) and \(k-1\) powers of \(C\). Each power is obtained by exponentiation by squaring, so each one costs \(O(\log F_k)\) multiplications of \(2\times 2\) matrices, and each such multiplication has constant cost.

The preprocessing phase is therefore \(O(k\log F_k)\). The final accumulation over all ordered pairs \((i,j)\) with \(2\le i,j\le k\) is \(O(k^2)\). Memory usage is \(O(k)\) for the Fibonacci values, the cached matrices, and the cached vectors. For the actual target \(k=50\), these costs are very small.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=940
2. Fibonacci numbers: Wikipedia - Fibonacci number
3. Recurrence relations: Wikipedia - Recurrence relation
4. Exponentiation by squaring: Wikipedia - Exponentiation by squaring
5. Companion matrices: Wikipedia - Companion matrix