Support faster Cholesky Decomposition of Toeplitz matrices

I recently discovered that cholesky decomposition of toeplitz matrices (useful for GPs evaluated on a grid) can be done with O(n^2) rather than the standard cholesky O(n^3). I found an implementation here, where an academic has posted C, C++, Fortran, MATLAB & Python code (also mirrored in separate repositories here). I took a look at the python and it wasn't particularly optimized (used loops instead of vectorized alternatives), so in case it's useful here's the vectorized version I came up with:

import numpy as np

def cholesky_toeplitz(r):
    """
    Cholesky decomposition of a Toeplitz matrix.
    :param r: the first row of the Toeplitz matrix
    :return: the lower triangular matrix L
    """
    n = len(r)
    l = np.zeros([n, n])
    l[:, 0] = r
    g = np.zeros([2, n])
    g[0, 1:] = r[0:(n-1)]
    g[1, 1:] = r[1:]
    for i in range ( 1, n ):
        rho = - g[1,i] / g[0,i]
        gam = np.sqrt ( ( 1.0 - rho ) * ( 1.0 + rho ) )
        alf = g[0, i:n]
        bet = g[1, i:n]
        temp_g_0 = (alf + rho * bet) / gam
        temp_g_1 = (rho * alf + bet) / gam
        g[0, i:n] = temp_g_0
        g[1, i:n] = temp_g_1
        l[i:n, i] = g[0, i:n]
        g[0, i+1:n] = g[0, i:n-1]
        g[0,i] = 0.0
    return(l)

And in some light testing, it does seem to outperform scipy.linalg.cholesky() when n is large (& n.b. scipy uses LAPACK, so my comparison was native python vs compiled fortran).

Start with the linked Toeplitz Cholesky implementation and the provided cholesky_toeplitz function; compare its assumptions and numerical behavior with scipy.linalg.cholesky(). Identify the appropriate Stan Math entry point and tests, then verify that the supported Toeplitz decomposition has the intended performance and correctness.