summaryrefslogtreecommitdiffstats
path: root/flow/gsl/gslmath.c
diff options
context:
space:
mode:
Diffstat (limited to 'flow/gsl/gslmath.c')
-rw-r--r--flow/gsl/gslmath.c1085
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);
+}