1 files changed, 113 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qpsrt.c b/gsl-1.9/integration/qpsrt.c
new file mode 100644
index 0000000..519bab4
--- /dev/null
+++ b/gsl-1.9/integration/qpsrt.c
@@ -0,0 +1,113 @@
+/* integration/qpsrt.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 inline void 
+qpsrt (gsl_integration_workspace * workspace);
+
+static inline
+void qpsrt (gsl_integration_workspace * workspace)
+{
+  const size_t last = workspace->size - 1;
+  const size_t limit = workspace->limit;
+
+  double * elist = workspace->elist;
+  size_t * order = workspace->order;
+
+  double errmax ;
+  double errmin ;
+  int i, k, top;
+
+  size_t i_nrmax = workspace->nrmax;
+  size_t i_maxerr = order[i_nrmax] ;
+  
+  /* Check whether the list contains more than two error estimates */
+
+  if (last < 2) 
+    {
+      order[0] = 0 ;
+      order[1] = 1 ;
+      workspace->i = i_maxerr ;
+      return ;
+    }
+
+  errmax = elist[i_maxerr] ;
+
+  /* This part of the routine is only executed if, due to a difficult
+     integrand, subdivision increased the error estimate. In the normal
+     case the insert procedure should start after the nrmax-th largest
+     error estimate. */
+
+  while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]]) 
+    {
+      order[i_nrmax] = order[i_nrmax - 1] ;
+      i_nrmax-- ;
+    } 
+
+  /* Compute the number of elements in the list to be maintained in
+     descending order. This number depends on the number of
+     subdivisions still allowed. */
+  
+  if(last < (limit/2 + 2)) 
+    {
+      top = last ;
+    }
+  else
+    {
+      top = limit - last + 1;
+    }
+  
+  /* Insert errmax by traversing the list top-down, starting
+     comparison from the element elist(order(i_nrmax+1)). */
+  
+  i = i_nrmax + 1 ;
+  
+  /* The order of the tests in the following line is important to
+     prevent a segmentation fault */
+
+  while (i < top && errmax < elist[order[i]])
+    {
+      order[i-1] = order[i] ;
+      i++ ;
+    }
+  
+  order[i-1] = i_maxerr ;
+  
+  /* Insert errmin by traversing the list bottom-up */
+  
+  errmin = elist[last] ;
+  
+  k = top - 1 ;
+  
+  while (k > i - 2 && errmin >= elist[order[k]])
+    {
+      order[k+1] = order[k] ;
+      k-- ;
+    }
+  
+  order[k+1] = last ;
+
+  /* Set i_max and e_max */
+
+  i_maxerr = order[i_nrmax] ;
+  
+  workspace->i = i_maxerr ;
+  workspace->nrmax = i_nrmax ;
+}
+
+