00001 #include "solve_poly_root.h"
00002
00003
00004 #define RADIX 2
00005 #define RADIX2 (RADIX*RADIX)
00006
00007 void
00008 balance_companion_matrix (double *m, size_t nc)
00009 {
00010 int not_converged = 1;
00011
00012 double row_norm = 0;
00013 double col_norm = 0;
00014
00015 while (not_converged)
00016 {
00017 size_t i, j;
00018 double g, f, s;
00019
00020 not_converged = 0;
00021
00022 for (i = 0; i < nc; i++)
00023 {
00024
00025
00026 if (i != nc - 1)
00027 {
00028 col_norm = fabs (MAT (m, i + 1, i, nc));
00029 }
00030 else
00031 {
00032 col_norm = 0;
00033
00034 for (j = 0; j < nc - 1; j++)
00035 {
00036 col_norm += fabs (MAT (m, j, nc - 1, nc));
00037 }
00038 }
00039
00040
00041
00042 if (i == 0)
00043 {
00044 row_norm = fabs (MAT (m, 0, nc - 1, nc));
00045 }
00046 else if (i == nc - 1)
00047 {
00048 row_norm = fabs (MAT (m, i, i - 1, nc));
00049 }
00050 else
00051 {
00052 row_norm = (fabs (MAT (m, i, i - 1, nc))
00053 + fabs (MAT (m, i, nc - 1, nc)));
00054 }
00055
00056 if (col_norm == 0 || row_norm == 0)
00057 {
00058 continue;
00059 }
00060
00061 g = row_norm / RADIX;
00062 f = 1;
00063 s = col_norm + row_norm;
00064
00065 while (col_norm < g)
00066 {
00067 f *= RADIX;
00068 col_norm *= RADIX2;
00069 }
00070
00071 g = row_norm * RADIX;
00072
00073 while (col_norm > g)
00074 {
00075 f /= RADIX;
00076 col_norm /= RADIX2;
00077 }
00078
00079 if ((row_norm + col_norm) < 0.95 * s * f)
00080 {
00081 not_converged = 1;
00082
00083 g = 1 / f;
00084
00085 if (i == 0)
00086 {
00087 MAT (m, 0, nc - 1, nc) *= g;
00088 }
00089 else
00090 {
00091 MAT (m, i, i - 1, nc) *= g;
00092 MAT (m, i, nc - 1, nc) *= g;
00093 }
00094
00095 if (i == nc - 1)
00096 {
00097 for (j = 0; j < nc; j++)
00098 {
00099 MAT (m, j, i, nc) *= f;
00100 }
00101 }
00102 else
00103 {
00104 MAT (m, i + 1, i, nc) *= f;
00105 }
00106 }
00107 }
00108 }
00109 }