5 #include "tnt_array1d.h"
6 #include "tnt_array2d.h"
7 #include "tnt_math_utils.h"
110 for (
int j = 0; j < n; j++) {
116 for (
int i = n-1; i > 0; i--) {
122 for (
int k = 0; k < i; k++) {
123 scale = scale + abs(d[k]);
127 for (
int j = 0; j < i; j++) {
136 for (
int k = 0; k < i; k++) {
148 for (
int j = 0; j < i; j++) {
154 for (
int j = 0; j < i; j++) {
157 g = e[j] + V[j][j] * f;
158 for (
int k = j+1; k <= i-1; k++) {
165 for (
int j = 0; j < i; j++) {
169 Real hh = f / (h + h);
170 for (
int j = 0; j < i; j++) {
173 for (
int j = 0; j < i; j++) {
176 for (
int k = j; k <= i-1; k++) {
177 V[k][j] -= (f * e[k] + g * d[k]);
188 for (
int i = 0; i < n-1; i++) {
193 for (
int k = 0; k <= i; k++) {
194 d[k] = V[k][i+1] / h;
196 for (
int j = 0; j <= i; j++) {
198 for (
int k = 0; k <= i; k++) {
199 g += V[k][i+1] * V[k][j];
201 for (
int k = 0; k <= i; k++) {
206 for (
int k = 0; k <= i; k++) {
210 for (
int j = 0; j < n; j++) {
227 for (
int i = 1; i < n; i++) {
234 Real eps = pow(2.0,-52.0);
235 for (
int l = 0; l < n; l++) {
239 tst1 = max(tst1,abs(d[l]) + abs(e[l]));
244 if (abs(e[m]) <= eps*tst1) {
262 Real p = (d[l+1] - g) / (2.0 * e[l]);
263 Real r = hypot(p,1.0);
267 d[l] = e[l] / (p + r);
268 d[l+1] = e[l] * (p + r);
271 for (
int i = l+2; i < n; i++) {
285 for (
int i = m-1; i >= l; i--) {
295 p = c * d[i] - s * g;
296 d[i+1] = h + s * (c * g + s * d[i]);
300 for (
int k = 0; k < n; k++) {
302 V[k][i+1] = s * V[k][i] + c * h;
303 V[k][i] = c * V[k][i] - s * h;
306 p = -s * s2 * c3 * el1 * e[l] / dl1;
312 }
while (abs(e[l]) > eps*tst1);
320 for (
int i = 0; i < n-1; i++) {
323 for (
int j = i+1; j < n; j++) {
332 for (
int j = 0; j < n; j++) {
353 for (
int m = low+1; m <= high-1; m++) {
358 for (
int i = m; i <= high; i++) {
359 scale = scale + abs(H[i][m-1]);
366 for (
int i = high; i >= m; i--) {
367 ort[i] = H[i][m-1]/scale;
368 h += ort[i] * ort[i];
380 for (
int j = m; j < n; j++) {
382 for (
int i = high; i >= m; i--) {
386 for (
int i = m; i <= high; i++) {
391 for (
int i = 0; i <= high; i++) {
393 for (
int j = high; j >= m; j--) {
397 for (
int j = m; j <= high; j++) {
401 ort[m] = scale*ort[m];
408 for (
int i = 0; i < n; i++) {
409 for (
int j = 0; j < n; j++) {
410 V[i][j] = (i == j ? 1.0 : 0.0);
414 for (
int m = high-1; m >= low+1; m--) {
415 if (H[m][m-1] != 0.0) {
416 for (
int i = m+1; i <= high; i++) {
419 for (
int j = m; j <= high; j++) {
421 for (
int i = m; i <= high; i++) {
422 g += ort[i] * V[i][j];
425 g = (g / ort[m]) / H[m][m-1];
426 for (
int i = m; i <= high; i++) {
427 V[i][j] += g * ort[i];
438 void cdiv(Real xr, Real xi, Real yr, Real yi) {
440 if (abs(yr) > abs(yi)) {
443 cdivr = (xr + r*xi)/d;
444 cdivi = (xi - r*xr)/d;
448 cdivr = (r*xr + xi)/d;
449 cdivi = (r*xi - xr)/d;
469 Real eps = pow(2.0,-52.0);
471 Real p=0,q=0,r=0,s=0,z=0,t,w,x,y;
476 for (
int i = 0; i < nn; i++) {
477 if ((i < low) || (i > high)) {
481 for (
int j = max(i-1,0); j < nn; j++) {
482 norm = norm + abs(H[i][j]);
495 s = abs(H[l-1][l-1]) + abs(H[l][l]);
499 if (abs(H[l][l-1]) < eps * s) {
509 H[n][n] = H[n][n] + exshift;
517 }
else if (l == n-1) {
518 w = H[n][n-1] * H[n-1][n];
519 p = (H[n-1][n-1] - H[n][n]) / 2.0;
522 H[n][n] = H[n][n] + exshift;
523 H[n-1][n-1] = H[n-1][n-1] + exshift;
545 r = sqrt(p * p+q * q);
551 for (
int j = n-1; j < nn; j++) {
553 H[n-1][j] = q * z + p * H[n][j];
554 H[n][j] = q * H[n][j] - p * z;
559 for (
int i = 0; i <= n; i++) {
561 H[i][n-1] = q * z + p * H[i][n];
562 H[i][n] = q * H[i][n] - p * z;
567 for (
int i = low; i <= high; i++) {
569 V[i][n-1] = q * z + p * V[i][n];
570 V[i][n] = q * V[i][n] - p * z;
595 w = H[n][n-1] * H[n-1][n];
602 for (
int i = low; i <= n; i++) {
605 s = abs(H[n][n-1]) + abs(H[n-1][n-2]);
620 s = x - w / ((y - x) / 2.0 + s);
621 for (
int i = low; i <= n; i++) {
638 p = (r * s - w) / H[m+1][m] + H[m][m+1];
639 q = H[m+1][m+1] - z - r - s;
641 s = abs(p) + abs(q) + abs(r);
648 if (abs(H[m][m-1]) * (abs(q) + abs(r)) <
649 eps * (abs(p) * (abs(H[m-1][m-1]) + abs(z) +
650 abs(H[m+1][m+1])))) {
656 for (
int i = m+2; i <= n; i++) {
665 for (
int k = m; k <= n-1; k++) {
666 int notlast = (k != n-1);
670 r = (notlast ? H[k+2][k-1] : 0.0);
671 x = abs(p) + abs(q) + abs(r);
681 s = sqrt(p * p + q * q + r * r);
689 H[k][k-1] = -H[k][k-1];
700 for (
int j = k; j < nn; j++) {
701 p = H[k][j] + q * H[k+1][j];
703 p = p + r * H[k+2][j];
704 H[k+2][j] = H[k+2][j] - p * z;
706 H[k][j] = H[k][j] - p * x;
707 H[k+1][j] = H[k+1][j] - p * y;
712 for (
int i = 0; i <= min(n,k+3); i++) {
713 p = x * H[i][k] + y * H[i][k+1];
715 p = p + z * H[i][k+2];
716 H[i][k+2] = H[i][k+2] - p * r;
718 H[i][k] = H[i][k] - p;
719 H[i][k+1] = H[i][k+1] - p * q;
724 for (
int i = low; i <= high; i++) {
725 p = x * V[i][k] + y * V[i][k+1];
727 p = p + z * V[i][k+2];
728 V[i][k+2] = V[i][k+2] - p * r;
730 V[i][k] = V[i][k] - p;
731 V[i][k+1] = V[i][k+1] - p * q;
744 for (n = nn-1; n >= 0; n--) {
753 for (
int i = n-1; i >= 0; i--) {
756 for (
int j = l; j <= n; j++) {
757 r = r + H[i][j] * H[j][n];
768 H[i][n] = -r / (eps * norm);
776 q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
777 t = (x * s - z * r) / q;
779 if (abs(x) > abs(z)) {
780 H[i+1][n] = (-r - w * t) / x;
782 H[i+1][n] = (-s - y * t) / z;
789 if ((eps * t) * t > 1) {
790 for (
int j = i; j <= n; j++) {
791 H[j][n] = H[j][n] / t;
804 if (abs(H[n][n-1]) > abs(H[n-1][n])) {
805 H[n-1][n-1] = q / H[n][n-1];
806 H[n-1][n] = -(H[n][n] - p) / H[n][n-1];
808 cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q);
814 for (
int i = n-2; i >= 0; i--) {
818 for (
int j = l; j <= n; j++) {
819 ra = ra + H[i][j] * H[j][n-1];
820 sa = sa + H[i][j] * H[j][n];
840 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
841 vi = (d[i] - p) * 2.0 * q;
842 if ((vr == 0.0) && (vi == 0.0)) {
843 vr = eps * norm * (abs(w) + abs(q) +
844 abs(x) + abs(y) + abs(z));
846 cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
849 if (abs(x) > (abs(z) + abs(q))) {
850 H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x;
851 H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x;
853 cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q);
861 t = max(abs(H[i][n-1]),abs(H[i][n]));
862 if ((eps * t) * t > 1) {
863 for (
int j = i; j <= n; j++) {
864 H[j][n-1] = H[j][n-1] / t;
865 H[j][n] = H[j][n] / t;
875 for (
int i = 0; i < nn; i++) {
876 if (i < low || i > high) {
877 for (
int j = i; j < nn; j++) {
885 for (
int j = nn-1; j >= low; j--) {
886 for (
int i = low; i <= high; i++) {
888 for (
int k = low; k <= min(j,high); k++) {
889 z = z + V[i][k] * H[k][j];
910 for (
int j = 0; (j < n) && issymmetric; j++) {
911 for (
int i = 0; (i < n) && issymmetric; i++) {
912 issymmetric = (A[i][j] == A[j][i]);
917 for (
int i = 0; i < n; i++) {
918 for (
int j = 0; j < n; j++) {
933 for (
int j = 0; j < n; j++) {
934 for (
int i = 0; i < n; i++) {
1012 for (
int i = 0; i < n; i++) {
1013 for (
int j = 0; j < n; j++) {
1019 }
else if (e[i] < 0) {