8 #define MAX(a, b) ((a)>(b)?(a):(b))
17 static double hypot2(
double x,
double y) {
23 static void tred2(
double V[
n][
n],
double d[
n],
double e[
n]) {
32 for (j = 0; j <
n; j++) {
38 for (i =
n-1; i > 0; i--) {
44 for (k = 0; k < i; k++) {
45 scale = scale + fabs(
d[k]);
49 for (j = 0; j < i; j++) {
58 for (k = 0; k < i; k++) {
70 for (j = 0; j < i; j++) {
76 for (j = 0; j < i; j++) {
79 g = e[j] + V[j][j] *
f;
80 for (k = j+1; k <= i-1; k++) {
87 for (j = 0; j < i; j++) {
92 for (j = 0; j < i; j++) {
95 for (j = 0; j < i; j++) {
98 for (k = j; k <= i-1; k++) {
99 V[k][j] -= (
f * e[k] + g *
d[k]);
110 for (i = 0; i <
n-1; i++) {
115 for (k = 0; k <= i; k++) {
116 d[k] = V[k][i+1] / h;
118 for (j = 0; j <= i; j++) {
120 for (k = 0; k <= i; k++) {
121 g += V[k][i+1] * V[k][j];
123 for (k = 0; k <= i; k++) {
128 for (k = 0; k <= i; k++) {
132 for (j = 0; j <
n; j++) {
142 static void tql2(
double V[
n][
n],
double d[
n],
double e[
n]) {
150 double g,p,r,dl1,h,
f,tst1,eps;
151 double c,c2,c3,el1,
s,s2;
153 for (i = 1; i <
n; i++) {
160 eps = pow(2.0,-52.0);
161 for (l = 0; l <
n; l++) {
165 tst1 =
MAX(tst1,fabs(
d[l]) + fabs(e[l]));
168 if (fabs(e[m]) <= eps*tst1) {
185 p = (
d[l+1] - g) / (2.0 * e[l]);
190 d[l] = e[l] / (p + r);
191 d[l+1] = e[l] * (p + r);
194 for (i = l+2; i <
n; i++) {
208 for (i = m-1; i >= l; i--) {
218 p = c *
d[i] -
s * g;
219 d[i+1] = h +
s * (c * g +
s *
d[i]);
223 for (k = 0; k <
n; k++) {
225 V[k][i+1] =
s * V[k][i] + c * h;
226 V[k][i] = c * V[k][i] -
s * h;
229 p = -
s * s2 * c3 * el1 * e[l] / dl1;
235 }
while (fabs(e[l]) > eps*tst1);
243 for (i = 0; i <
n-1; i++) {
246 for (j = i+1; j <
n; j++) {
255 for (j = 0; j <
n; j++) {
267 for (i = 0; i <
n; i++) {
268 for (j = 0; j <
n; j++) {