Tail recursive formulation of the Legendre polynomial relation The recursive formula for Legendre polynomials is widely known:
$(n + 1) P_{n+1}(x) 
= (2n + 1) x P_{n}(x) - n P_{n-1}(x).$
Let us rewrite the above as follows for convenience:
$P_{n}(x) 
= \frac{2n - 1}{n} x P_{n-1}(x) - \frac{n - 1}{n} P_{n-2}(x).$
If we assume a fixed $x$ we can rewrite the expression in the following form:
$f(n) 
= a(n) x f(n-1) - b(n) f(n-2),$
with $f(n) \equiv P_{n}(x),$ $a(n) \equiv \frac{2n - 1}{n}$ and $b(n) \equiv \frac{n - 1}{n}.$ The edge cases are $f(0) = 1$ and $f(1) = x.$
There are many ways of computing the Legendre polynomials with regular sums, most of which stem from the Rodriguez formula. Numerically this would entail regular loops or sums of lists. I am curious, however, if the above "regular"-recursive definition for $f(n)$ can be rewritten tail-recursively à la Fibonacci sequence to get a reasonably efficient recursive numerical routine.

In the process of writing this question I believe I have drawn up a working proof, so as per Meta suggestions I added the solution-verification tag and added my work as an answer
 A: Let the $\,2\times2\,$ matrix
$$ A_k(x):=\begin{bmatrix} \frac{2k+1}kx & \frac1k \\ -k &0 \end{bmatrix}. \tag{1}$$
Then if $\,n>0,\,$ $\,P_n(x)\,$ is derived from a matrix product
$$ P_n(x) = [x\;\; 1]\left(\prod_{k=1}^{n-1}A_k(x)\right)
\begin{bmatrix} \frac1n\\  0\end{bmatrix}. \tag{2}$$ There are several ways to compute the matrix product including
tail recursively. The analogous version for the Fibonacci sequence is
$$ F_n = [1\;\; 0]\left(\prod_{k=1}^{n-1}
 \begin{bmatrix} 1&1 \\ 1&0 \end{bmatrix} \right)\begin{bmatrix} 1\\  0\end{bmatrix}.$$
A: If this works the way it does for the Fibonacci sequence, we can define a helper function with a counter integer and two accumulator values:
$f(k) 
= g\big(k \, \big| \, f(0), f(1)\big) 
= g\big(k \, \big| \, 1, x\big),$
and the edge case for $g$ would be
$g\big(0 \, \big| \, f(n), f(n+1)\big) = f(n), 
\quad \text{or} 
\quad g\big(1 \, \big| \, f(n), f(n+1)\big) = f(n+1).$
For normal behavior we can work our way up from $g\big(2 \, \big| \, 1, x\big)$:
$g\big(2 \, \big| \, 1, x\big) 
= g\big(1 \, \big| \, x, a(2) x \cdot x - b(2) \cdot 1\big),$
from which we can deduce
$g\big(k \, \big| \, f(n), f(n+1)\big) 
= g\big(k - 1 \, \big| \, f(n+1), a(k) x f(n+1) - b(k) f(n)\big),$
which is nice and $\mathcal{O}(k).$
