Bremsstrahlung

I really wish I'd left myself a few more comments about how this code works:

void trid_inplace_lopt3(vector<double>& x){
    const unsigned int numL=16;
    static double l[numL]={1./4.,4./15.,15/56.,56./209.,209./780.,780./2911.,2911./10864,
    10864./40545.,40545./151316.,151316./564719.,564719./2107560.,2107560./7865521.,
    7865521./29354524.,29354524./109552575.,109552575./408855776.,408855776./1525870529.};
    const double lLim=2.-sqrt(3.);

    const unsigned int N=x.size();
    unsigned int i=1;
    if(N<=numL){
        for(; i<N; i++)
            x[i]=x[i]-l[i-1]*x[i-1];
        x[N-1]*=l[N-1];
        for(i=N-1; i>0; i--)
            x[i-1]=l[i-1]*(x[i-1]-x[i]);
    }
    else{
        for(; i<numL; i++)
            x[i]=x[i]-l[i-1]*x[i-1];
        for(; i<N; i++)
            x[i]=x[i]-lLim*x[i-1];
        x[N-1]*=lLim;
        for(i=N-1; i>numL; i--)
            x[i-1]=lLim*(x[i-1]-x[i]);
        for(; i>0; i--)
            x[i-1]=l[i-1]*(x[i-1]-x[i]);
    }
}

From context I'm pretty sure that it solves the tridiagonal matrix equation that appears in finding the coefficients of a cubic interpolating spline using an LU decomposition. I recall that at the time I wrote it I thought that the sequence of constants was pretty cool, and I still do, but of course I've lost the paper on which I derived it.