Do you have a link to the maths vingette?

My implementation looks like:

```
vector diagSPD_exp(real alpha, real rho, real L, int m) {
vector[m] lambda = (linspaced_vector(m, 1, m)*pi()/(2*L))^2;
vector[m] diag = 2*alpha*rho/(1+rho^2*sqrt(lambda)^2);
diag = diag.^(0.5);
return diag;
}
```

testing both with `alpha = 1`

, `rho = 0.5`

, `L = 1`

, `m = 8`

I got:

`[0.847367,0.748398,0.677581,0.623686,0.580895,0.545854,0.516474,0.491379] `

(yours)

vs

`[0.786439,0.537029,0.390683,0.303314,0.246773,0.207584,0.178955,0.157177] `

(mine)

I’m actively using this code in my work so I want to make sure I didn’t make any mistakes in my derivation :)

Edit: I have narrowed down the difference to how the \lambda^2 term is coded. E.g. modifying:

```
vector diagSPD_Matern12(real alpha, real rho, real L, int M) {
vector[M] indices = linspaced_vector(M, 1, M);
real factor = 2;
vector[M] denom = rho * ((1 / rho)^2 + (pi() / 2 / L)^2 * indices^2);
return alpha * sqrt(factor * inv(denom));
}
```

gives identical

`[0.786439,0.537029,0.390683,0.303314,0.246773,0.207584,0.178955,0.157177] `

(or vice vesa if you change `sqrt(lambda)^2`

to `sqrt(lambda)`

in the other function)

That being said, I need to dig a bit deeper to see which is correct.