text: #include <cmath>
  #include "nr.h"
  using namespace std;
  
  void NR::tqli(Vec_IO_DP &d, Vec_IO_DP &e, Mat_IO_DP &z)
  {
  	int m,l,iter,i,k;
  	DP s,r,p,g,f,dd,c,b;
  
  	int n=d.size();
  	for (i=1;i<n;i++) e[i-1]=e[i];
  	e[n-1]=0.0;
  	for (l=0;l<n;l++) {
  		iter=0;
  		do {
  			for (m=l;m<n-1;m++) {
  				dd=fabs(d[m])+fabs(d[m+1]);
  				if (fabs(e[m])+dd == dd) break;
  			}
  			if (m != l) {
  				if (iter++ == 30) nrerror("Too many iterations in tqli");
  				g=(d[l+1]-d[l])/(2.0*e[l]);
  				r=pythag(g,1.0);
  				g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
  				s=c=1.0;
  				p=0.0;
  				for (i=m-1;i>=l;i--) {
  					f=s*e[i];
  					b=c*e[i];
  					e[i+1]=(r=pythag(f,g));
  					if (r == 0.0) {
  						d[i+1] -= p;
  						e[m]=0.0;
  						break;
  					}
  					s=f/r;
  					c=g/r;
  					g=d[i+1]-p;
  					r=(d[i]-g)*s+2.0*c*b;
  					d[i+1]=g+(p=s*r);
  					g=c*r-b;
  					// Next loop can be omitted if eigenvectors not wanted
  					for (k=0;k<n;k++) {
  						f=z[k][i+1];
  						z[k][i+1]=s*z[k][i]+c*f;
  						z[k][i]=c*z[k][i]-s*f;
  					}
  				}
  				if (r == 0.0 && i >= l) continue;
  				d[l] -= p;
  				e[l]=g;
  				e[m]=0.0;
  			}
  		} while (m != l);
  	}
  }
  