1 files changed, 127 insertions, 0 deletions
diff --git a/gsl-1.9/poly/balance.c b/gsl-1.9/poly/balance.c
new file mode 100644
index 0000000..a1f7916
--- /dev/null
+++ b/gsl-1.9/poly/balance.c
@@ -0,0 +1,127 @@
+/* poly/balance.c
+ * 
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+static void balance_companion_matrix (double *m, size_t n);
+
+#define RADIX 2
+#define RADIX2 (RADIX*RADIX)
+
+static void
+balance_companion_matrix (double *m, size_t nc)
+{
+  int not_converged = 1;
+
+  double row_norm = 0;
+  double col_norm = 0;
+
+  while (not_converged)
+    {
+      size_t i, j;
+      double g, f, s;
+
+      not_converged = 0;
+
+      for (i = 0; i < nc; i++)
+        {
+          /* column norm, excluding the diagonal */
+
+          if (i != nc - 1)
+            {
+              col_norm = fabs (MAT (m, i + 1, i, nc));
+            }
+          else
+            {
+              col_norm = 0;
+
+              for (j = 0; j < nc - 1; j++)
+                {
+                  col_norm += fabs (MAT (m, j, nc - 1, nc));
+                }
+            }
+
+          /* row norm, excluding the diagonal */
+
+          if (i == 0)
+            {
+              row_norm = fabs (MAT (m, 0, nc - 1, nc));
+            }
+          else if (i == nc - 1)
+            {
+              row_norm = fabs (MAT (m, i, i - 1, nc));
+            }
+          else
+            {
+              row_norm = (fabs (MAT (m, i, i - 1, nc)) 
+                          + fabs (MAT (m, i, nc - 1, nc)));
+            }
+
+          if (col_norm == 0 || row_norm == 0)
+            {
+              continue;
+            }
+
+          g = row_norm / RADIX;
+          f = 1;
+          s = col_norm + row_norm;
+
+          while (col_norm < g)
+            {
+              f *= RADIX;
+              col_norm *= RADIX2;
+            }
+
+          g = row_norm * RADIX;
+
+          while (col_norm > g)
+            {
+              f /= RADIX;
+              col_norm /= RADIX2;
+            }
+
+          if ((row_norm + col_norm) < 0.95 * s * f)
+            {
+              not_converged = 1;
+
+              g = 1 / f;
+
+              if (i == 0)
+                {
+                  MAT (m, 0, nc - 1, nc) *= g;
+                }
+              else
+                {
+                  MAT (m, i, i - 1, nc) *= g;
+                  MAT (m, i, nc - 1, nc) *= g;
+                }
+
+              if (i == nc - 1)
+                {
+                  for (j = 0; j < nc; j++)
+                    {
+                      MAT (m, j, i, nc) *= f;
+                    }
+                }
+              else
+                {
+                  MAT (m, i + 1, i, nc) *= f;
+                }
+            }
+        }
+    }
+}