28 if (R == 1) { zeta [0] = (beta-alpha)/(alpha+beta+2);
return; }
29 using namespace Eigen;
30 Matrix<T,Dynamic,1>
diag(R), subdiag(R-1);
31 diag[0] = (beta-alpha)/(alpha+beta+2);
32 for (
size_t r = 1; r < R; r++) {
33 diag[r] = (sqr(beta) - sqr(alpha))/((alpha+beta+2*r)*(alpha+beta+2*r+2));
35 subdiag[0] = 2/(alpha+beta+2)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
36 for (
size_t r = 2; r < R; r++) {
37 subdiag[r-1] = 2/(alpha+beta+2*r)
38 *sqrt(r*(alpha+r)*(beta+r)*(alpha+beta+r)/((alpha+beta+2*r-1)*(alpha+beta+2*r+1)));
40 SelfAdjointEigenSolver<Matrix<T,Dynamic,1> > es;
41 es.computeFromTridiagonal (
diag, subdiag, EigenvaluesOnly);
42 for (
size_t i = 0; i < R; ++i) zeta[i] = es.eigenvalues()(i);
43 if (alpha == beta && R%2 == 1) zeta[R/2] = 0.;