// link: eigen real complex vector c_vector i_vector poly
#include <fstream.h>
#include "real.h"
#include "complex.h"
#include "vector.h"
#include "c_vector.h"
#include "i_vector.h"
#include "poly.h"

/*   This program finds eigenvalues and, if desired, eigenvector components
     of a real square matrix of size n by n.  First the matrix is transformed
     so that it consists of one or more companion matrices in diagonal
     position. Then each companion matrix is transformed so that the
     eigenvalues are along the diagonal, with each set of k multiple
     eigenvalues having k-1 ones appearing just above them.
     If an accuracy check is passed, then eigenvalues and eigenvectors are 
     displayed, otherwise precision is increased and the process repeated.
     When eigenvectors are desired, the transforming matrix Q
     is kept track of, and its final form gives the eigenvector
     components.
*/
void diagonal_test();
void diagonalize_companion();
void find_companion();
string solution_string();
string Y_col_string(int);
string Y_col_string_real(int);
string Y_col_string_imag(int);
     
enum eigen_type{group_member, group_end, list_end, eigenvector,
     near_eigenvector};

int  d, distinct_count, eigenvectors_wanted, n, n1, num_blocks,
     test_trouble = FALSE, transform_trouble = FALSE;
real delta;
matrix A, A0; cmatrix B, I, Q, Y;
smatrix Ae;
cvector eigenvalue;
ivector companion_bounds, eigen_mult, v_type; imatrix row, kind;

int main()
{
string msg1 = "\n\
Program to find the eigenvalues and eigenvectors of a real square matrix A\
\n\n";
cout << msg1;

set_precision(30);
int_entry("Enter the number of rows of A  ", n, 1, 20, "", "20 rows is max\n");
n1 = n+1;
A = matrix(n,n); Q = cmatrix(n,n);
int i, j;
I = cmatrix(n,n); for (i=0; i<n; i++) I(i,i) = one;

string msg2 = "\
Enter 0 for eigenvalues only\n\
  and 1 for both eigenvalues and eigenvector components   ";
int_entry(!msg2, eigenvectors_wanted, 0, 1);
  
msg2 = "\
Enter the number of decimal places desired for displayed results   ";
int_entry(!msg2, d, 0, 20, "", "20 decimal places is the maximum\n");
set_precision(d + n*2 + 5);
delta = ten ^ real(-(d+1)); // test value for near eigenvectors
real::list_tokens();
A.entry("enter a", Ae);

if (string::input_from_file == FALSE) {
     ofstream input_log("eigen.log");
     input_log << "// eigen.log\n";
     input_log << "// number of matrix rows\n" << n << "\n";
     input_log << "// yes or no for eigenvectors (1 or 0)\n";
     input_log << eigenvectors_wanted << "\n";
     input_log << "// number of decimal places\n" << d << "\n";
     for (i=0; i<n; i++) {
          for (j=0; j<n; j++) {
               input_log << "// a(" << i+1 << ", " << j+1 << ")\n";
               input_log << Ae(i,j) << "\n";}}
     input_log.close();}

cout << "Now computing eigenvalue ";
if (eigenvectors_wanted) cout << "and eigenvector ";
cout << "approximations\n"; cout.flush();
while (1) {
     if (eigenvectors_wanted) {
          Q = I.copy(); A0 = A.copy();}
     find_companion();
     diagonalize_companion();
     if (transform_trouble) {
          add_precision(12); transform_trouble = FALSE;}
     else {
          diagonal_test();
          if (test_trouble) {
               add_precision(12); test_trouble = FALSE;}
          else {
               if (test_failure) {
                    if (test_failure < 12) test_failure = 12;
                    add_precision();}
               else break;}}
     cout << "Precision increase to " << current_precision() << endl;
     A.evaluate(Ae);}
string disp = solution_string(); 
cout << disp;

ofstream record("eigen.prt");
int w;
record << msg1;
if (eigenvectors_wanted == FALSE) record << "  Only eigenvalues computed\n";
record << "The number of matrix rows = " << n << "\n";
if (n < 10) w = 1; else w = 2;
string::left_adjustment();
for (i=0; i<n; i++) {
     for (j=0; j<n; j++) {
          record << "a(" << ltos(i+1).setw(w) << "," << ltos(j+1).setw(w);
          record << ") = " << Ae(i,j) << "\n";}}
record << "Solution decimal places = " << d << " fixed-point\n";
record << disp;
record.close();

return 0;
}

void diagonal_test()
{
int h, i, j, k; complex z, z1;

// generate discrepancy matrix
matrix E(n,n);
for (h=0; h<distinct_count; h++) {
     k = 0; i = row(h, k); z = B(i,i);
     while (1) {
          for (j=0; j<i; j++) {
               z1 = B(i,j); E(i,j) = abs(z1.x) + abs(z1.y);}
          if (k != 0) {
               z1 = B(i,j) - z; E(i,j) = abs(z1.x) + abs(z1.y);}
          j++;
          if (kind(h,k) != list_end) {
               z1 = B(i,j) - one; E(i,j++) = abs(z1.x) + abs(z1.y);}
          for (; j<n; j++) {
               z1 = B(i,j); E(i,j) = abs(z1.x) + abs(z1.y);}
          if (kind(h,k) == list_end) break;
          i = row(h,++k);}}

// find max row sum of E
real r_max, r, r1;
for (i=0; i<n; i++) {
     r = zero;
     for (j=0; j<n; j++) r += E(i,j);
     r_max = maxi(r_max, r);}

// compute correctly ranged eigenvalues and max eigenvalue multiplicity
eigenvalue = cvector(distinct_count); int max_eigen_mult = 0;
for (k=0; k<distinct_count; k++) {
     h = eigen_mult(k);
     if (h > max_eigen_mult) max_eigen_mult = h;
     if (h == 1) r = r_max;
     else if (h == 2) r = sqrt(r_max);
     else r = r_max ^ (one / real(h));
     i = row(k,0);
     eigenvalue(k) = B(i,i) % r;}

// Check for overlapping eigenvalues; sets may have to be combined.
int k1, repeat;
do {
     repeat = FALSE;
     for (k=0; k<distinct_count; k++) {
          z = eigenvalue(k);
          for (k1=k+1; k1<distinct_count; k1++) {
               z1 = eigenvalue(k1);
               if (z != z1) continue;
               repeat = TRUE;
               v_type(k) = near_eigenvector;
               // add k1 list to k list
               eigenvalue(k) = z % (z1.wid() * two);
               eigen_mult(k) += eigen_mult(k1);
               i = 0;
               while (kind(k,i) != list_end) i++;
               kind(k,i) = group_end; j = 0;
               do {
                    row(k,++i) = row(k1,j);
                    kind(k,i) = kind(k1,j);}
               while (kind(k1,j++) != list_end);
               // move lists after k1 back
               for (i=k1+1; i<distinct_count; i++) {
                    j = i-1;
                    eigenvalue(j) = eigenvalue(i);
                    eigen_mult(j) = eigen_mult(i);
                    v_type(j) = v_type(i);
                    for (h=0; h<n; h++) {
                         row(j,h) = row(i,h);
                         kind(j,h) = kind(i,h);}}
               k1--; distinct_count--;}}}
while (repeat);

// test eigenvalues for accuracy
for (k=0; k<distinct_count; k++) eigenvalue(k).test(1,d);
if (test_failure) {
     if (max_eigen_mult > 1) {
          if (test_failure > 50) test_failure = 50;
          test_failure *= max_eigen_mult;}
     return;}
if (eigenvectors_wanted == FALSE) return;

// construct eigenvectors and near eigenvectors
int i1, i2, j1, l, last_group;
Y = cmatrix(n,n); vector X(n); cvector Y1(n);
for (k=0; k<distinct_count; k++) {
     if (v_type(k) == eigenvector) {
          j1 = row(k,0);
          for (i=0; i<n; i++) Y(i,j1) = Q(i,j1);
          z = eigenvalue(k); z1 = z.wid(); r = z1.x + z1.y + r_max;
          for (i=0; kind(k,i) != list_end; i++) {
               X(row(k,i+1)) = r;}
          for (k1=0; k1<distinct_count; k1++) {
               if (k1 == k) continue;
               i = 0;
               while (1) {
                    j = row(k1,i); r = abs(z - B(j,j)); i1 = i;
                    while (kind(k1,i) == group_member) i++;
                    if (kind(k1,i) == list_end) last_group = TRUE;
                    else last_group = FALSE;
                    j = row(k1,i); X(j) = r_max / r; i2 = i;
                    for (; i>i1; i--) {
                         l = row(k1,i-1); X(l) = (X(j) + r_max) / r; j = l;}
                    if (last_group) break;
                    i = ++i2;}}
          for (i=0; i<n; i++) {
               r = zero;
               for (j=0; j<n; j++) {
                    if (j == j1) continue;
                    z1 = Q(i,j) * X(j); r += abs(z1.x) + abs(z1.y);}
               Y(i,j1) = Y(i,j1) % r;}
          r = - one;
          for (i=0; i<n; i++) {
               r1 = abs_square(Y(i,j1));
               if (r1 > r) {
                    r = r1; j = i;}} 
          z1 = Y(j,j1); Y(j,j1) = one;
          if (z1 == zero) {
               test_trouble = TRUE; return;}
          for (i=0; i<n; i++) if (i != j) Y(i,j1) /= z1;
          for (i=0; i<n; i++) Y(i,j1).test(1,d);}
     else { // near_eigenvector case
          l = 0; z = eigenvalue(k);
          while (1) {
               j1 = row(k,l);
               for (i=0; i<n; i++) Y(i,j1) = Q(i,j1);
               r = -one;
               for (i=0; i<n; i++) {
                    r1 = abs_square(Y(i,j1)); 
                    if (r1 > r) {
                         r = r1; j = i;}} 
               z1 = Y(j,j1); Y(j,j1) = one;
               if (z1 == zero) {
                    test_trouble = TRUE; return;}
               for (i=0; i<n; i++) if (i != j) Y(i,j1) /= z1;
               for (i=0; i<n; i++) {
                    z1 = zero;
                    for (j=0; j<n; j++) z1 += Y(j,j1) * A0(i,j);
                    Y1(i) = z1 - z * Y(i,j1);}
               r = zero;
               for (i=0; i<n; i++) r += abs_square(Y1(i));
               if (sqrt(r) >= delta) {
                    test_trouble = TRUE; return;}
               for (i=0; i<n; i++) Y(i,j1).test(1,d);
               while (kind(k,l) == group_member) l++;
               if (kind(k,l) == list_end) break; else l++;}}}
}
               
void diagonalize_companion()
{
int deg, deg1, h, h_hi, h_lo, h_max, i, j, j_row, k, m, num_roots;
complex multiplier, term, z; real sum, temp, temp1;
vector p, q(n); cvector roots; ivector multiplicity;
row = imatrix(n,n); kind = imatrix(n,n); eigen_mult = ivector(n);
v_type = ivector(n);
cmatrix V(n,n), W;
distinct_count = 0;
for (h=0; h<num_blocks; h++) {
     h_lo = companion_bounds(h); h_hi = companion_bounds(h+1);
     h_max = h_hi - 1; j_row = h_lo - 1;
     deg = h_hi - h_lo; p = vector(deg); deg1 = deg - 1;
     for (i=h_lo; i<h_hi; i++) p(i-h_lo) = -A(i,h_max);
     i = find_roots(p, roots, multiplicity);
     if (i == TRUE) { // find_roots returns TRUE if roots not found
          transform_trouble = TRUE; return;}
     num_roots = roots.num_comps();
     for (i=0; i<num_roots; i++) {
          m = multiplicity(i); eigen_mult(distinct_count) = m; j_row += m;
          v_type(distinct_count) = eigenvector; z = roots(i);
          if (m > 1) for (k=1; k<deg1; k++) q(k) = one; 
          multiplier = one;
          for (k=0; k<deg1; k++) {
               V(j_row, h_lo+k) = multiplier; multiplier *= z;}
          V(j_row--, h_lo+deg1) = multiplier;
          for (j=1; j<m; j++) {
               multiplier = one; temp = one; term = one;
               for (k=j; k<deg1; k++) {
                    temp1 = temp + q(k); q(k) = temp; temp = temp1;
                    V(j_row, h_lo+k) = term;
                    multiplier *= z; term = multiplier * temp;}
               V(j_row--, h_lo+deg1) = term;}
           for (j=0; j<m; j++) {
                row(distinct_count, j) = ++j_row;
               kind(distinct_count, j) = group_member;}
          kind(distinct_count++, j-1) = list_end;}}
B = cmatrix(n,n); // Further work with A may be in the complex field
for (i=0; i<n; i++) for (j=0; j<n; j++) B(i,j).x = A(i,j);
B = V * B;
cmatrix::failure_report_wanted = TRUE;
W = V.inverse();
if (cmatrix::failure) {
     transform_trouble = TRUE; return;}
B = B * W;
if (eigenvectors_wanted) Q = Q * W;
}

void find_companion()
{
int block_start = 0, block_end = 0, h, i, j, k, be1, s, t;
real divisor, multiplier, sum, temp;
companion_bounds = ivector(n1);
companion_bounds(0) = 0; num_blocks = 0;
while (1) { // create companion block
     while (1) { // attempt to add column to block
          be1 = block_end + 1;
          if (be1 == n) {
               block_end++; break;}
          for (i=be1; i<n; i++) {
               if (A(i, block_end) != zero) break;}
          if (i < n) { // non-zero pivot element found
               if (i > be1) { // transform to move pivot element
                    for (k=0; k<n; k++) swap(A(be1,k), A(i,k));// row exchange
                    for (k=0; k<n; k++) swap(A(k,be1), A(k,i));// col exchange
                    if (eigenvectors_wanted) {
                         for (k=0; k<n; k++) swap(Q(k,be1), Q(k,i));}}
               // transform to get another 0 - 1 col
               for (h=0; h<n; h++) { // post-multiply A
                    sum = zero;
                    for (k=0; k<n; k++) {
                         sum += A(h,k) * A(k,block_end);}
                    A(h,be1) = sum;}
               if (eigenvectors_wanted) {
                    for (h=0; h<n; h++) { // post-multiply Q
                         sum = zero;
                         for (k=0; k<n; k++) {
                              sum += Q(h,k).x * A(k,block_end);}
                         Q(h,be1).x = sum;}}
               divisor = A(be1, block_end); A(be1, block_end) = one;
               for (k=0; k<n; k++) { // pre-multiply A
                    if (k != be1) { 
                         multiplier = A(k, block_end) / divisor;
                         A(k, block_end) = zero;
                         for (j=0; j<n; j++) {
                              if (j != block_end) {
                                   A(k,j) -= multiplier * A(be1,j);}}}}
               for (k=0; k<n; k++) {
                    if (k != block_end) {
                         A(be1,k) /= divisor;}}
               block_end++;}
          else { // no non-zero pivot; try to make matrix 0 right of block
               for (k=block_end; k>block_start; k--) {
                    s = k-1;
                    for (j=0; j<n; j++) { // pre-multiply A
                         sum = zero;
                         for (t=be1; t<n; t++) sum += A(k,t) * A(t,j);
                         A(s,j) += sum;}
                    if (eigenvectors_wanted) {
                         for (j=be1; j<n; j++) { // post-multiply Q
                              multiplier = A(k,j);
                              for (t=0; t<n; t++) {
                                   Q(t,j) -= Q(t,s) * multiplier;}}}
                    for (j=be1; j<n; j++) { // post-multiply A
                         multiplier = A(k,j); A(k,j) = zero;
                         for (t=0; t<n; t++) {
                              if (t != k) {
                                   A(t,j) -= multiplier * A(t,s);}}}}
               for (k=be1; k<n; k++) {
                    if (A(block_start,k) != zero) break;}
               if (k < n) { // move element to block_start col
                    for (i=0; i<n; i++) { // col cycle A
                         temp = A(i,k);
                         for (j=k; j>block_start; j--) A(i,j) = A(i,j-1);
                         A(i,block_start) = temp;}
                    if (eigenvectors_wanted) {
                         for (i=0; i<n; i++) { // col cycle Q
                              temp = Q(i,k).x;
                              for (j=k; j>block_start; j--) {
                                   Q(i,j).x = Q(i,j-1).x;}
                              Q(i,block_start).x = temp;}}
                    for (i=0; i<n; i++) { // row cycle A
                         temp = A(k,i);
                         for (j=k; j>block_start; j--) A(j,i) = A(j-1,i);
                         A(block_start,i) = temp;}
                    block_end = block_start;}
               else { // companion block constructed
                    block_end++; break;}}} 
     companion_bounds(++num_blocks) = block_end;
     if (block_end == n) break; else block_start = block_end;}
}

string solution_string()
{
// Sort eigenvalues by real part or, if real parts equal, by imaginary part
int i, j, k, repeat, t, w;
do {
     repeat = FALSE;
     for (i=1; i<distinct_count; i++) {
          j = i - 1;
          if (eigenvalue(j).x < eigenvalue(i).x ||
          eigenvalue(j).x == eigenvalue(i).x &&
          eigenvalue(j).y < eigenvalue(i).y) {
               swap(eigenvalue(i), eigenvalue(j));
               t = eigen_mult(j); eigen_mult(j) = eigen_mult(i);
               eigen_mult(i) = t;
               for (k=0; k<n; k++) {
                    t = row(i,k); row(i,k) = row(j,k);
                    row(j,k) = t;
                    t = kind(i,k); kind(i,k) = kind(j,k); kind(j,k) = t;}
               repeat = TRUE;}}}
while (repeat);

int lengthr = 0, lengthi = 0, fieldr = 14, fieldi = 14;
for (i=0; i<distinct_count; i++) {
     j = eigenvalue(i).x.str(1,d).len(); if (j > lengthr) lengthr = j;
     j = eigenvalue(i).y.str(1,d).len(); if (j > lengthi) lengthi = j;
}

if (lengthr > fieldr) fieldr = lengthr;
if (lengthi > fieldi) fieldi = lengthi;

int total_length = 2 + 3 + fieldr + 2 + fieldi + 2; string out, null;
out = string("\n") + null.setw(total_length/2 - 5) 
      + string("Distinct Eigenvalues\n");

int frontsp_r, backsp_r, frontsp_i, backsp_i;
frontsp_r = (fieldr - 9)/2; backsp_r = fieldr - frontsp_r - 9;
frontsp_i = (fieldi - 14)/2; backsp_i = fieldi - frontsp_i - 14;
out += null.setw(fieldr + fieldi + 11) + string("Apparent\n");
out += string("No.  ") + null.setw(frontsp_r) + string("Real Part");
out += null.setw(backsp_r + 2 + frontsp_i) + string("Imaginary Part");
out += null.setw(backsp_i) + string("  multiplicity\n");

if (lengthr < fieldr) {
     frontsp_r = (fieldr - lengthr)/2;
     backsp_r = fieldr - frontsp_r - lengthr;}
else {frontsp_r = 0; backsp_r = 0;}
if (lengthi < fieldi) {
     frontsp_i = (fieldi - lengthi)/2;
     backsp_i = fieldi - frontsp_i - lengthi;}
else {frontsp_i = 0; backsp_i = 0;}

w = 2;
string::left_adjustment();
for (i=0; i<distinct_count; i++) {
     out += ltos(i+1).setw(w) + string("   ") + null.setw(frontsp_r);
     out += eigenvalue(i).x.str(1, d).setw(lengthr);
     out += null.setw(backsp_r + 2 + frontsp_i);
     out += eigenvalue(i).y.str(1, d).setw(lengthi);
     out += null.setw(backsp_i + 7);
     out += ltos(eigen_mult(i)).setw(w) + string("\n");
}

if (eigenvectors_wanted) {
     int first_near_vector = TRUE, num_vectors;
     out += string("\n                        Eigenvectors\n");
     out += string("All vectors displayed have a 1 component of largest ");
     out += string("magnitude\n");
     for (i=0; i<distinct_count; i++) {
          k = 0; num_vectors = 1;
          while (kind(i,k) != list_end) {
               if (kind(i,k++) == group_end) num_vectors++;}
          out += string("\nEigenvalue No. ") + ltos(i+1);
          if (eigen_mult(i) == 1) out += string(" eigenvector");
          else {
               if (num_vectors == 1) out += string(": 1 eigenvector");
               else {
                    out += string(": ") + ltos(num_vectors);
                    out += string(" near-eigenvectors found");
                    if (first_near_vector) {
                         out += string("\n Note: A near-eigenvector is a ");
                         out += string("vector X such that the length\n");
                         out += string("of (AX - eigenvalue * X) is less ");
                         out += string("than 10^-(d+1)\n where d is the "); 
                         out += string("number of decimal places requested, ");
                         out += ltos(d) + string(" on this run.");
                         first_near_vector = FALSE;}}}
          out += string("\n");
          k = 0;
          do {
               out += Y_col_string(row(i,k));
               while (kind(i,k) == group_member) k++;}
          while(kind(i,k++) != list_end);
     }
}

// drf5n 6 June 1998 additions to output matricies of eigen...

out += string("\n") + null.setw(total_length/2 - 5) 
      + string("Eigenvalue Matrix (Complex)\n");

if (lengthr < fieldr) {
     frontsp_r = (fieldr - lengthr)/2;
     backsp_r = fieldr - frontsp_r - lengthr;}
else {frontsp_r = 0; backsp_r = 0;}
if (lengthi < fieldi) {
     frontsp_i = (fieldi - lengthi)/2;
     backsp_i = fieldi - frontsp_i - lengthi;}
else {frontsp_i = 0; backsp_i = 0;}

w = 2;
string::left_adjustment();
for (i=0; i<distinct_count; i++) {
  for (j=0; j<eigen_mult(i); j++) {
     out += null.setw(frontsp_r);
     out += eigenvalue(i).x.str(1, d).setw(lengthr);
     out += null.setw(backsp_r + 2 + frontsp_i);
     out += eigenvalue(i).y.str(1, d).setw(lengthi);
     out += null.setw(backsp_i + 7);
     out += string("\n");
  }
}

if (eigenvectors_wanted) {
     int first_near_vector = TRUE, num_vectors;
     out += string("\nEigenvector Matrix (Real) (Transposed)\n");
     for (i=0; i<distinct_count; i++) {
          k = 0; num_vectors = 1;
          while (kind(i,k) != list_end) {
               if (kind(i,k++) == group_end) num_vectors++;}
          k = 0;
          do {
               out += Y_col_string_real(row(i,k));
               while (kind(i,k) == group_member) k++;
               out += string("\n");
          }
          while(kind(i,k++) != list_end);
     }
     out += string("\nEigenvector Matrix (Imaginary) (Transposed)\n");
     for (i=0; i<distinct_count; i++) {
          k = 0; num_vectors = 1;
          while (kind(i,k) != list_end) {
               if (kind(i,k++) == group_end) num_vectors++;}
          k = 0;
          do {
               out += Y_col_string_imag(row(i,k));
               while (kind(i,k) == group_member) k++;
               out += string("\n");
          }
          while(kind(i,k++) != list_end);
     }
}

return out;
}

string Y_col_string_real(int col)
{
    int backspace, fieldwidth = 14, frontspace, i, printwidth;
    string zero_field = zero.str(1, d);

    printwidth = zero_field.len() + 1; // 1 space for possible negative sign
    if (printwidth > fieldwidth) fieldwidth = printwidth;
    frontspace = (fieldwidth - 9) / 2; 
    backspace = fieldwidth - 9 - frontspace;
    string Y_col, null;
    Y_col = null.setw(3+frontspace) ;
    frontspace = (fieldwidth - printwidth) / 2;
    backspace = fieldwidth - frontspace - printwidth;
    for (i=0; i<n; i++) {
         Y_col += null.setw(3+frontspace);
         Y_col += Y(i,col).x.str(1, d).setw(printwidth);
         Y_col += null.setw(backspace+3+frontspace);
    }
    return Y_col;
}

string Y_col_string_imag(int col)
{
    int backspace, fieldwidth = 14, frontspace, i, printwidth;
    string zero_field = zero.str(1, d);

    printwidth = zero_field.len() + 1; // 1 space for possible negative sign
    if (printwidth > fieldwidth) fieldwidth = printwidth;
    frontspace = (fieldwidth - 9) / 2; 
    backspace = fieldwidth - 9 - frontspace;
    string Y_col, null;
    Y_col = null.setw(3+frontspace) ;
    frontspace = (fieldwidth - printwidth) / 2;
    backspace = fieldwidth - frontspace - printwidth;
    for (i=0; i<n; i++) {
         Y_col += null.setw(3+frontspace);
         Y_col += Y(i,col).y.str(1, d).setw(printwidth);
         Y_col += null.setw(backspace+3+frontspace);
    }
    return Y_col;
}

string Y_col_string(int col)
{
int backspace, fieldwidth = 14, frontspace, i, printwidth;
string zero_field = zero.str(1, d);
printwidth = zero_field.len() + 1; // 1 space for possible negative sign
if (printwidth > fieldwidth) fieldwidth = printwidth;
frontspace = (fieldwidth - 9) / 2; backspace = fieldwidth - 9 - frontspace;
string Y_col, null;
Y_col = string("Component") + null.setw(3+frontspace) + string("Real Part");
Y_col += null.setw(backspace+3+frontspace) + string("Imag Part\n");
frontspace = (fieldwidth - printwidth) / 2;
backspace = fieldwidth - frontspace - printwidth;
for (i=0; i<n; i++) {
     Y_col += string("    ") + ltos(i+1).setw(2) + string("   ");
     Y_col += null.setw(3+frontspace);
     Y_col += Y(i,col).x.str(1, d).setw(printwidth);
     Y_col += null.setw(backspace+3+frontspace);
     Y_col += Y(i,col).y.str(1, d).setw(printwidth);
     Y_col += string("\n");}
return Y_col;
}