Copy of aRts for Trinity modifications

git-svn-id: svn://anonsvn.kde.org/home/kde/branches/trinity/dependencies/arts@1070145 283d02a7-25f6-0310-bc7c-ecb5cbfe19da

1 files changed, 1085 insertions, 0 deletions

diff --git a/flow/gsl/gslmath.c b/flow/gsl/gslmath.c
new file mode 100644
index 0000000..97ac675
--- /dev/null
+++ b/flow/gsl/gslmath.c
@@ -0,0 +1,1085 @@
+/* GSL - Generic Sound Layer
+ * Copyright (C) 2001 Stefan Westerfeld and Tim Janik
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2 of the License, or (at your option) any later version.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General
+ * Public License along with this library; if not, write to the
+ * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
+ * Boston, MA 02111-1307, USA.
+ */
+#include "gslmath.h"
+
+#include <string.h>
+#include <stdio.h>
+
+
+#define RING_BUFFER_LENGTH	(16)
+#define	PRINTF_DIGITS		"1270"
+#define	FLOAT_STRING_SIZE	(2048)
+
+
+/* factorization constants: 2^(1/12), ln(2^(1/12)) and 2^(1/(12*6))
+ * retrived with:
+ #include <stl.h>
+ #include <complex.h>
+ typedef long double ld;
+ 
+ int main (void)
+ {
+ ld r, l;
+ 
+ cout.precision(256);
+ 
+ r = pow ((ld) 2, (ld) 1 / (ld) 12);
+ cout << "2^(1/12) =\n";
+ cout << "2^" << (ld) 1 / (ld) 12 << " =\n";
+ cout << r << "\n";
+ 
+ l = log (r);
+ cout << "ln(2^(1/12)) =\n";
+ cout << "ln(" << r << ") =\n";
+ cout << l << "\n";
+ 
+ r = pow ((ld) 2, (ld) 1 / (ld) 72);
+ cout << "2^(1/72) =\n";
+ cout << "2^" << (ld) 1 / (ld) 72 << " =\n";
+ cout << r << "\n";
+ 
+ return 0;
+ }
+*/
+
+
+/* --- prototypes --- */
+static void		zrhqr	(double a[], int m, double rtr[], double rti[]);
+static double		rf	(double x, double y, double z);
+static double		ellf	(double phi, double ak);
+static void		sncndn	(double uu, double emmc, double *sn_p, double *cn_p, double *dn_p);
+static void		sncndnC	(GslComplex uu, GslComplex emmc, GslComplex *sn_p, GslComplex *cn_p, GslComplex *dn_p);
+static GslComplex	rfC	(GslComplex x, GslComplex y, GslComplex z);
+
+
+/* --- functions --- */
+static inline char*
+pretty_print_double (char  *str,
+		     double d)
+{
+  char *s= str;
+  
+  sprintf (s, "%."PRINTF_DIGITS"f", d);
+  while (*s) s++;
+  while (s[-1] == '0' && s[-2] != '.')
+    s--;
+  *s = 0;
+  return s;
+}
+
+char*
+gsl_complex_list (unsigned int n_points,
+		  GslComplex  *points,
+		  const char  *indent)
+{
+  static unsigned int rbi = 0;
+  static char* rbuffer[RING_BUFFER_LENGTH] = { NULL, };
+  char *s, *tbuffer = g_newa (char, (FLOAT_STRING_SIZE * 2 * n_points));
+  unsigned int i;
+  
+  rbi++; if (rbi >= RING_BUFFER_LENGTH) rbi -= RING_BUFFER_LENGTH;
+  if (rbuffer[rbi] != NULL)
+    g_free (rbuffer[rbi]);
+  s = tbuffer;
+  for (i = 0; i < n_points; i++)
+    {
+      *s = 0;
+      if (indent)
+	strcat (s, indent);
+      while (*s) s++;
+      s = pretty_print_double (s, points[i].re);
+      *s++ = ' ';
+      s = pretty_print_double (s, points[i].im);
+      *s++ = '\n';
+    }
+  *s++ = 0;
+  rbuffer[rbi] = g_strdup (tbuffer);
+  return rbuffer[rbi];
+}
+
+char*
+gsl_complex_str (GslComplex c)
+{
+  static unsigned int rbi = 0;
+  static char* rbuffer[RING_BUFFER_LENGTH] = { NULL, };
+  char *s, tbuffer[FLOAT_STRING_SIZE * 2];
+  
+  rbi++; if (rbi >= RING_BUFFER_LENGTH) rbi -= RING_BUFFER_LENGTH;
+  if (rbuffer[rbi] != NULL)
+    g_free (rbuffer[rbi]);
+  s = tbuffer;
+  *s++ = '{';
+  s = pretty_print_double (s, c.re);
+  *s++ = ',';
+  *s++ = ' ';
+  s = pretty_print_double (s, c.im);
+  *s++ = '}';
+  *s++ = 0;
+  rbuffer[rbi] = g_strdup (tbuffer);
+  return rbuffer[rbi];
+}
+
+char*
+gsl_poly_str (unsigned int degree,
+	      double      *a,
+	      const char  *var)
+{
+  static unsigned int rbi = 0;
+  static char* rbuffer[RING_BUFFER_LENGTH] = { NULL, };
+  char *s, *tbuffer = g_newa (char, degree * FLOAT_STRING_SIZE);
+  unsigned int i;
+  
+  if (!var)
+    var = "x";
+  rbi++; if (rbi >= RING_BUFFER_LENGTH) rbi -= RING_BUFFER_LENGTH;
+  if (rbuffer[rbi] != NULL)
+    g_free (rbuffer[rbi]);
+  s = tbuffer;
+  *s++ = '(';
+  s = pretty_print_double (s, a[0]);
+  for (i = 1; i <= degree; i++)
+    {
+      *s++ = '+';
+      *s = 0; strcat (s, var); while (*s) s++;
+      *s++ = '*';
+      *s++ = '(';
+      s = pretty_print_double (s, a[i]);
+    }
+  while (i--)
+    *s++ = ')';
+  *s++ = 0;
+  rbuffer[rbi] = g_strdup (tbuffer);
+  return rbuffer[rbi];
+}
+
+char*
+gsl_poly_str1 (unsigned int degree,
+	       double      *a,
+	       const char  *var)
+{
+  static unsigned int rbi = 0;
+  static char* rbuffer[RING_BUFFER_LENGTH] = { NULL, };
+  char *s, *tbuffer = g_newa (char, degree * FLOAT_STRING_SIZE);
+  unsigned int i, need_plus = 0;
+  
+  if (!var)
+    var = "x";
+  rbi++; if (rbi >= RING_BUFFER_LENGTH) rbi -= RING_BUFFER_LENGTH;
+  if (rbuffer[rbi] != NULL)
+    g_free (rbuffer[rbi]);
+  s = tbuffer;
+  *s++ = '(';
+  if (a[0] != 0.0)
+    {
+      s = pretty_print_double (s, a[0]);
+      need_plus = 1;
+    }
+  for (i = 1; i <= degree; i++)
+    {
+      if (a[i] == 0.0)
+	continue;
+      if (need_plus)
+	{
+	  *s++ = ' ';
+	  *s++ = '+';
+	  *s++ = ' ';
+	}
+      if (a[i] != 1.0)
+	{
+	  s = pretty_print_double (s, a[i]);
+	  *s++ = '*';
+	}
+      *s = 0;
+      strcat (s, var);
+      while (*s) s++;
+      if (i > 1)
+	{
+	  *s++ = '*';
+	  *s++ = '*';
+	  sprintf (s, "%u", i);
+	  while (*s) s++;
+	}
+      need_plus = 1;
+    }
+  *s++ = ')';
+  *s++ = 0;
+  rbuffer[rbi] = g_strdup (tbuffer);
+  return rbuffer[rbi];
+}
+
+void
+gsl_complex_gnuplot (const char  *file_name,
+		     unsigned int n_points,
+		     GslComplex  *points)
+{
+  FILE *fout = fopen (file_name, "w");
+  
+  fputs (gsl_complex_list (n_points, points, ""), fout);
+  fclose (fout);
+}
+
+double
+gsl_temp_freq (double kammer_freq,
+	       int    halftone_delta)
+{
+  double factor;
+  
+  factor = pow (GSL_2_POW_1_DIV_12, halftone_delta);
+  
+  return kammer_freq * factor;
+}
+
+void
+gsl_poly_from_re_roots (unsigned int degree,
+			double      *a,
+			GslComplex  *roots)
+{
+  unsigned int i;
+  
+  /* initialize polynomial */
+  a[1] = 1;
+  a[0] = -roots[0].re;
+  /* monomial factor multiplication */
+  for (i = 1; i < degree; i++)
+    {
+      unsigned int j;
+      
+      a[i + 1] = a[i];
+      for (j = i; j >= 1; j--)
+	a[j] = a[j - 1] - a[j] * roots[i].re;
+      a[0] *= -roots[i].re;
+    }
+}
+
+void
+gsl_cpoly_from_roots (unsigned int degree,
+		      GslComplex  *c,
+		      GslComplex  *roots)
+{
+  unsigned int i;
+  
+  /* initialize polynomial */
+  c[1].re = 1;
+  c[1].im = 0;
+  c[0].re = -roots[0].re;
+  c[0].im = -roots[0].im;
+  /* monomial factor multiplication */
+  for (i = 1; i < degree; i++)
+    {
+      GslComplex r = gsl_complex (-roots[i].re, -roots[i].im);
+      unsigned int j;
+      
+      c[i + 1] = c[i];
+      for (j = i; j >= 1; j--)
+	c[j] = gsl_complex_add (c[j - 1], gsl_complex_mul (c[j], r));
+      c[0] = gsl_complex_mul (c[0], r);
+    }
+}
+
+void
+gsl_poly_complex_roots (unsigned int degree,
+			double      *a,		/* [0..degree] (degree+1 elements) */
+			GslComplex  *roots)	/* [degree] */
+{
+  double *roots_re = g_newa (double, 1 + degree);
+  double *roots_im = g_newa (double, 1 + degree);
+  unsigned int i;
+  
+  zrhqr (a, degree, roots_re, roots_im);
+  for (i = 0; i < degree; i++)
+    {
+      roots[i].re = roots_re[1 + i];
+      roots[i].im = roots_im[1 + i];
+    }
+}
+
+double
+gsl_ellip_rf (double x,
+	      double y,
+	      double z)
+{
+  return rf (x, y, z);
+}
+
+double
+gsl_ellip_F (double phi,
+	     double ak)
+{
+  return ellf (phi, ak);
+}
+
+double
+gsl_ellip_sn (double u,
+	      double emmc)
+{
+  double sn;
+  
+  sncndn (u, emmc, &sn, NULL, NULL);
+  return sn;
+}
+
+double
+gsl_ellip_asn (double y,
+	       double emmc)
+{
+  return y * rf (1.0 - y * y, 1.0 - y * y * (1.0 - emmc), 1.0);
+}
+
+GslComplex
+gsl_complex_ellip_asn (GslComplex y,
+		       GslComplex emmc)
+{
+  return gsl_complex_mul (y,
+			  rfC (gsl_complex_sub (gsl_complex (1.0, 0),
+						gsl_complex_mul (y, y)),
+			       gsl_complex_sub (gsl_complex (1.0, 0),
+						gsl_complex_mul3 (y, y, gsl_complex_sub (gsl_complex (1.0, 0),
+											 emmc))),
+			       gsl_complex (1.0, 0)));
+}
+
+GslComplex
+gsl_complex_ellip_sn (GslComplex u,
+		      GslComplex emmc)
+{
+  GslComplex sn;
+  
+  sncndnC (u, emmc, &sn, NULL, NULL);
+  return sn;
+}
+
+double
+gsl_bit_depth_epsilon (guint n_bits)
+{
+  /* epsilon for various bit depths, based on significance of one bit,
+   * minus fudge. created with:
+   * { echo "scale=40"; for i in `seq 1 32` ; do echo "1/2^$i - 10^-($i+1)" ; done } | bc | sed 's/$/,/'
+   */
+  static const double bit_epsilons[] = {
+    .4900000000000000000000000000000000000000,
+    .2490000000000000000000000000000000000000,
+    .1249000000000000000000000000000000000000,
+    .0624900000000000000000000000000000000000,
+    .0312490000000000000000000000000000000000,
+    .0156249000000000000000000000000000000000,
+    .0078124900000000000000000000000000000000,
+    .0039062490000000000000000000000000000000,
+    .0019531249000000000000000000000000000000,
+    .0009765624900000000000000000000000000000,
+    .0004882812490000000000000000000000000000,
+    .0002441406249000000000000000000000000000,
+    .0001220703124900000000000000000000000000,
+    .0000610351562490000000000000000000000000,
+    .0000305175781249000000000000000000000000,
+    .0000152587890624900000000000000000000000,
+    .0000076293945312490000000000000000000000,
+    .0000038146972656249000000000000000000000,
+    .0000019073486328124900000000000000000000,
+    .0000009536743164062490000000000000000000,
+    .0000004768371582031249000000000000000000,
+    .0000002384185791015624900000000000000000,
+    .0000001192092895507812490000000000000000,
+    .0000000596046447753906249000000000000000,
+    .0000000298023223876953124900000000000000,
+    .0000000149011611938476562490000000000000,
+    .0000000074505805969238281249000000000000,
+    .0000000037252902984619140624900000000000,
+    .0000000018626451492309570312490000000000,
+    .0000000009313225746154785156249000000000,
+    .0000000004656612873077392578124900000000,
+    .0000000002328306436538696289062490000000,
+  };
+
+  return bit_epsilons[CLAMP (n_bits, 1, 32) - 1];
+}
+
+
+/* --- Numerical Receipes --- */
+#define	gsl_complex_rmul(scale, c)	gsl_complex_scale (c, scale)
+#define	ONE				gsl_complex (1.0, 0)
+#define	SIGN(a,b)	((b) >= 0.0 ? fabs (a) : -fabs(a))
+static inline int IMAX (int i1, int i2) { return i1 > i2 ? i1 : i2; }
+static inline double DMIN (double d1, double d2) { return d1 < d2 ? d1 : d2; }
+static inline double DMAX (double d1, double d2) { return d1 > d2 ? d1 : d2; }
+static inline double DSQR (double d) { return d == 0.0 ? 0.0 : d * d; }
+#define nrerror(error)	g_error ("NR-ERROR: %s", (error))
+static inline double* vector (long nl, long nh)
+     /* allocate a vector with subscript range v[nl..nh] */
+{
+  double *v = g_new (double, nh - nl + 1 + 1);
+  return v - nl + 1;
+}
+static inline void free_vector (double *v, long nl, long nh)
+{
+  g_free (v + nl - 1);
+}
+static inline double** matrix (long nrl, long nrh, long ncl, long nch)
+     /* allocate a matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+  long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
+  double **m = g_new (double*, nrow + 1);
+  m += 1;
+  m -= nrl;
+  m[nrl] = g_new (double, nrow * ncol + 1);
+  m[nrl] += 1;
+  m[nrl] -= ncl;
+  for (i = nrl + 1; i <= nrh; i++)
+    m[i] = m[i - 1] + ncol;
+  return m;
+}
+static inline void free_matrix (double **m, long nrl, long nrh, long ncl, long nch)
+{
+  g_free (m[nrl] + ncl - 1);
+  g_free (m + nrl - 1);
+}
+
+static void
+poldiv (double u[], int n, double v[], int nv, double q[], double r[])
+     /* Given the n+1 coefficients of a polynomial of degree n in u[0..n], and the nv+1 coefficients
+	of another polynomial of degree nv in v[0..nv], divide the polynomial u by the polynomial
+	v ("u"/"v") giving a quotient polynomial whose coefficients are returned in q[0..n], and a
+	remainder polynomial whose coefficients are returned in r[0..n]. The elements r[nv..n]
+	and q[n-nv+1..n] are returned as zero. */
+{
+  int k,j;
+  
+  for (j=0;j<=n;j++) {
+    r[j]=u[j];
+    q[j]=0.0;
+  }for (k=n-nv;k>=0;k--) {
+    q[k]=r[nv+k]/v[nv];
+    for (j=nv+k-1;j>=k;j--) r[j] -= q[k]*v[j-k];
+  }for (j=nv;j<=n;j++) r[j]=0.0;
+}
+
+#define	MAX_ITER_BASE	9	/* TIMJ: was 3 */
+#define	MAX_ITER_FAC	20	/* TIMJ: was 10 */
+static void
+hqr (double **a, int n, double wr[], double wi[])
+     /* Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n]. On input a can be
+	exactly as output from elmhes §11.5; on output it is destroyed. The real and imaginary parts
+	of the eigenvalues are returned in wr[1..n] and wi[1..n], respectively. */
+{
+  int nn,m,l,k,j,its,i,mmin;
+  double z,y,x,w,v,u,t,s,r,q,p,anorm;
+  r=q=p=0; /* TIMJ: silence compiler */
+  
+  anorm=0.0;                                  /* Compute matrix norm for possible use in lo- */
+  for (i=1;i<=n;i++)                          /* cating single small subdiagonal element. */
+    for (j=IMAX (i-1,1);j<=n;j++)
+      anorm += fabs (a[i][j]);
+  nn=n;
+  t=0.0;                                      /* Gets changed only by an exceptional shift. */
+  while (nn >= 1) {                           /* Begin search for next eigenvalue. */
+    its=0;
+    do {for (l=nn;l>=2;l--) {                 /* Begin iteration: look for single small subdi- */
+      s=fabs (a[l-1][l-1])+fabs (a[l][l]);      /* agonal element. */
+      if (s == 0.0) s=anorm;
+      if ((double)(fabs (a[l][l-1]) + s) == s) break;
+    }
+    x=a[nn][nn];
+    if (l == nn) {                    /* One root found. */
+      wr[nn]=x+t;
+      wi[nn--]=0.0;
+    } else {
+      y=a[nn-1][nn-1];
+      w=a[nn][nn-1]*a[nn-1][nn];
+      if (l == (nn-1)) {            /* Two roots found... */
+	p=0.5*(y-x);
+	q=p*p+w;
+	z=sqrt (fabs (q));
+	x += t;
+	if (q >= 0.0) {           /* ...a real pair. */
+	  z=p+SIGN (z,p);
+	  wr[nn-1]=wr[nn]=x+z;
+	  if (z) wr[nn]=x-w/z;
+	  wi[nn-1]=wi[nn]=0.0;
+	} else {                  /* ...a complex pair. */
+	  wr[nn-1]=wr[nn]=x+p;
+	  wi[nn-1]= -(wi[nn]=z);
+	}
+	nn -= 2;
+      } else {                      /* No roots found. Continue iteration. */
+	if (its == MAX_ITER_BASE * MAX_ITER_FAC)
+	  nrerror ("Too many iterations in hqr");
+	if (its && !(its%MAX_ITER_FAC)) {                /* Form exceptional shift. */
+	  t += x;
+	  for (i=1;i<=nn;i++) a[i][i] -= x;
+	  s=fabs (a[nn][nn-1])+fabs (a[nn-1][nn-2]);
+	  y=x=0.75*s;
+	  w = -0.4375*s*s;
+	}
+	++its;
+	for (m=(nn-2);m>=l;m--) {                    /* Form shift and then look for */
+	  z=a[m][m];                                 /* 2 consecutive small sub- */
+	  r=x-z;                                     /* diagonal elements. */
+	  s=y-z;
+	  p=(r*s-w)/a[m+1][m]+a[m][m+1];           /* Equation (11.6.23). */
+	  q=a[m+1][m+1]-z-r-s;
+	  r=a[m+2][m+1];
+	  s=fabs (p)+fabs (q)+fabs (r);             /* Scale to prevent overflow or */
+	  p /= s;                                    /* underflow. */
+	  q /= s;
+	  r /= s;
+	  if (m == l) break;
+	  u=fabs (a[m][m-1])*(fabs (q)+fabs (r));
+	  v=fabs (p)*(fabs (a[m-1][m-1])+fabs (z)+fabs (a[m+1][m+1]));
+	  if ((double)(u+v) == v)
+	    break;          /* Equation (11.6.26). */
+	}
+	for (i=m+2;i<=nn;i++) {
+	  a[i][i-2]=0.0;
+	  if (i != (m+2))
+	    a[i][i-3]=0.0;
+	}
+	for (k=m;k<=nn-1;k++) {
+	  /* Double QR step on rows l to nn and columns m to nn. */
+	  if (k != m) {
+	    p=a[k][k-1];                      /* Begin setup of Householder */
+	    q=a[k+1][k-1];                    /* vector. */
+	    r=0.0;
+	    if (k != (nn-1)) r=a[k+2][k-1];
+	    if ((x=fabs (p)+fabs (q)+fabs (r)) != 0.0) {
+	      p /= x;                     /* Scale to prevent overflow or */
+	      q /= x;                     /* underflow. */
+	      r /= x;
+	    }
+	  }
+	  if ((s=SIGN (sqrt (p*p+q*q+r*r),p)) != 0.0) {
+	    if (k == m) {
+	      if (l != m)
+		a[k][k-1] = -a[k][k-1];
+	    } else
+	      a[k][k-1] = -s*x;
+	    p += s;                           /* Equations (11.6.24). */
+	    x=p/s;
+	    y=q/s;
+	    z=r/s;
+	    q /= p;
+	    r /= p;
+	    for (j=k;j<=nn;j++) {             /* Row modification. */
+	      p=a[k][j]+q*a[k+1][j];
+	      if (k != (nn-1)) {
+		p += r*a[k+2][j];
+		a[k+2][j] -= p*z;
+	      }
+	      a[k+1][j] -= p*y;
+	      a[k][j] -= p*x;
+	    }
+	    mmin = nn<k+3 ? nn : k+3;
+	    for (i=l;i<=mmin;i++) {           /* Column modification. */
+	      p=x*a[i][k]+y*a[i][k+1];
+	      if (k != (nn-1)) {
+		p += z*a[i][k+2];
+		a[i][k+2] -= p*r;
+	      }a[i][k+1] -= p*q;
+	      a[i][k] -= p;
+	    }
+	  }
+	}
+      }
+    }
+    } while (l < nn-1);
+  }
+}
+
+#define RADIX 2.0
+static void
+balanc (double **a, int n)
+     /* Given a matrix a[1..n][1..n], this routine replaces it by a balanced matrix with identical
+	eigenvalues. A symmetric matrix is already balanced and is unaffected by this procedure. The
+	parameter RADIX should be the machine's floating-point radix. */
+{
+  int last,j,i;
+  double s,r,g,f,c,sqrdx;
+  
+  sqrdx=RADIX*RADIX;
+  last=0;
+  while (last == 0) {
+    last=1;
+    for (i=1;i<=n;i++) {        /* Calculate row and column norms. */
+      r=c=0.0;
+      for (j=1;j<=n;j++)
+	if (j != i) {
+	  c += fabs (a[j][i]);
+	  r += fabs (a[i][j]);
+	}
+      if (c && r) {     /*  If both are nonzero, */
+	g=r/RADIX;
+	f=1.0;
+	s=c+r;
+	while (c<g) {   /* find the integer power of the machine radix that */
+	  f *= RADIX;   /*  comes closest to balancing the matrix. */
+	  c *= sqrdx;
+	}
+	g=r*RADIX;
+	while (c>g) {
+	  f /= RADIX;
+	  c /= sqrdx;
+	}
+	if ((c+r)/f < 0.95*s) {
+	  last=0;
+	  g=1.0/f;
+	  for (j=1;j<=n;j++)
+	    a[i][j] *= g;  /*  Apply similarity transformation */
+	  for (j=1;j<=n;j++) a[j][i] *= f;
+	}
+      }
+    }
+  }
+}
+
+#define MAX_DEGREE 50
+
+static void
+zrhqr (double a[], int m, double rtr[], double rti[])
+     /* Find all the roots of a polynomial with real coefficients, E(i=0..m) a(i)x^i, given the degree m
+	and the coefficients a[0..m]. The method is to construct an upper Hessenberg matrix whose
+	eigenvalues are the desired roots, and then use the routines balanc and hqr. The real and
+	imaginary parts of the roots are returned in rtr[1..m] and rti[1..m], respectively. */
+{
+  int j,k;
+  double **hess,xr,xi;
+  
+  hess=matrix (1,MAX_DEGREE,1,MAX_DEGREE);
+  if (m > MAX_DEGREE || a[m] == 0.0 || /* TIMJ: */ fabs (a[m]) < 1e-15 )
+    nrerror ("bad args in zrhqr");
+  for (k=1;k<=m;k++)	            /* Construct the matrix. */
+    {
+      hess[1][k] = -a[m-k]/a[m];
+      for (j=2;j<=m;j++)
+	hess[j][k]=0.0;
+      if (k != m)
+	hess[k+1][k]=1.0;
+    }
+  balanc (hess,m);                  /* Find its eigenvalues. */
+  hqr (hess,m,rtr,rti);
+  if (0)		/* TIMJ: don't need sorting */
+    for (j=2;j<=m;j++)
+      {				    /* Sort roots by their real parts by straight insertion. */
+	xr=rtr[j];
+	xi=rti[j];
+	for (k=j-1;k>=1;k--)
+	  {
+	    if (rtr[k] <= xr)
+	      break;
+	    rtr[k+1]=rtr[k];
+	    rti[k+1]=rti[k];
+	  }
+	rtr[k+1]=xr;
+	rti[k+1]=xi;
+      }
+  free_matrix (hess,1,MAX_DEGREE,1,MAX_DEGREE);
+}
+
+
+#define EPSS 2.0e-16	/* TIMJ, was(float): 1.0e-7 */
+#define MR 8
+#define MT 100		/* TIMJ: was: 10 */
+#define MAXIT (MT*MR)
+/* Here EPSS is the estimated fractional roundoff error. We try to break (rare) limit cycles with
+   MR different fractional values, once every MT steps, for MAXIT total allowed iterations. */
+
+static void
+laguer (GslComplex a[], int m, GslComplex *x, int *its)
+     /* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial         mi=0 a[i]xi,
+	and given a complex value x, this routine improves x by Laguerre's method until it converges,
+	within the achievable roundoff limit, to a root of the given polynomial. The number of iterations
+	taken is returned as its. */
+{
+  int iter,j;
+  double abx,abp,abm,err;
+  GslComplex dx,x1,b,d,f,g,h,sq,gp,gm,g2;
+  static double frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
+  /* Fractions used to break a limit cycle. */
+  
+  for (iter=1;iter<=MAXIT;iter++)
+    {        /* Loop over iterations up to allowed maximum. */
+      *its=iter;
+      b=a[m];
+      err=gsl_complex_abs (b);
+      d=f=gsl_complex (0.0,0.0);
+      abx=gsl_complex_abs (*x);
+      for (j=m-1;j>=0;j--)
+	{            /* Efficient computation of the polynomial and */
+	  f=gsl_complex_add (gsl_complex_mul (*x,f),d);         /* its first two derivatives. */
+	  d=gsl_complex_add (gsl_complex_mul (*x,d),b);
+	  b=gsl_complex_add (gsl_complex_mul (*x,b),a[j]);
+	  err=gsl_complex_abs (b)+abx*err;
+	}
+      err *= EPSS;
+      /* Estimate of roundoff error in evaluating polynomial. */
+      if (gsl_complex_abs (b) <= err)
+	return;               /* We are on the root. */
+      g=gsl_complex_div (d,b);                              /* The generic case: use Laguerre's formula. */
+      g2=gsl_complex_mul (g,g);
+      h=gsl_complex_sub (g2,gsl_complex_rmul (2.0,gsl_complex_div (f,b)));
+      sq=gsl_complex_sqrt (gsl_complex_rmul ((double) (m-1),gsl_complex_sub (gsl_complex_rmul ((double) m,h),g2)));
+      gp=gsl_complex_add (g,sq);
+      gm=gsl_complex_sub (g,sq);
+      abp=gsl_complex_abs (gp);
+      abm=gsl_complex_abs (gm);
+      if (abp < abm)
+	gp=gm;
+      dx=((DMAX (abp,abm) > 0.0 ? gsl_complex_div (gsl_complex ((double) m,0.0),gp)
+	   : gsl_complex_rmul (1+abx,gsl_complex (cos ((double)iter),sin ((double)iter)))));
+      x1=gsl_complex_sub (*x,dx);
+      if (x->re == x1.re && x->im == x1.im)
+	return;                 /* Converged. */
+      if (iter % MT) *x=x1;
+      else *x=gsl_complex_sub (*x,gsl_complex_rmul (frac[iter/MT],dx));
+      /* Every so often we take a fractional step, to break any limit cycle (itself a rare occurrence). */
+    }
+  nrerror ("too many iterations in laguer");
+  /* Very unusual - can occur only for complex roots. Try a different starting guess for the root. */
+}
+
+/* Here is a driver routine that calls laguer in succession for each root, performs
+   the deflation, optionally polishes the roots by the same Laguerre method - if you
+   are not going to polish in some other way - and finally sorts the roots by their real
+   parts. (We will use this routine in Chapter 13.) */
+
+#define EPS	4.0e-15		/* TIMJ, was(float): 2.0e-6 */
+#define MAXM 100
+/* A small number, and maximum anticipated value of m. */
+
+static void
+zroots (GslComplex a[], int m, GslComplex roots[], int polish)
+     /* Given the degree m and the m+1 complex coefficients a[0..m] of the polynomial mi=0 a (i)xi,
+	this routine successively calls laguer and finds all m complex roots in roots[1..m]. The
+	boolean variable polish should be input as true (1) if polishing (also by Laguerre's method)
+	is desired, false (0) if the roots will be subsequently polished by other means. */
+{
+  int i,its,j,jj;
+  GslComplex x,b,c,ad[MAXM];
+  
+  for (j=0;j<=m;j++) ad[j]=a[j];              /* Copy of coefficients for successive deflation. */
+  for (j=m;j>=1;j--)                          /* Loop over each root to be found. */
+    {
+      x=gsl_complex (0.0,0.0);                  /* Start at zero to favor convergence to small- */
+      laguer (ad,j,&x,&its);                     /* est remaining root, and find the root. */
+      if (fabs (x.im) <= 2.0*EPS*fabs (x.re))
+	x.im=0.0;
+      roots[j]=x;
+      b=ad[j];                                  /* Forward deflation. */
+      for (jj=j-1;jj>=0;jj--)
+	{
+	  c=ad[jj];
+	  ad[jj]=b;
+	  b=gsl_complex_add (gsl_complex_mul (x,b),c);
+	}
+    }
+  if (polish)
+    for (j=1;j<=m;j++)                        /* Polish the roots using the undeflated coeffi- */
+      laguer (a,m,&roots[j],&its);            /* cients. */
+  for (j=2;j<=m;j++)                          /* Sort roots by their real parts by straight insertion */
+    {
+      x=roots[j];
+      for (i=j-1;i>=1;i--) {
+	if (roots[i].re <= x.re)
+	  break;
+	roots[i+1]=roots[i];
+      }
+      roots[i+1]=x;
+    }
+}
+
+#define ITMAX 20	/* At most ITMAX iterations. */
+#define TINY 2.0-15	/* TIMJ, was (float): 1.0e-6 */
+
+static void
+qroot (double p[], int n, double *b, double *c, double eps)
+     /* Given n+1 coefficients p[0..n] of a polynomial of degree n, and trial values for the coefficients
+	of a quadratic factor x*x+b*x+c, improve the solution until the coefficients b,c change by less
+	than eps. The routine poldiv §5.3 is used. */
+{
+  int iter;
+  double sc,sb,s,rc,rb,r,dv,delc,delb;
+  double *q,*qq,*rem;
+  double d[3];
+  
+  q=vector (0,n);
+  qq=vector (0,n);
+  rem=vector (0,n);
+  d[2]=1.0;
+  for (iter=1;iter<=ITMAX;iter++)
+    {
+      d[1]=(*b);
+      d[0]=(*c);
+      poldiv (p,n,d,2,q,rem);
+      s=rem[0];                                        /* First division r,s. */
+      r=rem[1];
+      poldiv (q,(n-1),d,2,qq,rem);
+      sb = -(*c)*(rc = -rem[1]);                       /* Second division partial r,s with respect to */
+      rb = -(*b)*rc+(sc = -rem[0]);                    /* c. */
+      dv=1.0/(sb*rc-sc*rb);                            /* Solve 2x2 equation. */
+      delb=(r*sc-s*rc)*dv;
+      delc=(-r*sb+s*rb)*dv;
+      *b += (delb=(r*sc-s*rc)*dv);
+      *c += (delc=(-r*sb+s*rb)*dv);
+      if ((fabs (delb) <= eps*fabs (*b) || fabs (*b) < TINY)
+	  && (fabs (delc) <= eps*fabs (*c) || fabs (*c) < TINY))
+	{
+	  free_vector (rem,0,n);                      /* Coefficients converged. */
+	  free_vector (qq,0,n);
+	  free_vector (q,0,n);
+	  return;
+	}
+    }
+  nrerror ("Too many iterations in routine qroot");
+}
+
+#define SNCNDN_CA 0.0003                      /* The accuracy is the square of SNCNDN_CA. */
+static void
+sncndn (double uu, double emmc, double *sn_p, double *cn_p, double *dn_p)
+     /* Returns the Jacobian elliptic functions sn(u, kc), cn(u, kc), and dn(u, kc). Here uu = u, while
+	emmc = k2c. */
+{
+  double a,b,c,d,emc,u,sn,cn,dn;
+  double em[14],en[14];
+  int i,ii,l,bo;
+  d=0; /* TIMJ: shutup compiler */
+  
+  emc=emmc;
+  u=uu;
+  if (emc) {
+    bo=(emc < 0.0);
+    if (bo) {
+      d=1.0-emc;
+      emc /= -1.0/d;
+      u *= (d=sqrt(d));
+    }a=1.0;
+    dn=1.0;
+    for (i=1;i<=13;i++) {
+      l=i;
+      em[i]=a;
+      en[i]=(emc=sqrt(emc));
+      c=0.5*(a+emc);
+      if (fabs(a-emc) <= SNCNDN_CA*a) break;
+      emc *= a;
+      a=c;
+    }u *= c;
+    sn=sin(u);
+    cn=cos(u);
+    if (sn) {
+      a=cn/sn;
+      c *= a;
+      for (ii=l;ii>=1;ii--) {
+	b=em[ii];
+	a *= c;
+	c *= dn;
+	dn=(en[ii]+a)/(b+a);
+	a=c/b;
+      }a=1.0/sqrt(c*c+1.0);
+      sn=(sn >= 0.0 ? a : -a);
+      cn=c*sn;
+    }if (bo) {
+      a=dn;
+      dn=cn;
+      cn=a;
+      sn /= d;
+    }
+  } else {
+    cn=1.0/cosh(u);
+    dn=cn;
+    sn=tanh(u);
+  }
+  if (sn_p)
+    *sn_p = sn;
+  if (cn_p)
+    *cn_p = cn;
+  if (dn_p)
+    *dn_p = dn;
+}
+
+static void
+sncndnC (GslComplex uu, GslComplex emmc, GslComplex *sn_p, GslComplex *cn_p, GslComplex *dn_p)
+{
+  GslComplex a,b,c,d,emc,u,sn,cn,dn;
+  GslComplex em[14],en[14];
+  int i,ii,l,bo;
+  
+  emc=emmc;
+  u=uu;
+  if (emc.re || emc.im) /* gsl_complex_abs (emc)) */
+    {
+      /* bo=gsl_complex_abs (emc) < 0.0; */
+      bo=emc.re < 0.0;
+      if (bo) {
+	d=gsl_complex_sub (ONE, emc);
+	emc = gsl_complex_div (emc, gsl_complex_div (gsl_complex (-1.0, 0), d));
+	d = gsl_complex_sqrt (d);
+	u = gsl_complex_mul (u, d);
+      }
+      a=ONE; dn=ONE;
+      for (i=1;i<=13;i++) {
+	l=i;
+	em[i]=a;
+	emc = gsl_complex_sqrt (emc);
+	en[i]=emc;
+	c = gsl_complex_mul (gsl_complex (0.5, 0), gsl_complex_add (a, emc));
+	if (gsl_complex_abs (gsl_complex_sub (a, emc)) <=
+	    gsl_complex_abs (gsl_complex_mul (gsl_complex (SNCNDN_CA, 0), a)))
+	  break;
+	emc = gsl_complex_mul (emc, a);
+	a=c;
+      }
+      u = gsl_complex_mul (u, c);
+      sn = gsl_complex_sin (u);
+      cn = gsl_complex_cos (u);
+      if (sn.re) /* gsl_complex_abs (sn)) */
+	{
+	  a= gsl_complex_div (cn, sn);
+	  c = gsl_complex_mul (c, a);
+	  for (ii=l;ii>=1;ii--) {
+	    b = em[ii];
+	    a = gsl_complex_mul (a, c);
+	    c = gsl_complex_mul (c, dn);
+	    dn = gsl_complex_div (gsl_complex_add (en[ii], a), gsl_complex_add (b, a));
+	    a = gsl_complex_div (c, b);
+	  }
+	  a = gsl_complex_div (ONE, gsl_complex_sqrt (gsl_complex_add (ONE, gsl_complex_mul (c, c))));
+	  if (sn.re >= 0.0) /* gsl_complex_arg (sn) >= 0.0) */
+	    sn = a;
+	  else
+	    {
+	      sn.re = -a.re;
+	      sn.im = a.im;
+	    }
+	  cn = gsl_complex_mul (c, sn);
+	}
+      if (bo) {
+	a=dn;
+	dn=cn;
+	cn=a;
+	sn = gsl_complex_div (sn, d);
+      }
+    } else {
+      cn=gsl_complex_div (ONE, gsl_complex_cosh (u));
+      dn=cn;
+      sn=gsl_complex_tanh (u);
+    }
+  if (sn_p)
+    *sn_p = sn;
+  if (cn_p)
+    *cn_p = cn;
+  if (dn_p)
+    *dn_p = dn;
+}
+
+#define RF_ERRTOL	0.0025		/* TIMJ, was(float): 0.08 */
+#define RF_TINY		2.2e-307	/* TIMJ, was(float): 1.5e-38 */
+#define RF_BIG		1.5e+307	/* TIMJ, was(float):  3.0e37 */
+#define RF_THIRD	(1.0/3.0)
+#define RF_C1		(1.0/24.0)
+#define RF_C2		0.1
+#define RF_C3		(3.0/44.0)
+#define RF_C4		(1.0/14.0)
+
+static double
+rf (double x, double y, double z)
+     /* Computes Carlson's elliptic integral of the first kind, RF (x, y, z). x, y, and z must be nonneg-
+	ative, and at most one can be zero. RF_TINY must be at least 5 times the machine underflow limit,
+	RF_BIG at most one fifth the machine overflow limit. */
+{
+  double alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt;
+  
+  if (1 /* TIMJ: add verbose checks */)
+    {
+      if (DMIN (DMIN (x, y), z) < 0.0)
+	nrerror ("rf: x,y,z have to be positive");
+      if (DMIN (DMIN (x + y, x + z), y + z) < RF_TINY)
+	nrerror ("rf: only one of x,y,z may be 0");
+      if (DMAX (DMAX (x, y), z) > RF_BIG)
+	nrerror ("rf: at least one of x,y,z is too big");
+    }
+  if (DMIN(DMIN(x,y),z) < 0.0 || DMIN(DMIN(x+y,x+z),y+z) < RF_TINY ||
+      DMAX(DMAX(x,y),z) > RF_BIG)
+    nrerror("invalid arguments in rf");
+  xt=x;
+  yt=y;
+  zt=z;
+  do {
+    sqrtx=sqrt(xt);
+    sqrty=sqrt(yt);
+    sqrtz=sqrt(zt);
+    alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
+    xt=0.25*(xt+alamb);
+    yt=0.25*(yt+alamb);
+    zt=0.25*(zt+alamb);
+    ave=RF_THIRD*(xt+yt+zt);
+    delx=(ave-xt)/ave;
+    dely=(ave-yt)/ave;
+    delz=(ave-zt)/ave;
+  } while (DMAX(DMAX(fabs(delx),fabs(dely)),fabs(delz)) > RF_ERRTOL);
+  e2=delx*dely-delz*delz;
+  e3=delx*dely*delz;
+  return (1.0+(RF_C1*e2-RF_C2-RF_C3*e3)*e2+RF_C4*e3)/sqrt(ave);
+}
+
+static GslComplex
+rfC (GslComplex x, GslComplex y, GslComplex z)
+{
+  GslComplex alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt;
+  GslComplex RFC_C1 = {1.0/24.0, 0}, RFC_C2 = {0.1, 0}, RFC_C3 = {3.0/44.0, 0}, RFC_C4 = {1.0/14.0, 0};
+  
+  if (DMIN (DMIN (gsl_complex_abs (x), gsl_complex_abs (y)), gsl_complex_abs (z)) < 0.0)
+    nrerror ("rf: x,y,z have to be positive");
+  if (DMIN (DMIN (gsl_complex_abs (x) + gsl_complex_abs (y), gsl_complex_abs (x) + gsl_complex_abs (z)),
+	    gsl_complex_abs (y) + gsl_complex_abs (z)) < RF_TINY)
+    nrerror ("rf: only one of x,y,z may be 0");
+  if (DMAX (DMAX (gsl_complex_abs (x), gsl_complex_abs (y)), gsl_complex_abs (z)) > RF_BIG)
+    nrerror ("rf: at least one of x,y,z is too big");
+  xt=x;
+  yt=y;
+  zt=z;
+  do {
+    sqrtx = gsl_complex_sqrt (xt);
+    sqrty = gsl_complex_sqrt (yt);
+    sqrtz = gsl_complex_sqrt (zt);
+    alamb = gsl_complex_add (gsl_complex_mul (sqrtx, gsl_complex_add (sqrty, sqrtz)), gsl_complex_mul (sqrty, sqrtz));
+    xt =    gsl_complex_mul (gsl_complex (0.25, 0), gsl_complex_add (xt, alamb));
+    yt =    gsl_complex_mul (gsl_complex (0.25, 0), gsl_complex_add (yt, alamb));
+    zt =    gsl_complex_mul (gsl_complex (0.25, 0), gsl_complex_add (zt, alamb));
+    ave =   gsl_complex_mul (gsl_complex (RF_THIRD, 0), gsl_complex_add3 (xt, yt, zt));
+    delx =  gsl_complex_div (gsl_complex_sub (ave, xt), ave);
+    dely =  gsl_complex_div (gsl_complex_sub (ave, yt), ave);
+    delz =  gsl_complex_div (gsl_complex_sub (ave, zt), ave);
+    /* } while (DMAX (DMAX (fabs (delx.re), fabs (dely.re)), fabs (delz.re)) > RF_ERRTOL); */
+  } while (DMAX (DMAX (gsl_complex_abs (delx), gsl_complex_abs (dely)), gsl_complex_abs (delz)) > RF_ERRTOL);
+  e2 = gsl_complex_sub (gsl_complex_mul (delx, dely), gsl_complex_mul (delz, delz));
+  e3 = gsl_complex_mul3 (delx, dely, delz);
+  return gsl_complex_div (gsl_complex_add3 (gsl_complex (1.0, 0),
+					    gsl_complex_mul (e2,
+							     gsl_complex_sub3 (gsl_complex_mul (RFC_C1, e2),
+									       RFC_C2,
+									       gsl_complex_mul (RFC_C3, e3))),
+					    gsl_complex_mul (RFC_C4, e3)),
+			  gsl_complex_sqrt (ave));
+}
+
+
+static double
+ellf (double phi, double ak)
+     /* Legendre elliptic integral of the 1st kind F(phi, k), evaluated using Carlson's function RF.
+	The argument ranges are 0 <= phi <= pi/2, 0 <=  k*sin(phi) <= 1. */
+{
+  double s=sin(phi);
+  return s*rf(DSQR(cos(phi)),(1.0-s*ak)*(1.0+s*ak),1.0);
+}