1 files changed, 132 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qc25c.c b/gsl-1.9/integration/qc25c.c
new file mode 100644
index 0000000..94f4b1f
--- /dev/null
+++ b/gsl-1.9/integration/qc25c.c
@@ -0,0 +1,132 @@
+/* integration/qc25c.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.
+ */
+
+struct fn_cauchy_params
+{
+  gsl_function *function;
+  double singularity;
+};
+
+static double fn_cauchy (double t, void *params);
+
+static void compute_moments (double cc, double *moment);
+
+static void
+qc25c (gsl_function * f, double a, double b, double c, 
+       double *result, double *abserr, int *err_reliable);
+
+static void
+qc25c (gsl_function * f, double a, double b, double c, 
+       double *result, double *abserr, int *err_reliable)
+{
+  double cc = (2 * c - b - a) / (b - a);
+
+  if (fabs (cc) > 1.1)
+    {
+      double resabs, resasc;
+
+      gsl_function weighted_function;
+      struct fn_cauchy_params fn_params;
+
+      fn_params.function = f;
+      fn_params.singularity = c;
+
+      weighted_function.function = &fn_cauchy;
+      weighted_function.params = &fn_params;
+
+      gsl_integration_qk15 (&weighted_function, a, b, result, abserr,
+                            &resabs, &resasc);
+      
+      if (*abserr == resasc)
+        {
+          *err_reliable = 0;
+        }
+      else 
+        {
+          *err_reliable = 1;
+        }
+
+      return;
+    }
+  else
+    {
+      double cheb12[13], cheb24[25], moment[25];
+      double res12 = 0, res24 = 0;
+      size_t i;
+      gsl_integration_qcheb (f, a, b, cheb12, cheb24);
+      compute_moments (cc, moment);
+
+      for (i = 0; i < 13; i++)
+        {
+          res12 += cheb12[i] * moment[i];
+        }
+
+      for (i = 0; i < 25; i++)
+        {
+          res24 += cheb24[i] * moment[i];
+        }
+
+      *result = res24;
+      *abserr = fabs(res24 - res12) ;
+      *err_reliable = 0;
+
+      return;
+    }
+}
+
+static double
+fn_cauchy (double x, void *params)
+{
+  struct fn_cauchy_params *p = (struct fn_cauchy_params *) params;
+  gsl_function *f = p->function;
+  double c = p->singularity;
+  return GSL_FN_EVAL (f, x) / (x - c);
+}
+
+static void
+compute_moments (double cc, double *moment)
+{
+  size_t k;
+
+  double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
+  double a1 = 2 + a0 * cc;
+
+  moment[0] = a0;
+  moment[1] = a1;
+
+  for (k = 2; k < 25; k++)
+    {
+      double a2;
+
+      if ((k % 2) == 0)
+        {
+          a2 = 2.0 * cc * a1 - a0;
+        }
+      else
+        {
+          const double km1 = k - 1.0;
+          a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
+        }
+
+      moment[k] = a2;
+
+      a0 = a1;
+      a1 = a2;
+    }
+}