begin comment eigenvalues of a real symmetric matrix
     by the QR method. Algorithm 253, P.A. Businger, CACM 8 (1965) 217;
   integer n;
   n ≔ read;
   begin integer i,j;
      real array a[1:n,1:n];
      procedure symmetric QR1(n,g);
         value n; integer n; array g;
           comment uses Housholders.s method and the QR algorithm to
           find all n eigenvalues of the real symmetric matrix whose lower
           triangular part is given in the array g[1:n,1:n]. The computed
           eigenvalues are stored as the diagonal elements g[i,i]. The
           original contents of the lower triangular part of g are lost during
           the computation whereas the strictly upper triagular part of g
           is left untouched;
      begin
         real procedure sum(i,m,n,a);
            value m,n; integer i,m,n; real a;
         begin real s; s ≔ 0;
            for i ≔ m step 1 until n do s ≔ s + a;
            sum ≔ s
         end sum;
         real procedure max (a,b);
            value a,b; real a,b;
            max ≔ if a > b then a else b;
         procedure Housholder tridiagonalization 1(n,g,a,bq,norm);
            value n; integer n; array g,a,bq; real norm;
              comment nonlocal real procedure sum, max;
              comment reduces the given real symmetric n by n matrix g
              to tridiagonal form using n - 2 elementary orthogonal trans-
              formations (I-2ww.) = (I-gamma uu.). Only the lower tri
              angular part of g need be given. The diagonal elements and
              the squares of the subdiagonal elements of the reduced matrix
              are stored in a[1:n] and bq[1:n-1] respectively. norm is set
              equal to the infinity norm of the reduced matrix. The columns
              of the strictly lower triagular part of g are replaced by the
              nonzero portions of the vectors u;
         begin integer i,j,k;
            real t,absb,alpha,beta,gamma,sigma;
            array p[2:n];
            norm ≔ absb ≔ 0;
            for k ≔ 1 step 1 until n - 2 do
               begin a[k] ≔ g[k,k];
                  sigma ≔ bq[k] ≔ sum(i,k+1,n,g[i,k]⭡2);
                  t ≔ absb + abs(a[k]); absb ≔ sqrt(sigma);
                  norm ≔ max(norm,t+absb);
                  if sigma ≠ 0 then
                  begin alpha ≔ g[k+1,k];
                     beta ≔ if alpha < 0 then absb else - absb;
                     gamma ≔ 1 / (sigma-alpha×beta); g[k+1,k] ≔ alpha - beta;
                     for i ≔ k + 1 step 1 until n do
                          p[i] ≔ gamma ×
                          (sum(j,k+1,i,g[i,j]×g[j,k]) + sum(j,i+1,n,g[j,i]×g[j,k]));
                     t ≔ 0·5 × gamma × sum(i,k+1,n,g[i,k]×p[i]);
                     for i ≔ k + 1 step 1 until n do
                          p[i] ≔ p[i] - t×g[i,k];
                     for i ≔ k + 1 step 1 until n do
                          for j ≔ k + 1 step 1 until i do
                             g[i,j] ≔ g[i,j] - g[i,k]×p[j] - p[i]×g[j,k]
                  end
               end k;
            a[n-1] ≔ g[n-1,n-1]; bq[n-1] ≔ g[n,n-1]⭡2;
            a[n] ≔ g[n,n]; t ≔ abs(g[n,n-1]);
            norm ≔ max(norm,absb+abs(a[n-1])+t);
            norm ≔ max(norm,t+abs(a[n]))
         end Housholder tridiagonalization 1;
         integer i,k,m,m1;
         real norm,epsq,lambda,mu,sq1,sq2,u,pq,gamma,t;
         array a[1:n],bq[0:n-1];
         Housholder tridiagonalization 1(n,g,a,bq,norm);
         epsq ≔ 2·25÷-22×norm⭡2;
         comment The tollerance used in the QR iteration depends
           on the square of the relative machine precision. Here 2.25%-22
           is used which is appropriate for a machine with a 36-bit
           mantissa;
         mu ≔ 0; m ≔ n;
         inspect: if m = 0
           then goto return else i ≔ k ≔ m1 ≔ m - 1;
         bq[0] ≔ 0;
         if bq[k] < epsq then
         begin g[m,m] ≔ a[m]; mu ≔ 0; m ≔ k;
            goto inspect
         end;
         for i ≔ i - 1 while bq[i] > epsq do k ≔ i;
         if k = m1 then
         begin comment treat 2 * 2 block separately;
            mu ≔ a[m1]×a[m] - bq[m1]; sq1 ≔ a[m1] + a[m];
            sq2 ≔ sqrt((a[m1]-a[m])⭡2+4×bq[m1]);
            lambda ≔ 0·5×(if sq1>0 then sq1+sq2 else sq1-sq2);
            g[m1,m1] ≔ lambda; g[m,m] ≔ mu / lambda;
            mu ≔ 0; m ≔ m - 2; goto inspect
         end;
         lambda ≔ if abs(a[m]-mu) < 0·5×abs(a[m])
           then a[m] + 0·5×sqrt(bq[m1]) else 0·0;
         mu ≔ a[m]; sq1 ≔ sq2 ≔ u ≔ 0;
         for i ≔ k step 1 until m1 do
            begin comment shortcut single QR iteration;
               gamma ≔ a[i] - lambda - u;
               pq ≔ if sq1 ≠ 1
                 then gamma⭡2/(1-sq1) else (1-sq2)×bq[i-1];
               t ≔ pq + bq[i]; bq[i-1] ≔ sq1 × t; sq2 ≔ sq1;
               sq1 ≔ bq[i] / t; u ≔ sq1 × (gamma+a[i+1]-lambda);
               a[i] ≔ gamma + u + lambda
            end i;
         gamma ≔ a[m] - lambda - u;
         bq[m1] ≔ sq1 ×
           (if sq1≠1 then gamma⭡2/(1-sq1) else (1-sq2)×bq[m1]);
         a[m] ≔ gamma + lambda; goto inspect;
         return:
      end symmetric QR 1;
      for i ≔ 1 step 1 until n do
           for j ≔ 1 step 1 until i do
              a[i,j] ≔ read;
      symmetric QR 1(n,a);
      for i ≔ 1 step 1 until n do
         begin NLCR; ABSFIXT(2,0,i); PRINT(a[i,i]) end
   end
end