function pythag(a, b){
	
	var absa,absb;
    absa=abs(a);
    absb=abs(b);
    if (absa > absb) return absa*sqrt(1.0+(absb/absa)*(absb/absa));
    else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+(absa/absb)*(absa/absb)));
	
}

function SIGN(a, b){
	return(b >= 0 ? abs(a) : -abs(a));
}


function TRED2(A, D, E){
	var N = array_length(A);
	var L, K, J, I, scale, HH, H, G, F;
	for(I = N - 1; I > 0; I--){
		L = I - 1;
		H = 0;
		scale = 0;
		if(L > 0){
			for(K = 0; K <= L; K++){ scale += abs(A[I][K]); }
			if(scale == 0){
				E[I] = A[I][L];
			}else{
				for(K = 0; K <= L; K++){
					A[I][K] /= scale;
					H += A[I][K] * A[I][K];
				}
				F = A[I][L];
				G = (F >= 0 ? -sqrt(H) : sqrt(H));
				E[I] = scale * G;
				H -= F * G;
				A[I][L] = F - G;
				F = 0;
				for(J = 0; J <= L; J++){
					G = 0;
					for(K = 0; K <= J; K++){ G += A[J][K] * A[I][K]; }
					for(K = J + 1; K <= L; K++){ G += A[K][J] * A[I][K]; }
					E[J] = G / H;
					F += E[J] * A[I][J];
				}
				HH = F / (H + H);
				for(J = 0; J <= L; J++){
					F = A[I][J];
					G = E[J] - HH * F;
					E[J] = G;
					for(K = 0; K <= J; K++){ A[J][K] -= (F * E[K] + G * A[I][K]); }
				}
			}
		}else{
			E[I] = A[I][L];
		}
		D[I] = H;
	}
	for(I = 0; I < N; I++){
		D[I] = A[I][I];
	}
}

function TQLI(D, E){
	var N = array_length(D);
	var M, L, iter, I, K, S, R, P, G, F, DD, C, B;
	for(I = 1; I < N; I++){
		E[I - 1] = E[I];
	}
	E[N - 1] = 0;
	for(L = 0; L < N; L++){
		iter = 0;
		do{
			for(M = L; M < N - 1; M++){
				DD = abs(D[M]) + abs(D[M + 1]);
				if(abs(E[M]) + DD == DD){ break; }
			}
			if(M != L){
				if(iter++ == 30){ show_debug_message("TOO MANY ITERATIONS!"); }
				G = (D[L + 1] - D[L])/(2 * E[L]);
				R = pythag(G, 1);
				G = D[M] - D[L] + E[L]/(G + SIGN(R, G));
				S = 1;
				C = 1;
				P = 0;
				for(I = M - 1; I >= L; I--){
					F = S * E[I];
					B = C * E[I];
					R = pythag(F, G);
					E[I + 1] = R;
					if(R == 0){
						D[I + 1] -= P;
						E[M] = 0;
						break;
					}
					S = F/R;
					C = G/R;
					G = D[I + 1] - P;
					R = (D[I] - G) * S + 2 * C * B;
					P = S * R;
					D[I + 1] = G + P;
					G = C * R - B;
				}
				if((R == 0 and I >= L)) continue;
					D[L] -= P;
					E[L] = G;
					E[M] = 0;
				
			}
		}until(M == L);
	}
}