text:
#ifndef __SIMQ_H
#define __SIMQ_H

#include "Queue.h"
#include "statistics.h"
#include "input.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* Defining non-queue function prototypes */
Customer* createCustomer(gsl_rng **);
gsl_rng * randomGenerator();
void validateArgs(int, char **);
#endif

ext: h


text:
#include <pygsl/solver.h>
#include <gsl/gsl_multifit_nlin.h>

const char  * filename = __FILE__;
PyObject *module = NULL;
static const char multifit_f_type_name[] = "F-MultiFitSolver";
static const char multifit_fdf_type_name[] = "FdF-MultiFitSolver";

static int
PyGSL_multifit_function_wrap(const gsl_vector *x, void *params, gsl_vector *f)
{
     PyGSL_solver *self = (PyGSL_solver *) params;
     return PyGSL_function_wrap_Op_On(x, f, self->cbs[0], self->args, x->size, f->size,
				      __FUNCTION__);
}

static int
PyGSL_multifit_function_wrap_df(const gsl_vector *x, void *params, gsl_matrix *df)
{
     PyGSL_solver *self = (PyGSL_solver *) params;
     /* size 1 or size 2 from matrix ? */
     return PyGSL_function_wrap_Op_Opn(x, df, self->cbs[1], self->args, df->size1, x->size,
				       __FUNCTION__);
}

static int
PyGSL_multifit_function_wrap_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *df)
{
     PyGSL_solver *self = (PyGSL_solver *) params;
     /* size 1 or size 2 from matrix ? */
     return PyGSL_function_wrap_Op_On_Opn(x, f, df, self->cbs[2], self->args, f->size, x->size,
					  __FUNCTION__);
}

static PyObject *
PyGSL_multifit_fdfsolver_set(PyGSL_solver *self, PyObject *pyargs, PyObject *kw)
{

     gsl_multifit_function_fdf * c_sys;
     struct pygsl_solver_n_set info = {1, NULL, (set_m_t)gsl_multifit_fdfsolver_set};
     PyObject * tmp;

     FUNC_MESS_BEGIN();
     if(self->c_sys == NULL){
	  if((c_sys = calloc(1, sizeof(gsl_multifit_function_fdf))) == NULL){
	       PyGSL_ERROR_NULL("Could not allocate the memory for the c_sys", GSL_ENOMEM);
	  }
	  c_sys->n=self->problem_dimensions[1];
	  c_sys->p=self->problem_dimensions[0];
	  c_sys->f  = PyGSL_multifit_function_wrap;
	  c_sys->df = PyGSL_multifit_function_wrap_df;
	  c_sys->fdf = PyGSL_multifit_function_wrap_fdf;
	  c_sys->params=(void*)self;
	  info.c_sys = c_sys;
     }else{
	  info.c_sys = self->c_sys;
     }
     tmp =  PyGSL_solver_n_set(self, pyargs, kw, &info);
     if(tmp == NULL){
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__ - 2);
     }
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_multifit_fdfsolver_position(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_vec(self, args, (ret_vec)gsl_multifit_fdfsolver_position);
}
static gsl_multifit_fdfsolver *
PyGSL_get_multifit_solver(PyGSL_solver *self)
{
     FUNC_MESS_BEGIN();
     assert(PyGSL_solver_check(self));
     FUNC_MESS_END();
     return (gsl_multifit_fdfsolver *) (self->solver);
}

static PyObject*
PyGSL_multifit_fdfsolver_x(PyGSL_solver *self, PyObject *args)
{
     return (PyObject *) PyGSL_copy_gslvector_to_pyarray(PyGSL_get_multifit_solver(self)->x);
}

static PyObject*
PyGSL_multifit_fdfsolver_dx(PyGSL_solver *self, PyObject *args)
{
     return (PyObject *) PyGSL_copy_gslvector_to_pyarray(PyGSL_get_multifit_solver(self)->dx);
}

static PyObject*
PyGSL_multifit_fdfsolver_f(PyGSL_solver *self, PyObject *args)
{
     return (PyObject *) PyGSL_copy_gslvector_to_pyarray(PyGSL_get_multifit_solver(self)->f);
}

static PyObject*
PyGSL_multifit_fdfsolver_J(PyGSL_solver *self, PyObject *args)
{
     return (PyObject *) PyGSL_copy_gslmatrix_to_pyarray(PyGSL_get_multifit_solver(self)->J);
}

static PyObject*
PyGSL_multifit_fdfsolver_test_delta(PyGSL_solver *self, PyObject *args)
{
     int flag;
     double epsabs, epsrel;
     gsl_multifit_fdfsolver *s = self->solver;
     if(!PyArg_ParseTuple(args, "dd", &epsabs, &epsrel))
	  return NULL;
     flag = gsl_multifit_test_delta(s->dx, s->x, epsabs, epsrel);
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}

static PyObject*
PyGSL_multifit_fdfsolver_test_gradient(PyGSL_solver *self, PyObject *args)
{
     int flag;
     double epsabs;
     gsl_vector *g = NULL;
     gsl_multifit_fdfsolver *s = self->solver;

     if(!PyArg_ParseTuple(args, "d", &epsabs))
	  return NULL;
     flag = gsl_multifit_gradient(s->J, s->f, g);
     if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS)
	  return NULL;
     flag = gsl_multifit_test_gradient(g, epsabs);
     return PyGSL_ERROR_FLAG_TO_PYINT(flag);
}

static PyMethodDef PyGSL_multifit_fmethods[] = {
     {NULL, NULL, 0, NULL}
};

static PyMethodDef PyGSL_multifit_fdfmethods[] = {
     {"J",    (PyCFunction)PyGSL_multifit_fdfsolver_J,   METH_NOARGS, NULL},
     {"dx",   (PyCFunction)PyGSL_multifit_fdfsolver_dx,   METH_NOARGS, NULL},
     {"x",    (PyCFunction)PyGSL_multifit_fdfsolver_x,   METH_NOARGS, NULL},
     {"position",    (PyCFunction)PyGSL_multifit_fdfsolver_position,   METH_NOARGS, NULL},
     {"f",    (PyCFunction)PyGSL_multifit_fdfsolver_f,    METH_NOARGS, NULL},
     {"set",  (PyCFunction)PyGSL_multifit_fdfsolver_set,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"test_delta",  (PyCFunction)PyGSL_multifit_fdfsolver_test_delta,  METH_VARARGS, NULL},
     {"test_gradient",  (PyCFunction)PyGSL_multifit_fdfsolver_test_gradient,  METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

const struct _SolverStatic
multifit_solver_f   = {{ (void_m_t) gsl_multifit_fsolver_free,
			 /* gsl_multifit_fsolver_restart */  (void_m_t) NULL,
			 (name_m_t) gsl_multifit_fsolver_name,
			 (int_m_t) gsl_multifit_fsolver_iterate},
		       1, PyGSL_multifit_fmethods, multifit_f_type_name},
multifit_solver_fdf = {{(void_m_t) gsl_multifit_fdfsolver_free,
		     /* gsl_multifit_fdfsolver_restart (void_m_t) */ NULL,
		     (name_m_t) gsl_multifit_fdfsolver_name,
		     (int_m_t)  gsl_multifit_fdfsolver_iterate},
		       3, PyGSL_multifit_fdfmethods, multifit_fdf_type_name};

static PyObject*
PyGSL_multifit_f_init(PyObject *self, PyObject *args,
		      const gsl_multifit_fsolver_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_multifit_fsolver_alloc,
			      &multifit_solver_f};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 2);
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_multifit_fdf_init(PyObject *self, PyObject *args,
		      const gsl_multifit_fdfsolver_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_multifit_fdfsolver_alloc,
			      &multifit_solver_fdf};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 2);
     FUNC_MESS_END();
     return tmp;
}

#define AMFIT_FDF(name)                                                  \
static PyObject* PyGSL_multifit_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_multifit_fdf_init(self, args,  gsl_multifit_fdfsolver_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}

AMFIT_FDF(lmsder)
AMFIT_FDF(lmder)

PyObject *
PyGSL_multifit_gradient(PyObject *self, PyObject *args)
{
  PyArrayObject *J_a = NULL, *f_a = NULL, *g_a = NULL;
  PyObject *J_o = NULL, *f_o = NULL;
  gsl_vector_view f;
  gsl_vector_view g;
  gsl_matrix_view J;

  PyGSL_array_index_t stride_recalc, dimension;
  int flag;

  if(!PyArg_ParseTuple(args, "OO:gsl_multifit_gradient", &J_o, &f_o)){
       return NULL;
  }

  J_a = PyGSL_matrix_check(J_o, -1, -1, PyGSL_DARRAY_CINPUT(1), NULL, NULL, NULL);
  if(J_a == NULL) goto fail;

  dimension = J_a->dimensions[0];
  /* Numpy calculates strides in bytes, gsl in basis type */
  f_a = PyGSL_vector_check(f_o, dimension, PyGSL_DARRAY_INPUT(2), &stride_recalc, NULL);
  if(f_a == NULL) goto fail;

  dimension = J_a->dimensions[1];
  g_a = (PyArrayObject *) PyGSL_New_Array(1, &dimension, PyArray_DOUBLE);
  if(g_a == NULL) goto fail;

  J = gsl_matrix_view_array((double *) J_a->data, J_a->dimensions[0], J_a->dimensions[1]);
  f = gsl_vector_view_array_with_stride((double *) f_a->data, stride_recalc, f_a->dimensions[0]);
  g = gsl_vector_view_array((double *) g_a->data, dimension);
  flag = gsl_multifit_gradient(&J.matrix, &f.vector, &g.vector);

  Py_DECREF(J_a);
  Py_DECREF(f_a);

  if((PyGSL_ERROR_FLAG(flag)) != GSL_SUCCESS)
       goto fail;

  return (PyObject * )g_a;

  fail :
    Py_XDECREF(J_a);
    Py_XDECREF(f_a);
    Py_XDECREF(g_a);
    return NULL;
}

PyObject *
PyGSL_multifit_covar(PyObject *self, PyObject *args)
{
  PyArrayObject *J_a = NULL, *C_a = NULL;
  PyObject *J_o = NULL;
  gsl_matrix_view J, C;
  PyGSL_array_index_t dimensions[2];
  int flag;
  double epsrel;

  if(!PyArg_ParseTuple(args, "Od:gsl_multifit_covar", &J_o, &epsrel)){
    return NULL;
  }

  J_a = PyGSL_matrix_check(J_o, -1, -1, PyGSL_DARRAY_CINPUT(1), NULL, NULL, NULL);
  if(J_a == NULL) goto fail;

  dimensions[0] = J_a->dimensions[1];
  dimensions[1] = J_a->dimensions[1];
  C_a = (PyArrayObject *) PyGSL_New_Array(2, dimensions, PyArray_DOUBLE);
  if(C_a == NULL) goto fail;

  J = gsl_matrix_view_array((double *) J_a->data, J_a->dimensions[0], J_a->dimensions[1]);
  C = gsl_matrix_view_array((double *) C_a->data, C_a->dimensions[0], C_a->dimensions[1]);

  flag = gsl_multifit_covar(&J.matrix, epsrel, &C.matrix);

  Py_DECREF(J_a);
  if((PyGSL_ERROR_FLAG(flag)) != GSL_SUCCESS)
       goto fail;
  return (PyObject * )C_a;

  fail :
    Py_XDECREF(J_a);
    Py_XDECREF(C_a);
    return NULL;
}

static PyObject *
PyGSL_multifit_test_delta(PyObject * self, PyObject * args)
{
     return PyGSL_solver_vvdd_i(self, args, gsl_multifit_test_delta);
}

static PyObject *
PyGSL_multifit_test_gradient(PyObject * self, PyObject * args)
{
     return PyGSL_solver_vd_i(self, args, gsl_multifit_test_gradient);
}

static PyMethodDef mMethods[] = {
     /* multifit solvers */
     {"lmder",          PyGSL_multifit_init_lmder,  METH_VARARGS, NULL},
     {"lmsder",          PyGSL_multifit_init_lmsder,  METH_VARARGS, NULL},
     /* multifit funcs */
     {"fit_test_delta",    PyGSL_multifit_test_delta,     METH_VARARGS, NULL},
     {"fit_test_gradient", PyGSL_multifit_test_gradient,  METH_VARARGS, NULL},
     {"gradient",          PyGSL_multifit_gradient,  METH_VARARGS, NULL},
     {"covar",             PyGSL_multifit_covar,  METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

static const char PyGSL_multifit_nlin_module_doc[] = "XXX Missing \n";
void
initmultifit_nlin(void)
{
     PyObject* m, *dict, *item;
     FUNC_MESS_BEGIN();

     m=Py_InitModule("multifit_nlin", mMethods);
     module = m;
     assert(m);
     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     init_pygsl()
     import_pygsl_solver();
     assert(PyGSL_API);

     if (!(item = PyString_FromString((char*)PyGSL_multifit_nlin_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }
     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     FUNC_MESS_END();
     return;

 fail:
     FUNC_MESS("FAIL");
     return;
}

ext: c

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* movstat/fill.c
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * 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 3 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.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_movstat.h>

/*
gsl_movstat_fill()
  Fill window for sample 'idx' from x using given end conditions

Inputs: endtype     - how to handle end points
        x           - input vector, length n
        idx         - index of center sample in window in [0,n-1]
        H           - number of samples left of center to include
        J           - number of samples right of center to include
        window      - (output) window of samples centered on x_{idx},
                      W_{idx}^{H,J}, length H + J + 1

Return: size of filled window (<= H + J + 1)
*/

size_t
gsl_movstat_fill(const gsl_movstat_end_t endtype, const gsl_vector * x, const size_t idx,
                 const size_t H, const size_t J, double * window)
{
  if (idx >= x->size)
    {
      GSL_ERROR_VAL ("window center index must be between 0 and n - 1", GSL_EDOM, 0);
    }
  else
    {
      const int n = (int) x->size;
      const int iidx = (int) idx;
      const int iH = (int) H;
      const int iJ = (int) J;
      int idx1, idx2, j;
      size_t window_size;

      if (endtype == GSL_MOVSTAT_END_TRUNCATE)
        {
          idx1 = GSL_MAX(iidx - iH, 0);
          idx2 = GSL_MIN(iidx + iJ, n - 1);
        }
      else
        {
          idx1 = iidx - iH;
          idx2 = iidx + iJ;
        }

      window_size = (size_t) (idx2 - idx1 + 1);

      /* fill sliding window */
      for (j = idx1; j <= idx2; ++j)
        {
          int widx = j - idx1;

          if (j < 0)
            {
              /* initial condition */
              if (endtype == GSL_MOVSTAT_END_PADZERO)
                window[widx] = 0.0;
              else if (endtype == GSL_MOVSTAT_END_PADVALUE)
                window[widx] = gsl_vector_get(x, 0);
            }
          else if (j >= n)
            {
              if (endtype == GSL_MOVSTAT_END_PADZERO)
                window[widx] = 0.0;
              else if (endtype == GSL_MOVSTAT_END_PADVALUE)
                window[widx] = gsl_vector_get(x, n - 1);
            }
          else
            {
              window[widx] = gsl_vector_get(x, j);
            }
        }

      return window_size;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* PETSC_DEPTRECATED causes parse erros in bindgen */
#define PETSC_DEPRECATED(arg)
/* Not sure how this is generated... */
#define PETSC_BITS_PER_BYTE 8
#include <petsc.h>

max_stars_repo_name: tkonolige/petsc-bindgen-rust

max_issues_repo_name: tkonolige/petsc-bindgen-rust

max_forks_repo_name: tkonolige/petsc-bindgen-rust


text:
/*
 *  CosmoCalcs.c
 *  Library for computing physical parameters in lambda CDM cosmologies
 *  FRLW universes.
 *
 *  Created by Sean Lake on 5/24/14.
 */

#include "CosmoCalcs.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <gsl/gsl_poly.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

#define MAX(a, b) \
	({ __typeof__ (a) x = (a); \
		__typeof__ (b) y = (b); \
		x > y ? x : y; })
#define MIN(a, b) \
	({ __typeof__ (a) x = (a); \
		__typeof__ (b) y = (b); \
		x <= y ? x : y; })
#define CLAMP(a, min, max) \
	({ __typeof__ (a) hilim = (max); \
		__typeof__ (a) lolim = (min); \
		__typeof__ (a) x = (a); \
		x = x <= hilim ? x : hilim; \
		x > lolim ? x : lolim; })
#define EXPEL(a, leftbound, rightbound) \
	({ __typeof__ (a) l = (leftbound); \
		__typeof__ (a) r = (rightbound); \
		__typeof__ (a) x = (a); \
		__typeof__ (a) Two = 2; \
		x >= r || x <= l ? x : ((Two * (x - leftbound) >= \
			rightbound - leftbound) ? rightbound : leftbound); })

#ifdef FP_FAST_FMA
	#define FMA_m(a, b, c) fma(a, b, c)
#else
	#define FMA_m(a, b, c) ( (a) * (b) + (c) )
#endif

//Function prototypes (because we don't export everything)
static double NotImplementedFunc(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr);
static double NotImplementedInvFunc(struct UniverseLCDM *uni, double Dc, unsigned int branchNum,
	double epsrel, unsigned int maxiter);

//interpolation functions
static void LtEdgeBases(double, double, double *restrict);
static inline void RtEdgeBases(double, double, double *restrict);
static void centerBases(double, double, double *restrict);
static void Dc0Derivatives(UniverseLCDM *, double, double *);
static void tL0Derivatives(UniverseLCDM *, double, double *);
static double MeanValueInterval(double, double *, double *, double *);

//integration functions
static inline double Einv(UniverseLCDM *, double);
static inline double dDc0_dz_gsl(double, void *);
static inline double dtL0_dz_gsl(double, void *);

//Cosmology functions
static inline double DcT0Flat(struct UniverseLCDM *, double, bool, double, double, size_t, double *);
static double DcT0NegK(struct UniverseLCDM *, double, bool, double, double, size_t, double *);
static double DcT0PosK(struct UniverseLCDM *, double, bool, double, double, size_t, double *);
static inline double Vc0Flat(struct UniverseLCDM *, double, bool, double, double, size_t, double *);
static double Vc0NegK(struct UniverseLCDM *, double, bool, double, double, size_t, double *);
static double Vc0PosK(struct UniverseLCDM *, double, bool, double, double, size_t, double *);

//Inverse cosmology functions
static double FindDa0Max(UniverseLCDM *, double, unsigned int);
static double DcT0InvFlat(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double DcT0InvNegK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double DcT0InvPosK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Da0InvFlat(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Da0InvNegK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Da0InvPosK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Dl0InvFlat(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Dl0InvNegK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Dl0InvPosK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Vc0InvFlat(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Vc0InvNegK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);
static double Vc0InvPosK(struct UniverseLCDM *, double, unsigned int, double, unsigned int);

//Constants
const double clight = 299792458.0;
const double MeterPerAu = 0.149598e12;
const double AuPerMeter = 0x1.d662b65385936p-38; //1.0 / MeterPerAu;
const double MeterPerParsec = 0x1.b68064bc05bf1p+54; //MeterPerAu / tan( pi / (3600.0 * 180.0) );
const double ParsecPerMeter = 0x1.2ae8abdcb93c6p-55; //1.0 / MeterPerParsec;
const double MeterPerMpc = 0x1.a2300a115feaep+74; //1e6 * MeterPerParsec;
const double MpcPerMeter = 0x1.396dbd4dbf44bp-75; //1.0 / MpcPerMeter;
const double SecPerYear = 31557600.0; //365.25 * 24.0 * 3600.0;
const double YearPerSec = 0x1.1032d78f1540bp-25; //1.0 / SecPerYear;
const double SecPerGyear = 31557600.0e9; //1e9 * SecPerYear;
const double GyearPerSec = 0x1.244561c42c152p-55; //1.0 / SecPerGyear;

const double Pi = 0x1.921fb54442d18p+1;

static double NotImplementedFunc(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	fprintf(stderr, "Error: calling a function that isn't implemented yet.");
	exit(EXIT_FAILURE);

	return NAN;
}

static double NotImplementedInvFunc(struct UniverseLCDM *uni, double Dc, unsigned int branchNum,
	double epsrel, unsigned int maxiter) {
	fprintf(stderr, "Error: calling an inverse function that isn't implemented yet.");
	exit(EXIT_FAILURE);

	return NAN;
}

/* Define the Hermite interpolation basis functions, with x the distance from
 * the central point in units of h, the distance between points.
 * Assumes equidistance points. results written to yvals, which must have length 4.*/

static void LtEdgeBases(double x, double h, double *restrict yvals) {
	double *res = yvals;
	double scale = x*x;
	scale *= scale;

	//zeroth derivative basis polynomial
	*res = FMA_m(FMA_m(FMA_m(105.0/32.0, x, -2.0), x, -385.0/32.0), x, 15.0/2.0);
	*res = FMA_m(FMA_m(FMA_m(*res, x, 495/32.0), x, -10.0), x, -231.0/32.0);
	*res = FMA_m(*res, x, 5.0) * scale;

	//first derivative basis polynomial
	++res;
	scale *= h;
	*res = FMA_m(FMA_m(FMA_m(41.0/32.0, x, -29.0/32.0), x, -145.0/32.0), x, 105.0/32.0);
	*res = FMA_m(FMA_m(FMA_m(*res, x, 175.0/32.0), x, -131.0/32.0), x, -71.0/32.0);
	*res = FMA_m(*res, x, 55.0/32.0) * scale;

	//second derivative basis polynomial
	++res;
	scale *= h;
	*res = FMA_m(FMA_m(FMA_m(3.0/16.0, x, -5.0/32.0), x, -5.0/8.0), x, 17.0/32.0);
	*res = FMA_m(FMA_m(FMA_m(*res, x, 11.0/16.0), x, -19.0/32.0), x, -1.0/4.0);
	*res = FMA_m(*res, x, 7.0/32.0) * scale;

	//third derivative basis polynomial
	++res;
	scale *= h;
	*res = FMA_m(FMA_m(FMA_m(1.0/96.0, x, -1.0/96.0), x, -1.0/32.0), x, 1.0/32.0);
	*res = FMA_m(FMA_m(FMA_m(*res, x, 1.0/32.0), x, -1.0/32.0), x, -1.0/96.0);
	*res = FMA_m(*res, x, 1.0/96.0) * scale;

	return;
}

static inline void RtEdgeBases(double x, double h, double *restrict yvals) {
	LtEdgeBases(-x, h, yvals);
	yvals[1] *= -1.0;
	yvals[3] *= -1.0;

	return;
}

static void centerBases(double x, double h, double *restrict yvals) {
	double *res = yvals;
	double x2 = x*x;
	double scale = h;

	//zeroth derivative basis polynomial
	*res = FMA_m(FMA_m(FMA_m(4.0, x2, -15.0), x2, 20.0), x2, -10.0);
	*res = FMA_m(*res, x2*x2, 1.0);

	//first derivative basis polynomial
	++res;
	*res = (*(res-1)) * x * h;

	//second derivative basis polynomial
	++res;
	scale *= h;
	*res = FMA_m(FMA_m(FMA_m(0.5, x2, -2.0), x2, 3.0), x2, -2.0);
	*res = FMA_m(*res, x2, 0.5) * x2 * scale;

	//third derivative basis polynomial
	++res;
	*res = (*(res-1)) * x * h / 3.0;

	return;
}

//Computes Einverse and its first two derivatives
static void Dc0Derivatives(UniverseLCDM *uni, double invscale, double *results) {
	double rootarg, E, scale;
	double *res, *coeff;

	if ( uni->flat ) {
		rootarg = FMA_m(uni->OmegaRelativistic, invscale, uni->OmegaMatter);
		rootarg = FMA_m(rootarg, invscale * invscale * invscale, uni->OmegaLambda );
	}
	else {
		rootarg = FMA_m(uni->OmegaRelativistic, invscale, uni->OmegaMatter);
		rootarg = FMA_m(rootarg, invscale, uni->OmegaCurvature);
		rootarg = FMA_m(rootarg, invscale * invscale, uni->OmegaLambda );
	}

	E = sqrt(rootarg);
	res = results;

	//First order derivative
	scale = E / rootarg;
	*res = scale;

	//Second order derivative
	++res;
	coeff = uni->Dc0_2ndDerCoeffs;
	*res = FMA_m(FMA_m(coeff[2], invscale, coeff[1]), invscale, coeff[0]) * invscale;
	scale /= rootarg;
	*res *= scale;

	//Third order derivative
	++res;
	coeff = uni->Dc0_3rdDerCoeffs;
	*res = FMA_m(FMA_m(FMA_m(coeff[6], invscale, coeff[5]), invscale, coeff[4]), invscale, coeff[3]);
	*res = FMA_m(FMA_m(FMA_m(*res, invscale, coeff[2]), invscale, coeff[1]), invscale, coeff[0]);
	scale /= rootarg;
	*res *= scale;

	return;
}

//Computes the first four derivatives of tL0
static void tL0Derivatives(UniverseLCDM *uni, double invscale, double *results) {
	double rootarg, E, scale, scalestep;
	double *res, *coeffs;

	if ( uni->flat ) {
		rootarg = FMA_m(uni->OmegaRelativistic, invscale, uni->OmegaMatter);
		rootarg = FMA_m(rootarg, invscale * invscale * invscale, uni->OmegaLambda );
	}
	else {
		rootarg = FMA_m(uni->OmegaRelativistic, invscale, uni->OmegaMatter);
		rootarg = FMA_m(rootarg, invscale, uni->OmegaCurvature);
		rootarg = FMA_m(rootarg, invscale * invscale, uni->OmegaLambda );
	}

	E = sqrt(rootarg);
	res = results;

	//First order derivative
	scalestep = invscale * rootarg;
	scale = E / scalestep;
	*res = scale;

	//Second order derivative
	++res;
	coeffs = uni->tL0_2ndDerCoeffs;
	*res = FMA_m(FMA_m(coeffs[4], invscale, coeffs[3]), invscale, coeffs[2]);
	*res = FMA_m(FMA_m(*res, invscale, coeffs[1]), invscale, coeffs[0]);
	scale /= scalestep;
	*res *= scale;

	//Third order derivative
	++res;
	coeffs = uni->tL0_3rdDerCoeffs;
	*res = FMA_m(FMA_m(FMA_m(coeffs[8], invscale, coeffs[7]), invscale, coeffs[6]), invscale, coeffs[5]);
	*res = FMA_m(FMA_m(FMA_m(*res, invscale, coeffs[4]), invscale, coeffs[3]), invscale, coeffs[2]);
	*res = FMA_m(FMA_m(*res, invscale, coeffs[1]), invscale, coeffs[0]);
	scale /= scalestep;
	*res *= scale;

	return;
}

//Used for computing the cached integral values
static double MeanValueInterval(double sampleSpacing, double *leftFuncDirs,
						double *centerFuncDirs, double *rightFuncDirs) {
	double result;
	double x = sampleSpacing * sampleSpacing;

	result = (-FMA_m(-32.0, centerFuncDirs[2], leftFuncDirs[2] + rightFuncDirs[2] ) \
		* (1.0/630.0));

	result = FMA_m(result, x, \
		FMA_m(32.0, centerFuncDirs[0], 5.0*(leftFuncDirs[0] + rightFuncDirs[0])) \
		* (1.0/42.0));

	return result;
}

static double Einv(UniverseLCDM *uni, double invscale) {
	double rootarg;

	if ( uni->flat ) {
		rootarg = FMA_m(uni->OmegaRelativistic, invscale, uni->OmegaMatter);
		rootarg = FMA_m(rootarg, invscale * invscale * invscale, uni->OmegaLambda );
	}
	else {
		rootarg = FMA_m(uni->OmegaRelativistic, invscale, uni->OmegaMatter);
		rootarg = FMA_m(rootarg, invscale, uni->OmegaCurvature);
		rootarg = FMA_m(rootarg, invscale * invscale, uni->OmegaLambda );
	}

	return sqrt(rootarg) / rootarg;
}

static double dDc0_dz_gsl(double ainv, void *p) {
	return Einv((UniverseLCDM *)p, ainv);
}

double dDc0_dz_Cosmic(UniverseLCDM *uni, double z) {
	return Einv(uni, 1.0 + z);
}

double dDc_dz_Cosmic(UniverseLCDM *uni, double z) {
	return Einv(uni, 1.0 + z) * uni->DH;
}

double dtL0_dz_Cosmic(UniverseLCDM *uni, double z) {
	double invscale = 1.0 + z;
	return Einv(uni, invscale) / invscale;
}

static double dtL0_dz_gsl(double ainv, void *p) {
	return Einv((UniverseLCDM *)p, ainv) / ainv;
}

double dtL_dz_Cosmic(UniverseLCDM *uni, double z) {
	double invscale = 1.0 + z;
	return Einv(uni, invscale) / invscale * uni->tH;
}

bool InitUniverse(UniverseLCDM *uni, double H0, double OmegaL,
					double OmegaM, double OmegaR, bool flat, bool GSLcache ) {
	double OmegaK;

	uni->F_DcT0 = &NotImplementedFunc;
	uni->F_Vc0 = &NotImplementedFunc;
	uni->F_DcT0Inv = &NotImplementedInvFunc;
	uni->F_Vc0Inv = &NotImplementedInvFunc;

	if ( H0 < 0.0 || OmegaM < 0.0 || OmegaR < 0.0 ) {
		double *f;

		fprintf(stderr, "InitUniverse: unphysical universe.%s",
			" All of the following must be positive:\n");
		fprintf(stderr, "H0 = %.4g, OmegaM=%.4g, OmegaR=%.4g\n", H0, OmegaM, OmegaR);

		uni->OmegaLambda = NAN;
		uni->OmegaMatter = NAN;
		uni->OmegaRelativistic = NAN;
		uni->OmegaCurvature = NAN;
		uni->H0 = NAN;
		uni->h = NAN;
		uni->tH = NAN;
		uni->DH = NAN;
		uni->VH = NAN;
		uni->age0 = NAN;
		uni->zLastTurn = NAN;
		uni->zNextTurn = NAN;
		uni->CacheValid = false;

		for (f = uni->Dc0_2ndDerCoeffs; f < uni->Dc0_2ndDerCoeffs+3; ++f) {
			*f = NAN;
		}

		for (f = uni->Dc0_3rdDerCoeffs; f < uni->Dc0_3rdDerCoeffs+7; ++f) {
			*f = NAN;
		}

		for (f = uni->tL0_2ndDerCoeffs; f < uni->tL0_2ndDerCoeffs+5; ++f) {
			*f = NAN;
		}

		for (f = uni->tL0_3rdDerCoeffs; f < uni->tL0_3rdDerCoeffs+9; ++f) {
			*f = NAN;
		}

		for (f = uni->Dc0_TaylorCoeffs; f < uni->Dc0_TaylorCoeffs+5; ++f) {
			*f = NAN;
		}

		for (f = uni->tL0_TaylorCoeffs; f < uni->tL0_TaylorCoeffs+5; ++f) {
			*f = NAN;
		}

		return false;
	}

	uni->H0 = H0;
	uni->OmegaMatter = OmegaM;
	uni->OmegaRelativistic = OmegaR;
	uni->flat = flat;

	if ( flat == true ) {
		OmegaL = 1.0 - OmegaM - OmegaR;
		OmegaK = 0.0;
	}
	else {
		OmegaK = 1.0 - ( OmegaM + OmegaL + OmegaR );
	}
	uni->OmegaCurvature = OmegaK;
	uni->OmegaLambda = OmegaL;

	uni->h = H0 * 0.01 * 1e-3 * MeterPerMpc;
	uni->tH = 1.0 / H0;
	uni->DH = clight / H0;
	uni->VH = uni->DH * uni->DH * uni->DH;

	{
		//Find the turning points (if they exist)
		int polyorder;

		gsl_poly_complex_workspace *scratch;

		const double poly[5] = { OmegaL, 0.0, uni->OmegaCurvature, OmegaM, OmegaR };
		double roots[8];
		int err;
		double *iPart, *rPart;

		if (OmegaR > 0.0) {
			polyorder = 4;
		} else if (OmegaM > 0.0) {
			polyorder = 3;
		} else if (fabs(OmegaK) > NumericallyFlat) {
			polyorder = 2;
		} else {
			polyorder = 0;
		}

		if (polyorder > 0) {
			scratch = gsl_poly_complex_workspace_alloc( polyorder + 1 );
			err = gsl_poly_complex_solve( poly, polyorder + 1, scratch, roots );
			if ( err == GSL_EFAILED ) {
				fprintf(stderr, "GSL failed to find all of the candidate turning points%s",
					"for this cosmology:\n");
				fprintf(stderr, "OmegaL = %.4f, OmegaK = %.4f, OmegaM = %.4f, OmegaR = %.4f\n",
					OmegaL, uni->OmegaCurvature, OmegaM, OmegaR);

				return false;
			}
		}

		uni->zLastTurn = INFINITY;
		uni->zNextTurn = -1.0;

		for (rPart = roots; rPart < roots + 2*polyorder-1; rPart += 2) {
			iPart = rPart + 1;

			if (fabs(*iPart) < IMAGINARYTOL && *rPart > 0.0) {
				*rPart -= 1.0; //convert inverse scale to z

				if (*rPart < uni->zLastTurn && *rPart > 0.0) {
					uni->zLastTurn = *rPart;
				}
				else if (*rPart > uni->zNextTurn && *rPart <= 0.0) {
					uni->zNextTurn = *rPart;
				}
			}
		}
	}

	{
		//Init the derivative parameters
		double *f;
		double OmegaK = uni->OmegaCurvature;

		//Dc0 Second derivative parameters
		f = &(uni->Dc0_2ndDerCoeffs[0]);
		f[0] = - OmegaK;
		f[1] = - 1.5 * OmegaM;
		f[2] = -2.0 * OmegaR;

		//Dc0 Third derivative parameters
		f = uni->Dc0_3rdDerCoeffs;
		f[0] = - OmegaK * OmegaL;
		f[1] = -3.0 * OmegaL * OmegaM;
		f[2] = 2.0 * OmegaK * OmegaK - 6.0 * OmegaL * OmegaR;
		f[3] = 5.0 * OmegaK * OmegaM;
		f[4] = 5.0 * OmegaK * OmegaR + 3.75 * OmegaM * OmegaM;
		f[5] = 9.0 * OmegaM * OmegaR;
		f[6] = 6.0 * OmegaR * OmegaR;

		//tL0 Second derivative parameters
		f = uni->tL0_2ndDerCoeffs;
		f[0] = -OmegaL;
		f[1] = 0.0;
		f[2] = -2.0 * OmegaK;
		f[3] = -2.5 * OmegaM;
		f[4] = -3.0 * OmegaR;

		//tL0 Third derivative parameters
		f = uni->tL0_3rdDerCoeffs;
		f[0] = 2.0 * OmegaL * OmegaL;
		f[1] = 0.0;
		f[2] = 5.0 * OmegaK * OmegaL;
		f[3] = 4.0 * OmegaL * OmegaM;
		f[4] = 2.0 * OmegaL * OmegaR + 6.0 * OmegaK * OmegaK;
		f[5] = 14.0 * OmegaK * OmegaM;
		f[6] = 15.0 * OmegaK * OmegaR + 8.75 * OmegaM * OmegaM;
		f[7] = 20.0 * OmegaM * OmegaR;
		f[8] = 12.0 * OmegaR * OmegaR;

		//Dc0 Low redshift Taylor expansion parameters
		f = uni->Dc0_TaylorCoeffs;
		f[0] = 1.0;

		f[1] = (- OmegaR \
				- 0.5 * OmegaM \
				+ OmegaL \
				- 1.0) * 0.5;

		f[2] = (1.5 * OmegaR * OmegaR \
				+ 1.5 * OmegaR * OmegaM \
				- 3.0 * OmegaL * OmegaR \
				+ 0.5 * OmegaR \
				+ 0.375 * OmegaM * OmegaM \
				- 1.5 * OmegaL * OmegaM \
				+ 0.5 * OmegaM \
				+ 1.5 * OmegaL * OmegaL \
				- 2.5 * OmegaL \
				+ 1.0 ) / 3.0;

		f[3] = ( - 2.5 * OmegaR * OmegaR * OmegaR \
				- 3.75 * OmegaR * OmegaR * OmegaM \
				+ 7.5 * OmegaR * OmegaR * OmegaL \
				- 1.875 * OmegaR * OmegaM * OmegaM \
				+ 7.5 * OmegaR * OmegaM * OmegaL \
				- 0.75 * OmegaR * OmegaM \
				- 7.5 * OmegaR * OmegaL * OmegaL \
				+ 6.0 * OmegaR * OmegaL \
				- 0.5 * OmegaR \
				- 0.3125 * OmegaM * OmegaM * OmegaM \
				+ 1.875 * OmegaM * OmegaM * OmegaL \
				- 0.375 * OmegaM * OmegaM \
				- 3.75 * OmegaM * OmegaL * OmegaL \
				+ 3.75 * OmegaM * OmegaL \
				- 0.5 * OmegaM \
				+ 2.5 * OmegaL * OmegaL * OmegaL \
				- 6.0 * OmegaL * OmegaL \
				+ 4.5 * OmegaL \
				- 1.0 ) * 0.25;

		f[4] = ( 4.375 * OmegaR * OmegaR * OmegaR * OmegaR \
				+ 8.75 * OmegaM * OmegaR * OmegaR * OmegaR \
				- 17.5 * OmegaL * OmegaR * OmegaR * OmegaR \
				- 1.25 * OmegaR * OmegaR * OmegaR \
				+ 6.5625 * OmegaM * OmegaM * OmegaR * OmegaR \
				- 26.25 * OmegaM * OmegaL * OmegaR * OmegaR \
				+ 26.25 * OmegaL * OmegaL * OmegaR * OmegaR \
				- 11.25 * OmegaL * OmegaR * OmegaR \
				+ 0.375 * OmegaR * OmegaR \
				+ 2.1875 * OmegaM * OmegaM * OmegaM * OmegaR \
				- 13.125 * OmegaL * OmegaM * OmegaM * OmegaR \
				+ 0.9375 * OmegaM * OmegaM * OmegaR \
				+ 26.25 * OmegaL * OmegaL * OmegaM * OmegaR \
				- 15.0 * OmegaL * OmegaM * OmegaR \
				+ 0.75 * OmegaM * OmegaR \
				- 17.5 * OmegaL * OmegaL * OmegaL * OmegaR \
				+ 26.25 * OmegaL * OmegaL * OmegaR \
				- 9.75 * OmegaL * OmegaR \
				+ 0.5 * OmegaR \
				+ 0.2734375 * OmegaM * OmegaM * OmegaM * OmegaM \
				- 2.1875 * OmegaL * OmegaM * OmegaM * OmegaM \
				+ 0.3125 * OmegaM * OmegaM * OmegaM \
				+ 6.5625 * OmegaL * OmegaL * OmegaM * OmegaM \
				- 4.6875 * OmegaL * OmegaM * OmegaM \
				+ 0.375 * OmegaM * OmegaM \
				- 8.75 * OmegaL * OmegaL * OmegaL * OmegaM \
				+ 15.0 * OmegaL * OmegaL * OmegaM \
				- 6.75 * OmegaL * OmegaM \
				+ 0.5 * OmegaM \
				+ 4.375 * OmegaL * OmegaL * OmegaL * OmegaL \
				- 13.75 * OmegaL * OmegaL * OmegaL \
				+ 15.375 * OmegaL * OmegaL \
				- 7.0 * OmegaL \
				+ 1.0 ) / 5.0;

		//tL0 Low redshift Taylor expansion parameters
		f = uni->tL0_TaylorCoeffs;
		f[0] = 1.0;

		f[1] = ( - OmegaR \
				- 0.5 * OmegaM \
				+ OmegaL \
				- 2.0 ) * 0.5;

		f[2] = ( 0.5 * OmegaR * OmegaR \
				+ 0.5 * OmegaM * OmegaR \
				- OmegaL * OmegaR \
				+ 0.5 * OmegaR \
				+ 0.125 * OmegaM * OmegaM \
				- 0.5 * OmegaL * OmegaM \
				+ OmegaM / 3.0 \
				+ 0.5 * OmegaL * OmegaL \
				- 7.0 * OmegaL / 6.0 \
				+ 1.0 );

		f[3] = ( - 2.5 * OmegaR * OmegaR * OmegaR \
				- 3.75 * OmegaM * OmegaR * OmegaR \
				+ 7.5 * OmegaL * OmegaR * OmegaR \
				- 1.5 * OmegaR * OmegaR \
				- 1.875 * OmegaM * OmegaM * OmegaR \
				+ 7.5 * OmegaL * OmegaM * OmegaR \
				+ 2.25 * OmegaM * OmegaR \
				- 7.5 * OmegaL * OmegaL * OmegaR \
				+ 9.0 * OmegaL * OmegaR \
				- 2.0 * OmegaR \
				- 0.3125 * OmegaM * OmegaM * OmegaM \
				+ 1.875 * OmegaL * OmegaM * OmegaM \
				- 0.75 * OmegaM * OmegaM \
				- 3.75 * OmegaL * OmegaL * OmegaM \
				+ 5.25 * OmegaL * OmegaM \
				- 1.5 * OmegaM \
				+ 2.5 * OmegaL * OmegaL * OmegaL \
				- 7.5 * OmegaL * OmegaL \
				+ 8.0 * OmegaL \
				- 4.0 ) * 0.25;

		f[4] = ( 4.375 * OmegaR * OmegaR * OmegaR * OmegaR \
				+ 8.75 * OmegaM * OmegaR * OmegaR * OmegaR \
				- 17.5 * OmegaL * OmegaR * OmegaR * OmegaR \
				+ 1.25 * OmegaR * OmegaR * OmegaR \
				+ 6.5625 * OmegaM * OmegaM * OmegaR * OmegaR \
				- 26.25 * OmegaL * OmegaM * OmegaR * OmegaR \
				+ 3.75 * OmegaM * OmegaR * OmegaR \
				+ 26.25 * OmegaL * OmegaL * OmegaR * OmegaR \
				- 18.75 * OmegaL * OmegaR * OmegaR \
				+ 1.875 * OmegaR * OmegaR \
				+ 2.1875 * OmegaM * OmegaM * OmegaM * OmegaR \
				- 13.125 * OmegaL * OmegaM * OmegaM * OmegaR \
				+ 2.8125 * OmegaM * OmegaM * OmegaR \
				+ 26.25 * OmegaL * OmegaL * OmegaM * OmegaR \
				- 22.5 * OmegaL * OmegaM * OmegaR \
				+ 3.0 * OmegaM * OmegaR \
				- 17.5 * OmegaL * OmegaL * OmegaL * OmegaR \
				+ 33.75 * OmegaL * OmegaL * OmegaR \
				- 18.75 * OmegaL * OmegaR \
				+ 2.5 * OmegaR \
				+ 0.2734375 * OmegaM * OmegaM * OmegaM * OmegaM \
				- 2.1875 * OmegaL * OmegaM * OmegaM * OmegaM \
				+ 0.625 * OmegaM * OmegaM * OmegaM \
				+ 6.5625 * OmegaL * OmegaL * OmegaM * OmegaM \
				- 6.5625 * OmegaL * OmegaM * OmegaM \
				+ 1.125 * OmegaM * OmegaM \
				- 8.75 * OmegaL * OmegaL * OmegaL * OmegaM \
				+ 18.75 * OmegaL * OmegaL * OmegaM \
				- 12.0 * OmegaL * OmegaM \
				+ 2.0 * OmegaM \
				+ 4.375 * OmegaL * OmegaL * OmegaL * OmegaL \
				- 16.25 * OmegaL * OmegaL * OmegaL \
				+ 22.875 * OmegaL * OmegaL \
				- 15.0 * OmegaL \
				+ 5.0 ) / 5.0;
	}

	//Cache only works if there are no turning points in the cache interval
	if (uni->zLastTurn > zMaxInterp && uni->zLastTurn > zMinInterp && \
		uni->zNextTurn < zMaxInterp && uni->zNextTurn < zMinInterp ) {

		double errest;
		double cur_Dc0, cur_tL0;
		CachePoint *curPt_Dc0, *curPt_tL0;
		double curAinv = 1.0 + zMinInterp;
		int idx;
		int DcintegratorStatus, tLintegratorStatus;
		gsl_integration_workspace *scratch;
		gsl_function Dintegrand, tLintegrand;

		//Integrate to the left side of the cache using GSL
		scratch = gsl_integration_workspace_alloc( MaxGSLlevel );

		Dintegrand.function = &dDc0_dz_gsl;
		Dintegrand.params = (void *)uni;
		DcintegratorStatus = gsl_integration_qag( &Dintegrand, \
			1.0, curAinv, \
			GSLCacheEpsabs, GSLCacheEpsrel, MaxGSLlevel, GSL_INTEG_GAUSS51, \
			scratch, &cur_Dc0, &errest);

		tLintegrand.function = &dtL0_dz_gsl;
		tLintegrand.params = (void *)uni;
		tLintegratorStatus = gsl_integration_qag( &tLintegrand, \
			1.0, curAinv, \
			GSLCacheEpsabs, GSLCacheEpsrel, MaxGSLlevel, GSL_INTEG_GAUSS51, \
			scratch, &cur_tL0, &errest);

		if (DcintegratorStatus == GSL_EFAILED || tLintegratorStatus==GSL_EFAILED){
			uni->CacheValid = false;
		}

		else if (GSLcache) {
			double Dc0ders[3], tL0ders[3];

			uni->CacheValid = true;
			curPt_Dc0 = uni->Dc0Cache;
			curPt_tL0 = uni->tL0Cache;
			curPt_Dc0->ders[0] = cur_Dc0;
			curPt_tL0->ders[0] = cur_tL0;
			for (idx = 1; idx < InterpIntervals + 1; ++idx) {
				++curPt_Dc0;
				++curPt_tL0;
				curAinv = 1.0 + zMinInterp + zInterpCellWidth * (double) idx;

				DcintegratorStatus = gsl_integration_qag( &Dintegrand, \
				1.0, curAinv, \
				GSLCacheEpsabs, GSLCacheEpsrel, MaxGSLlevel, GSL_INTEG_GAUSS51, \
				scratch, &cur_Dc0, &errest);

				tLintegratorStatus = gsl_integration_qag( &tLintegrand, \
				1.0, curAinv, \
				GSLCacheEpsabs, GSLCacheEpsrel, MaxGSLlevel, GSL_INTEG_GAUSS51, \
				scratch, &cur_tL0, &errest);

				if (DcintegratorStatus == GSL_EFAILED || \
					tLintegratorStatus == GSL_EFAILED) {
					uni->CacheValid = false;
					break;
				}

				curPt_Dc0->ders[0] = cur_Dc0;
				curPt_tL0->ders[0] = cur_tL0;

				Dc0Derivatives(uni, curAinv, Dc0ders);
				tL0Derivatives(uni, curAinv, tL0ders);
				memcpy( &(curPt_Dc0->ders[1]), Dc0ders, 3 * sizeof(double) );
				memcpy( &(curPt_tL0->ders[1]), tL0ders, 3 * sizeof(double) );
			}
		}

		else {
			double deltaDc0, deltatL0;
			int intervalsPerCell = 4;
			double deltazinteg = zInterpCellWidth / (double) (intervalsPerCell << 1);
			double deltazinteg2 = 2.0 * deltazinteg;
			double Diffs_Dc0[9], Diffs_tL0[9];
			double *leftDiffs_Dc0, *centDiffs_Dc0, *rightDiffs_Dc0;
			double *leftDiffs_tL0, *centDiffs_tL0, *rightDiffs_tL0;
			double *temp;
			int cacheInt = 0;

			//Do a rolling integral using a 3 point 0-3 derivative rule to accumulate
			//the cache.
			leftDiffs_Dc0 = Diffs_Dc0;
			centDiffs_Dc0 = leftDiffs_Dc0 + 3;
			rightDiffs_Dc0 = centDiffs_Dc0 + 3;
			leftDiffs_tL0 = Diffs_tL0;
			centDiffs_tL0 = leftDiffs_tL0 + 3;
			rightDiffs_tL0 = centDiffs_tL0 + 3;

			Dc0Derivatives(uni, curAinv, rightDiffs_Dc0);
			tL0Derivatives(uni, curAinv, rightDiffs_tL0);

			curPt_Dc0 = uni->Dc0Cache;
			curPt_tL0 = uni->tL0Cache;
			curPt_Dc0->ders[0] = cur_Dc0;
			curPt_tL0->ders[0] = cur_tL0;
			memcpy( &(curPt_Dc0->ders[1]), rightDiffs_Dc0, 3 * sizeof(double) );
			memcpy( &(curPt_tL0->ders[1]), rightDiffs_tL0, 3 * sizeof(double) );

			deltaDc0 = 0.0;
			deltatL0 = 0.0;
			for (idx = 0; idx < InterpIntervals * intervalsPerCell; ++idx) {
				//swap left and right diffs pointers
				temp = leftDiffs_Dc0;
				leftDiffs_Dc0 = rightDiffs_Dc0;
				rightDiffs_Dc0 = temp;

				temp = leftDiffs_tL0;
				leftDiffs_tL0 = rightDiffs_tL0;
				rightDiffs_tL0 = temp;

				//Fill in the center and right Diffs pointers
				curAinv += deltazinteg;
				Dc0Derivatives(uni, curAinv, centDiffs_Dc0);
				tL0Derivatives(uni, curAinv, centDiffs_tL0);

				curAinv += deltazinteg;
				Dc0Derivatives(uni, curAinv, rightDiffs_Dc0);
				tL0Derivatives(uni, curAinv, rightDiffs_tL0);

				deltaDc0 = FMA_m(deltazinteg2, \
					MeanValueInterval(deltazinteg, \
						leftDiffs_Dc0, centDiffs_Dc0, rightDiffs_Dc0), deltaDc0);
				deltatL0 = FMA_m(deltazinteg2, \
					MeanValueInterval(deltazinteg, \
						leftDiffs_tL0, centDiffs_tL0, rightDiffs_tL0), deltatL0);

				cacheInt += 1;
				if (cacheInt == intervalsPerCell) {
					//Cache the right hand side.
					cacheInt = 0;

					cur_Dc0 += deltaDc0;
					cur_tL0 += deltatL0;
					deltaDc0 = 0.0;
					deltatL0 = 0.0;

					++curPt_Dc0;
					++curPt_tL0;

					curPt_Dc0->ders[0] = cur_Dc0;
					curPt_tL0->ders[0] = cur_tL0;
					memcpy( &(curPt_Dc0->ders[1]), rightDiffs_Dc0, 3 * sizeof(double) );
					memcpy( &(curPt_tL0->ders[1]), rightDiffs_tL0, 3 * sizeof(double) );
				}
			}

			uni->CacheValid = true;
		}

		gsl_integration_workspace_free( scratch );
	}

	else {
		fprintf(stderr, "Warning: InitUniverse cannot produce an integration cache\n%s",
			"if the cosmology has a turning point in the cache interval.\n");
		uni->CacheValid = false;
	}

	{ //Assign function pointers
		if (uni->flat || fabs(OmegaK) <= NumericallyFlat) {
			uni->F_DcT0 = &DcT0Flat;
			uni->F_Vc0	= &Vc0Flat;
			uni->F_DcT0Inv = &DcT0InvFlat;
			uni->F_Da0Inv = &Da0InvFlat;
			uni->F_Dl0Inv = &Dl0InvFlat;
			uni->F_Vc0Inv = &Vc0InvFlat;
		}

		else if (OmegaK < 0.0) {
			uni->F_DcT0 = &DcT0NegK;
			uni->F_Vc0 = &Vc0NegK;
			uni->F_DcT0Inv = &DcT0InvNegK;
			uni->F_Da0Inv = &Da0InvNegK;
			uni->F_Dl0Inv = &Dl0InvNegK;
			uni->F_Vc0Inv = &Vc0InvNegK;
		}

		else {
			uni->F_DcT0 = &DcT0PosK;
			uni->F_Vc0 = &Vc0PosK;
			uni->F_DcT0Inv = &DcT0InvPosK;
			uni->F_Da0Inv = &Da0InvPosK;
			uni->F_Dl0Inv = &Dl0InvPosK;
			uni->F_Vc0Inv = &Vc0InvPosK;
		}

	}

	{//Calculate vital statistics
		double dummy;
		if (isinf(uni->zLastTurn)) {
			uni->age0 = tL0_Cosmic(uni, INFINITY, true,
				1e-6, 1e-6, MaxGSLlevel, &dummy);
			uni->Dc0Max = Dc0_Cosmic(uni, INFINITY, true,
				1e-6, 1e-6, MaxGSLlevel, &dummy);

			uni->tL0LastTurn = NAN;
			uni->Dc0LastTurn = NAN;
			uni->Dc0InvD0max = uni->Dc0Max;
		}
		else {
			uni->age0 = INFINITY;
			uni->Dc0Max = INFINITY;

			uni->tL0LastTurn = tL0_Cosmic(uni, uni->zLastTurn, true,
				1e-6, 1e-6, MaxGSLlevel, &dummy);
			uni->Dc0LastTurn = Dc0_Cosmic(uni, uni->zLastTurn, true,
				1e-6, 1e-6, MaxGSLlevel, &dummy);
			uni->Dc0InvD0max = uni->Dc0LastTurn;
		}

		if (uni->zNextTurn <= -1.0) {
			uni->tL0NextTurn = INFINITY;
		}
		else {
			uni->tL0NextTurn = tL0_Cosmic(uni, uni->zNextTurn, true,
				1e-6, 1e-6, MaxGSLlevel, &dummy);
		}

		uni->z_DaMax = FindDa0Max(uni, 0x1.0p-20, 1024);

		if (uni->flat || fabs(OmegaK) <= NumericallyFlat) {
			uni->DcT0Max = uni->Dc0Max;

			uni->Vc0Universe = Vc0Flat(uni, uni->zLastTurn, true, \
				1e-6, 1e-6, MaxGSLlevel, &dummy );
			uni->Da0Max = Dc0_Cosmic(uni, uni->z_DaMax, true, \
				1e-6, 1e-6, 128, &dummy) / (1.0 + uni->z_DaMax);
		}
		else if (OmegaK < 0.0) {
			double N = sqrt(-OmegaK);
			uni->DcT0Max = uni->DH * N / (-OmegaK);
			uni->DcT0Max = MIN( uni->DcT0Max,
				sin( N * uni->DcT0Max ) * N / (-OmegaK) );

			uni->Vc0Universe = Vc0NegK(uni, uni->zLastTurn, true, \
				1e-6, 1e-6, MaxGSLlevel, &dummy );
			uni->Da0Max = DcT0NegK(uni, uni->z_DaMax, true, \
				1e-6, 1e-6, 128, &dummy) / (1.0 + uni->z_DaMax);
		}
		else { //OmegaK > 0.0
			double N = sqrt(OmegaK);
			uni->DcT0Max = sinh( N * uni->Dc0Max ) * N / OmegaK;

			uni->Vc0Universe = Vc0PosK(uni, uni->zLastTurn, true, \
				1e-6, 1e-6, MaxGSLlevel, &dummy );
			uni->Da0Max = DcT0PosK(uni, uni->z_DaMax, true, \
				1e-6, 1e-6, 128, &dummy) / (1.0 + uni->z_DaMax);
		}

	}

	return true;
}

UniverseLCDM *MakeUniverse(double H0, double OmegaL,
					double OmegaM, double OmegaR, bool flat, bool GSLcache) {
	UniverseLCDM *uni;
	bool unigood;

	if ( H0 < 0.0 || OmegaM < 0.0 || OmegaR < 0.0 ) {
		uni = NULL;
	}
	else {
		uni = (UniverseLCDM *) malloc( sizeof (UniverseLCDM) );
		unigood = InitUniverse(uni, H0, OmegaL, OmegaM, OmegaR, flat,\
								GSLcache);

		if (unigood == false) {
			free(uni);
			uni = NULL;
		}
	}

	return uni;
}

void FreeUniverse(UniverseLCDM *uni) {
	/*if (uni != NULL) {
		free(uni);
	}*/
	return;
}

double Dc0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double result, weight, zpos;

	if (z < -1.0) {
		fprintf(stderr, "Dc0_Cosmic: invalid redshift. Redshifts must be >= -1.0\n");
		fprintf(stderr, "z = %.16e\n", z);

		return NAN;
	}
	else if (z > uni->zLastTurn || z < uni->zNextTurn) {
		fprintf(stderr, "Dc0_Cosmic: invalid redshift. Cannot integrate past turning points.\n");
		fprintf(stderr, "z = %.16e, zNextTurn = %.16e, zLastTurn = %.16e\n",
			z, uni->zNextTurn, uni->zLastTurn );

		return NAN;
	}

	zpos = fabs(z);
	result = 0.0;
	weight = 1.0;

	if (zpos <= TaylorMaxZ) { //Perform the taylor series calculation
		double *coeffs = uni->Dc0_TaylorCoeffs;

		if (zpos > TaylorInterpMinZ) {
			weight = (TaylorMaxZ - zpos) * TaylorInterpInterval;
		}
		else {
			weight = 1.0;
		}

		result = FMA_m(FMA_m(coeffs[4], z, coeffs[3]), z, coeffs[2]);
		result = FMA_m(FMA_m(result, z, coeffs[1]), z, coeffs[0]) * z * weight;

		weight = 1.0 - weight;
	}

	*abserr = 0.0;
	if (fast && uni->CacheValid && z >= zMinInterp && z <= zMaxInterp \
		&& zpos > TaylorInterpMinZ) {
		double tempres, zCenter, x;
		double bases[4], *ders;
		int centeridx;
		CachePoint *curPt;

		//Use the interpolation cache
		centeridx = round((z - zMinInterp) / zInterpCellWidth);
		centeridx = CLAMP(centeridx, 1, InterpIntervals);

		zCenter = ((double) centeridx) * zInterpCellWidth + zMinInterp;
		x = (z - zCenter) / zInterpCellWidth;
		tempres = 0.0;

		curPt = uni->Dc0Cache + (centeridx - 1);
		ders = curPt->ders;
		LtEdgeBases(x, zInterpCellWidth, bases);
		tempres = FMA_m(bases[2], ders[2], FMA_m(bases[3], ders[3], tempres));
		tempres = FMA_m(bases[0], ders[0], FMA_m(bases[1], ders[1], tempres));

		++curPt;
		ders = curPt->ders;
		centerBases(x, zInterpCellWidth, bases);
		tempres = FMA_m(bases[2], ders[2], FMA_m(bases[3], ders[3], tempres));
		tempres = FMA_m(bases[0], ders[0], FMA_m(bases[1], ders[1], tempres));

		++curPt;
		ders = curPt->ders;
		RtEdgeBases(x, zInterpCellWidth, bases);
		tempres = FMA_m(bases[2], ders[2], FMA_m(bases[3], ders[3], tempres));
		tempres = FMA_m(bases[0], ders[0], FMA_m(bases[1], ders[1], tempres));

		result = FMA_m(weight, tempres, result);
		*abserr = tempres * pow(x, 12.0);
	}

	else if ( ((zpos > TaylorInterpMinZ && (!fast || !uni->CacheValid)) ||
		z > zMaxInterp || z < zMinInterp) ) {
		gsl_integration_workspace *scratch;
		gsl_function integrand;
		int DcintegratorStatus;
		double tempres, temperr;

		scratch = gsl_integration_workspace_alloc( limit );

		integrand.function = &dDc0_dz_gsl;
		integrand.params = (void *)uni;

		if (!isinf(z) && \
			(z >= uni->zLastTurn * 0x1.fffffff8p-1 || \
				z <= uni->zNextTurn * 0x1.fffffff8p-1 || \
				z <= -1.0 * 0x1.fffffff8p-1)) {
			DcintegratorStatus = gsl_integration_qags(&integrand, \
				1.0, 1.0 + z, \
				epsabs, epsrel, limit, \
				scratch, &tempres, &temperr);
		}

		else if (isinf(z)) {
			DcintegratorStatus = gsl_integration_qagiu(&integrand, \
				1.0, \
				epsabs, epsrel, limit, \
				scratch, &tempres, &temperr);
		}

		else {
			double z0 = 0.0;
			double Dc0_z0 = 0.0;

			if ( uni->CacheValid ) {
				if ( z > zMaxInterp ) {
					z0 = zMaxInterp;
					Dc0_z0 = uni->Dc0Cache[InterpIntervals - 1].ders[0];
				}
				else if ( z < zMinInterp ) {
					z0 = zMinInterp;
					Dc0_z0 = uni->Dc0Cache[0].ders[0];
				}
			}

			DcintegratorStatus = gsl_integration_qag(&integrand, \
				1.0 + z0, 1.0 + z, \
				epsabs, epsrel, limit, GSL_INTEG_GAUSS41, \
				scratch, &tempres, &temperr);

			tempres += Dc0_z0;
		}

		gsl_integration_workspace_free( scratch );

		if (DcintegratorStatus == GSL_EFAILED) {
			fprintf(stderr, "Dc0_Cosmic: gsl integrator failed.\n");
			*abserr = NAN;
			return NAN;
		}

		result = FMA_m(weight, tempres, result);
		*abserr = FMA_m(weight, temperr, *abserr);
	}

	return result;
}

double Dc_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	return Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr) * uni->DH;
}

double tL0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double result, weight, zpos;

	if (z < -1.0) {
		fprintf(stderr, "tL0_Cosmic: invalid redshift. Redshifts must be >= -1.0\n");
		fprintf(stderr, "z = %.16e\n", z);

		return NAN;
	}
	else if (z > uni->zLastTurn || z < uni->zNextTurn) {
		fprintf(stderr, "tL0_Cosmic: invalid redshift. Cannot integrate past turning points.\n");
		fprintf(stderr, "z = %.16e, zNextTurn = %.16e, zLastTurn = %.16e\n",
			z, uni->zNextTurn, uni->zLastTurn );

		return NAN;
	}

	zpos = fabs(z);
	result = 0.0;
	weight = 1.0;

	if (zpos <= TaylorMaxZ) { //Perform the taylor series calculation
		double *coeffs = uni->tL0_TaylorCoeffs;

		if (zpos > TaylorInterpMinZ) {
			weight = (TaylorMaxZ - zpos) * TaylorInterpInterval;
		}
		else {
			weight = 1.0;
		}

		result = FMA_m(FMA_m(coeffs[4], z, coeffs[3]), z, coeffs[2]);
		result = FMA_m(FMA_m(result, z, coeffs[1]), z, coeffs[0]) * z * weight;

		weight = 1.0 - weight;
	}

	*abserr = 0.0;

	if (fast && uni->CacheValid && z >= zMinInterp && z <= zMaxInterp \
		&& zpos > TaylorInterpMinZ) {
		double tempres, zCenter, x;
		double bases[4], *ders;
		int centeridx;
		CachePoint *curPt;

		//Use the interpolation cache
		centeridx = round((z - zMinInterp) / zInterpCellWidth);
		if (centeridx < 1) {
			centeridx = 1;
		}

		else if (centeridx >= InterpIntervals + 1) {
			centeridx = InterpIntervals;
		}

		zCenter = ((double) centeridx) * zInterpCellWidth + zMinInterp;
		x = (z - zCenter) / zInterpCellWidth;
		tempres = 0.0;

		curPt = uni->tL0Cache + centeridx - 1;
		ders = curPt->ders;
		LtEdgeBases(x, zInterpCellWidth, bases);
		tempres = FMA_m(bases[2], ders[2], FMA_m(bases[3], ders[3], tempres));
		tempres = FMA_m(bases[0], ders[0], FMA_m(bases[1], ders[1], tempres));

		++curPt;
		ders = curPt->ders;
		centerBases(x, zInterpCellWidth, bases);
		tempres = FMA_m(bases[2], ders[2], FMA_m(bases[3], ders[3], tempres));
		tempres = FMA_m(bases[0], ders[0], FMA_m(bases[1], ders[1], tempres));

		++curPt;
		ders = curPt->ders;
		RtEdgeBases(x, zInterpCellWidth, bases);
		tempres = FMA_m(bases[2], ders[2], FMA_m(bases[3], ders[3], tempres));
		tempres = FMA_m(bases[0], ders[0], FMA_m(bases[1], ders[1], tempres));

		result = FMA_m(weight, tempres, result);
		*abserr = tempres * pow(x, 12.0);
	}

	else if ( ((zpos > TaylorInterpMinZ && !fast) ||
		z > zMaxInterp || z < zMinInterp) ) {
		gsl_integration_workspace *scratch;
		gsl_function integrand;
		int tLintegratorStatus;
		double tempres, temperr;

		scratch = gsl_integration_workspace_alloc( limit );

		integrand.function = &dtL0_dz_gsl;
		integrand.params = (void *)uni;

		if (!isinf(z) && \
			(z >= uni->zLastTurn * 0x1.fffffff8p-1 || \
				z <= uni->zNextTurn * 0x1.fffffff8p-1 || \
				z <= -1.0 * 0x1.fffffff8p-1)) {
			tLintegratorStatus = gsl_integration_qags(&integrand, \
				1.0, 1.0 + z, \
				epsabs, epsrel, limit, \
				scratch, &tempres, &temperr);
		}

		else if (isinf(z)) {
			tLintegratorStatus = gsl_integration_qagiu(&integrand, \
				1.0, \
				epsabs, epsrel, limit, \
				scratch, &tempres, &temperr);
		}

		else {
			double z0 = 0.0;
			double tL0_z0 = 0.0;

			if ( uni->CacheValid ) {
				if ( z > zMaxInterp ) {
					z0 = zMaxInterp;
					tL0_z0 = uni->tL0Cache[InterpIntervals - 1].ders[0];
				}
				else if ( z < zMinInterp ) {
					z0 = zMinInterp;
					tL0_z0 = uni->tL0Cache[0].ders[0];
				}
			}

			tLintegratorStatus = gsl_integration_qag(&integrand, \
				1.0 + z0, 1.0 + z, \
				epsabs, epsrel, limit, GSL_INTEG_GAUSS41, \
				scratch, &tempres, &temperr);

			tempres += tL0_z0;
		}

		gsl_integration_workspace_free( scratch );

		if (tLintegratorStatus == GSL_EFAILED) {
			fprintf(stderr, "tL0_Cosmic: gsl integrator failed.\n");
			*abserr = NAN;
			return NAN;
		}

		result = FMA_m(weight, tempres, result);
		*abserr = FMA_m(weight, temperr, *abserr);
	}

	return result;
}

double tL_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	return tL0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr) * uni->tH;
}

double DcT_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double result;

	result = (*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr);

	result *= uni->DH;
	*abserr *= uni->DH;

	return result;
}

double DcT0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return (*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr);
}

static inline double DcT0Flat(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr);
}

static double DcT0NegK(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double result, nrm;
	double OmegaKpos = -uni->OmegaCurvature;

	result = Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr);
	nrm = sqrt(OmegaKpos);
	result = sin(nrm * result) * nrm / OmegaKpos;

	return result;
}

static double DcT0PosK(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double result, nrm;
	double OmegaKpos = uni->OmegaCurvature;

	result = Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr);
	nrm = sqrt(uni->OmegaCurvature);

	result = sinh(nrm * result) * nrm / uni->OmegaCurvature;

	return result;
}

double Da0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return ((*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr) \
		/ (1.0 + z));
}

double Da_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return ((*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr) \
		* uni->DH / (1.0 + z));
}

double Dl0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return ((*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr) \
		* (1.0 + z));
}

double Dl_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return ((*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr) \
		* (1.0 + z) * uni->DH);
}

double dVc_dzdOmega0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double D = (*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr);

	return D * D * Einv(uni, 1.0 + z);
}

double dVc_dzdOmega_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double D = (*(uni->F_DcT0))(uni, z, fast, epsabs, epsrel, limit, abserr);

	return D * D * Einv(uni, 1.0 + z) * uni->VH;
}

double Vc_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return ((*(uni->F_Vc0))(uni, z, fast, epsabs, epsrel, limit, abserr) \
		* uni->VH);
}

double Vc0_Cosmic(UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {

	return (*(uni->F_Vc0))(uni, z, fast, epsabs, epsrel, limit, abserr);
}

static double Vc0Flat(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double DcT = Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr);
	double result, scale;

	scale = 4.0*Pi * DcT * DcT;
	result = DcT * scale / 3.0;
	*abserr *= scale;

	return result;
}

static double Vc0NegK(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double Dc = Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr);
	double OmegaKpos = -uni->OmegaCurvature;
	double x, nrm, result, scale;

	nrm = sqrt(OmegaKpos);
	x = 2.0 * Dc * nrm;

	if ( fabs(x) < 0x1.cp-2 ) { //Use Taylor expansion where more accurate
		double x2 = x*x;
		double errscale;

		//Taylor expansion for sin(x) - x
		result = FMA_m(0x1.6124613a86d09p-33, x2, -0x1.ae64567f544e4p-26);
		result = FMA_m(FMA_m(result, x2, 0x1.71de3a556c734p-19), x2, -0x1.a01a01a01a01ap-13);
		result = FMA_m(FMA_m(result, x2, 0x1.1111111111111p-7), x2, -0x1.5555555555555p-3);
		result *= x2 * x;
	}

	else {
		result = sin(x) - x;
	}

	//First scaling of the abserr
	scale = sin(0.5 * x);
	scale *= 4.0 * Pi * scale / OmegaKpos;
	*abserr *= scale;

	//Final scaling of results
	result *= -Pi * nrm / (OmegaKpos * OmegaKpos);

	return result;
}

static double Vc0PosK(struct UniverseLCDM *uni, double z, bool fast,
	double epsabs, double epsrel, size_t limit, double *abserr) {
	double Dc = Dc0_Cosmic(uni, z, fast, epsabs, epsrel, limit, abserr);
	double OmegaKpos = uni->OmegaCurvature;
	double x, nrm, result, scale;

	nrm = sqrt(OmegaKpos);
	x = 2.0 * Dc * nrm;

	if ( fabs(x) < 0.5 ) { //Use Taylor expansion where more accurate
		double x2 = x*x;

		//Taylor expansion for sinh(x) - x
		result = FMA_m(0x1.6124613a86d09p-33, x2, 0x1.ae64567f544e4p-26);
		result = FMA_m(FMA_m(result, x2, 0x1.71de3a556c734p-19), x2, 0x1.a01a01a01a01ap-13);
		result = FMA_m(FMA_m(result, x2, 0x1.1111111111111p-7), x2, 0x1.5555555555555p-3);
		result *= x2 * x;
	}

	else {
		result = sinh(x) - x;
	}

	//First scaling of the abserr
	scale = sinh(0.5 * x);
	scale *= 4.0 * Pi * scale / OmegaKpos;
	*abserr *= scale;

	//Final scaling of results
	result *= Pi * nrm / (OmegaKpos * OmegaKpos);

	return result;
}

double Dc0Inv_Cosmic(UniverseLCDM *uni, double Dc0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, f, fp, fppoverfp, foverfp, delta, ainv, dummy;
	double c3, c2, c1;
	double epsrelDc0 = 0.0625 * epsrel;
	double epsabsDc0 = epsrelDc0 * Dc0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Dc0) || isnan(epsrel)) {
		return NAN;
	} else if (Dc0 < uni->Dc0InvD0max) {
		//do nothing
	}
	else if (!isnan(uni->Dc0LastTurn) && Dc0 == uni->Dc0LastTurn) {
		return uni->zLastTurn;
	}
	else if (Dc0 == uni->Dc0Max) {
		return INFINITY;
	}
	else {
		fprintf(stderr, "Dc0Inv_Cosmic: Dc0 outside the invertable region.\n");
		return NAN;
	}

	//Approximation for Dc0 = z / (1 + OmegaM * z )
	zCur = Dc0 / (1.0 - uni->OmegaMatter * Dc0);
	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Iterate using modified Halley's method
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;
		f = (Dc0_Cosmic(uni, zCur, true, epsabsDc0, epsrelDc0, MaxGSLlevel, &dummy) \
			- Dc0);
		ainv = 1.0 + zCur;
		fp = Einv(uni, ainv);

		fppoverfp = FMA_m(FMA_m(c3, ainv, c2), ainv, c1) * ainv * fp * fp;
		foverfp = f / fp;

		delta = -foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= -1.0) {
			if (zCur < -1.0) {
				zCur = fabs(zCur + 1.0) - 1.0;
			}
			else {
				zCur = -0.5;
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) <= epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double DcInv_Cosmic(UniverseLCDM *uni, double Dc,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double Dc0 = Dc / uni->DH;

	return Dc0Inv_Cosmic(uni, Dc0, branchNum, epsrel, maxiter);
}

double DcT0Inv_Cosmic(UniverseLCDM *uni, double DcT0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	return (*(uni->F_DcT0Inv))(uni, DcT0, branchNum, epsrel, maxiter);
}

static double DcT0InvFlat(struct UniverseLCDM *uni, double DcT0,
	unsigned int branchNum, double epsrel, unsigned int maxiter){
	return Dc0Inv_Cosmic(uni, DcT0, branchNum, epsrel, maxiter);
}

static double DcT0InvNegK(struct UniverseLCDM *uni, double DcT0, unsigned int branchNum,
	double epsrel, unsigned int maxiter){
	double x, nrm, Dc0;
	double OmegaKpos = -uni->OmegaCurvature;

	nrm = sqrt(OmegaKpos);
	x = asin(DcT0 * nrm);

	if (branchNum == 0) {
		//common case, do nothing
	}

	else if ((branchNum & 1) == 0) { //branchNum even
		x += copysign(2.0 * Pi * (double) branchNum, x);
	}

	else { //branchNum odd
		//Reflect across Pi/2
		x += 2.0 * copysign( 0.5 * Pi - fabs(x), x);

		//Add the remaining branch shifts
		x += copysign(2.0 * Pi * (double) (branchNum - 1), x);
	}

	Dc0 = x * nrm / OmegaKpos;

	return Dc0Inv_Cosmic(uni, Dc0, 0, epsrel, maxiter);
}

static double DcT0InvPosK(struct UniverseLCDM *uni, double DcT0, unsigned int branchNum,
	double epsrel, unsigned int maxiter){
	double x, nrm, Dc0;
	double OmegaKpos = uni->OmegaCurvature;

	nrm = sqrt(OmegaKpos);
	x = asinh(DcT0 * nrm);

	Dc0 = x * nrm / OmegaKpos;

	return Dc0Inv_Cosmic(uni, Dc0, 0, epsrel, maxiter);
}

double DcTInv_Cosmic(UniverseLCDM *uni, double DcT,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {

	return DcT0Inv_Cosmic(uni, DcT / uni->DH, branchNum, epsrel, maxiter);
}

double Da0Inv_Cosmic(UniverseLCDM *uni, double Da0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {

	return (*(uni->F_Da0Inv))(uni, Da0, branchNum, epsrel, maxiter);
}

double Da0InvFlat(struct UniverseLCDM *uni, double Da0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, f, ainv, a, Einv_, dummy;
	double fp, fppoverfp, foverfp, delta;
	double Dc0;
	double zMin, zMax;
	double c1, c2, c3;
	double epsrelDa0 = 0.0625 * epsrel;
	double epsabsDa0 = Da0 * epsrelDa0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Da0) || isnan(epsrel)) {
		return NAN;
	}
	else if (Da0 < uni->Da0Max) {
		if (branchNum == 0) {
			zMin = -1.0;
			zMax = uni->z_DaMax;
		}
		else {
			double z = uni->zLastTurn;
			double dummy;
			double Da0LastTurn = (DcT0Flat(uni, z, true, 1e-6, epsrel, 128, &dummy) \
				/ (1.0 + z));

			if (Da0 < Da0LastTurn) {
				return NAN;
			}
			else if (Da0 == Da0LastTurn) {
				return uni->zLastTurn;
			}

			zMin = uni->z_DaMax;
			zMax = uni->zLastTurn;
		}
	}
	else if (Da0 == uni->Da0Max) {
		return uni->z_DaMax;
	}
	else {
		fprintf(stderr, "Da0InvFlat: Da0 outside the invertable region.\n");
		printf("%.16e\n", uni->Da0Max);
		return NAN;
	}

	if (branchNum == 0) {
		zCur = 0.5 * (zMin + zMax);
	}
	else {
		zCur = 2.0 * zMin;
	}
	f = Da0;

	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Because Da often has a maximum, we proceed using Halley's iteration instead of Newtons.
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;

		ainv = 1.0 + zCur;
		a = 1.0 / ainv;
		Dc0 = Dc0_Cosmic(uni, zCur, true, epsabsDa0, epsrelDa0, 128, &dummy);
		Einv_ = Einv(uni, ainv);

		f = FMA_m(Dc0, a, - Da0);
		fp = FMA_m(-Dc0, a, Einv_) * a;
		foverfp = f / fp;

		fppoverfp = FMA_m(FMA_m(c3, ainv, c2), ainv, c1) * ainv / fp;
		fppoverfp = FMA_m(fppoverfp, Einv_ * Einv_ * Einv_, -2.0) * a;

		delta = - foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= zMin){
			if (zCur < zMin) {
				zCur = fabs(zCur - zMin) + zMin;
			}
			else {
				zCur = zMin + MIN(0.5, 0.125 * (zMax - zMin));
			}
			delta = zCur - zLast;
		}

		else if (zCur >= zMax) {
			if (zCur > zMax) {
				zCur = zMax - fabs(zMax - zCur);
			}
			else {
				zCur = zMax - MIN(0.5, 0.125 * (zMax - zMin));
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) <= epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double Da0InvNegK(struct UniverseLCDM *uni, double Da0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, f, ainv, a, Einv_, dummy;
	double fp, fppoverfp, foverfp, delta;
	double tanReduced, sinReduced, cosReduced, Dc0;
	double Wkpos = -uni->OmegaCurvature;
	double rootWk = sqrt(Wkpos);
	double zMin, zMax;
	double c1, c2, c3;
	double epsrelDa0 = 0.0625 * epsrel;
	double epsabsDa0 = Da0 * epsrelDa0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Da0) || isnan(epsrel)) {
		return NAN;
	}
	else if (Da0 < uni->Da0Max) {
		if (branchNum == 0) {
			zMin = -1.0;
			zMax = uni->z_DaMax;
		}
		else {
			double z = uni->zLastTurn;
			double dummy;
			double Da0LastTurn = (DcT0Flat(uni, z, true, 1e-6, epsrel, 128, &dummy) \
				/ (1.0 + z));

			if (Da0 < Da0LastTurn) {
				return NAN;
			}
			else if (Da0 == Da0LastTurn) {
				return uni->zLastTurn;
			}

			zMin = uni->z_DaMax;
			zMax = uni->zLastTurn;
		}
	}
	else if (Da0 == uni->Da0Max) {
		return uni->z_DaMax;
	}
	else {
		fprintf(stderr, "Da0Inv_Cosmic: Da0 outside the invertable region.\n");
		return NAN;
	}

	if (branchNum == 0) {
		zCur = 0.5 * (zMin + zMax);
	}
	else {
		zCur = 2.0 * zMin;
	}

	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Because Da often has a maximum, we proceed using Halley's iteration instead of Newtons.
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;

		ainv = 1.0 + zCur;
		a = 1.0 / ainv;
		Dc0 = Dc0_Cosmic(uni, zCur, true, epsabsDa0, epsrelDa0, 128, &dummy);
		Einv_ = Einv(uni, ainv);
		sinReduced = sin(Dc0 * rootWk)/rootWk;
		cosReduced = cos(Dc0 * rootWk);
		tanReduced = sinReduced / cosReduced;

		f =  FMA_m(sinReduced, a, -Da0);
		fp = FMA_m(-sinReduced, a, cosReduced * Einv_) * a;
		foverfp = f / fp;

		fppoverfp = FMA_m(FMA_m(c3, ainv, c2), ainv, c1) / fp;
		fppoverfp = FMA_m(fppoverfp, Einv_ * Einv_ * Einv_,
			FMA_m(Einv_ * Einv_, -Wkpos, 2.0 * a) * tanReduced);

		delta = - foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= zMin){
			if (zCur < zMin) {
				zCur = fabs(zCur - zMin) + zMin;
			}
			else {
				zCur = zMin + MIN(0.5, 0.125 * (zMax - zMin));
			}
			delta = zCur - zLast;
		}

		else if (zCur >= zMax) {
			if (zCur > zMax) {
				zCur = zMax - fabs(zMax - zCur);
			}
			else {
				zCur = zMax - MIN(0.5, 0.125 * (zMax - zMin));
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) <= epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double Da0InvPosK(struct UniverseLCDM *uni, double Da0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, f, ainv, a, Einv_, dummy;
	double fp, fppoverfp, foverfp, delta;
	double tanhReduced, sinhReduced, coshReduced, Dc0;
	double Wkpos = uni->OmegaCurvature;
	double rootWk = sqrt(Wkpos);
	double zMin, zMax;
	double c1, c2, c3;
	double epsrelDa0 = 0.0625 * epsrel;
	double epsabsDa0 = Da0 * epsrelDa0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Da0) || isnan(epsrel)) {
		return NAN;
	}
	else if (Da0 < uni->Da0Max) {
		if (branchNum == 0) {
			zMin = -1.0;
			zMax = uni->z_DaMax;
		}
		else {
			double z = uni->zLastTurn;
			double dummy;
			double Da0LastTurn = (DcT0Flat(uni, z, true, 1e-6, epsrel, 128, &dummy) \
				/ (1.0 + z));

			if (Da0 < Da0LastTurn) {
				return NAN;
			}
			else if (Da0 == Da0LastTurn) {
				return uni->zLastTurn;
			}

			zMin = uni->z_DaMax;
			zMax = uni->zLastTurn;
		}
	}
	else if (Da0 == uni->Da0Max) {
		return uni->z_DaMax;
	}
	else {
		fprintf(stderr, "Da0Inv_Cosmic: Da0 outside the invertable region.\n");
		return NAN;
	}

	if (branchNum == 0) {
		zCur = 0.5 * (zMin + zMax);
	}
	else {
		zCur = 2.0 * zMin;
	}

	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Because Da often has a maximum, we proceed using Halley's iteration instead of Newtons.
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;

		ainv = 1.0 + zCur;
		a = 1.0 / ainv;
		Dc0 = Dc0_Cosmic(uni, zCur, true, epsabsDa0, epsrelDa0, 128, &dummy);
		Einv_ = Einv(uni, ainv);
		sinhReduced = sinh(Dc0 * rootWk)/rootWk;
		coshReduced = cosh(Dc0 * rootWk);
		tanhReduced = tanh(Dc0 * rootWk)/rootWk;

		f =  FMA_m(sinhReduced, a, -Da0);
		fp = FMA_m(-sinhReduced, a, coshReduced * Einv_) * a;
		foverfp = f / fp;

		fppoverfp = FMA_m(FMA_m(c3, ainv, c2), ainv, c1) / fp;
		fppoverfp = FMA_m(fppoverfp, Einv_ * Einv_ * Einv_,
			FMA_m(Einv_ * Einv_, Wkpos, 2.0 * a) * tanhReduced);

		delta = - foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= zMin){
			if (zCur < zMin) {
				zCur = fabs(zCur - zMin) + zMin;
			}
			else {
				zCur = zMin + MIN(0.5, 0.125 * (zMax - zMin));
			}
			delta = zCur - zLast;
		}

		else if (zCur >= zMax) {
			if (zCur > zMax) {
				zCur = zMax - fabs(zMax - zCur);
			}
			else {
				zCur = zMax - MIN(0.5, 0.125 * (zMax - zMin));
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) <= epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double DaInv_Cosmic(UniverseLCDM *uni, double Da, unsigned int branchNum,
	double epsrel, unsigned int maxiter) {
	return Da0Inv_Cosmic(uni, Da / uni->DH, branchNum, epsrel, maxiter);
}

double Dl0Inv_Cosmic(UniverseLCDM *uni, double Dl0, unsigned int branchNum,
	double epsrel, unsigned int maxiter) {
	return (*(uni->F_Dl0Inv))(uni, Dl0, branchNum, epsrel, maxiter);
}

double Dl0InvFlat(struct UniverseLCDM *uni, double Dl0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, branchsign, f, ainv, Einv_, dummy;
	double fp, fppoverfp, foverfp, delta;
	double Dc0;
	double zMin, zMax;
	double c1, c2, c3;
	double epsrelDl0 = 0.0625 * epsrel;
	double epsabsDl0 = Dl0 * epsrelDl0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Dl0) || isnan(epsrel)) {
		return NAN;
	}
	else if (isinf(Dl0)) {
		return INFINITY;
	}
	else if ((!isnan(uni->Dc0LastTurn) && \
		Dl0 < uni->Dc0LastTurn * (1.0 + uni->zLastTurn)) || \
		isnan(uni->Dc0LastTurn)) {
		//do nothing
	}
	else {
		fprintf(stderr, "Dl0InvFlat: Unhandled case.\n");
		return NAN;
	}

	//Approximation for Dc0 = z * (1+z) / ( 1 + OmegaM * z )
	zCur = 0.5 * (1.0 - Dc0 * uni->OmegaMatter);
	zCur = sqrt(FMA_m(zCur, zCur, Dc0)) - zCur;
	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Iterate using modified Halley's method
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;

		ainv = 1.0 + zCur;
		Dc0 = Dc0_Cosmic(uni, zCur, true, epsabsDl0, epsrelDl0, 128, &dummy);
		Einv_ = Einv(uni, ainv);

		f = FMA_m(Dc0, ainv, -Dl0);
		fp = FMA_m(Einv_, ainv, Dc0);
		foverfp = f / fp;

		fppoverfp = FMA_m(FMA_m(FMA_m(c3, ainv, c2), ainv, c1),
			ainv * ainv * Einv_ * Einv_, 2.0) * Einv_;
		fppoverfp /= fp;

		delta = -foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= -1.0) {
			if (zCur < -1.0) {
				zCur = fabs(zCur + 1.0) - 1.0;
			}
			else {
				zCur = -0.5;
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) < epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double Dl0InvNegK(struct UniverseLCDM *uni, double Dl0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, branchsign, f, ainv, Einv_, dummy;
	double fp, fppoverfp, foverfp, delta;
	double tanReduced, sinReduced, cosReduced, Dc0;
	double Wkpos = -uni->OmegaCurvature;
	double rootWk = sqrt(Wkpos);
	double zMin, zMax;
	double c1, c2, c3;
	double epsrelDl0 = 0.0625 * epsrel;
	double epsabsDl0 = Dl0 * epsrelDl0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Dl0) || isnan(epsrel)) {
		return NAN;
	}
	else if (isinf(Dl0)) {
		return INFINITY;
	}
	else if ((!isnan(uni->Dc0LastTurn) && \
		Dl0 < uni->Dc0LastTurn * (1.0 + uni->zLastTurn)) || \
		isnan(uni->Dc0LastTurn)) {
		//do nothing
	}
	else {
		fprintf(stderr, "Dl0InvNegK: Unhandled case.\n");
		return NAN;
	}

	//Approximation for Dc0 = z * (1+z) / ( 1 + OmegaM * z )
	zCur = 0.5 * (1.0 - Dc0 * uni->OmegaMatter);
	zCur = sqrt(FMA_m(zCur, zCur, Dc0)) - zCur;
	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Iterate using modified Halley's method
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;

		ainv = 1.0 + zCur;
		Dc0 = Dc0_Cosmic(uni, zCur, true, epsabsDl0, epsrelDl0, 128, &dummy);
		Einv_ = Einv(uni, ainv);
		sinReduced = sin(Dc0 * rootWk)/rootWk;
		cosReduced = cos(Dc0 * rootWk);
		tanReduced = tan(Dc0 * rootWk)/rootWk;

		f = FMA_m(sinReduced, ainv, -Dl0);
		fp = FMA_m(Einv_ * cosReduced, ainv, sinReduced);
		foverfp = f / fp;

		fppoverfp = (FMA_m(FMA_m(c3, ainv, c2), ainv, c1) * \
			ainv * Einv_ * Einv_);
		fppoverfp /= FMA_m(Einv_, ainv, tanReduced);
		fppoverfp += FMA_m(tanReduced, ainv * Einv_ / -Wkpos, 2.0);
		fppoverfp *= Einv_;

		delta = -foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= -1.0) {
			if (zCur < -1.0) {
				zCur = fabs(zCur + 1.0) - 1.0;
			}
			else {
				zCur = -0.5;
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) < epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double Dl0InvPosK(struct UniverseLCDM *uni, double Dl0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, branchsign, f, ainv, Einv_, dummy;
	double fp, fppoverfp, foverfp, delta;
	double tanhReduced, sinhReduced, coshReduced, Dc0;
	double Wkpos = uni->OmegaCurvature;
	double rootWk = sqrt(Wkpos);
	double zMin, zMax;
	double c1, c2, c3;
	double epsrelDl0 = 0.0625 * epsrel;
	double epsabsDl0 = Dl0 * epsrelDl0;
	unsigned int iter;
	bool converged = false;

	if (isnan(Dl0) || isnan(epsrel)) {
		return NAN;
	}
	else if (isinf(Dl0)) {
		return INFINITY;
	}
	else if ((!isnan(uni->Dc0LastTurn) && \
		Dl0 < uni->Dc0LastTurn * (1.0 + uni->zLastTurn)) || \
		isnan(uni->Dc0LastTurn)) {
		//do nothing
	}
	else {
		fprintf(stderr, "Dl0InvPosK: Unhandled case.\n");
		return NAN;
	}

	//Approximation for Dc0 = z * (1+z) / ( 1 + OmegaM * z )
	zCur = 0.5 * (1.0 - Dc0 * uni->OmegaMatter);
	zCur = sqrt(FMA_m(zCur, zCur, Dc0)) - zCur;
	c3 = -2.0 * uni->OmegaRelativistic;
	c2 = -1.5 * uni->OmegaMatter;
	c1 = -uni->OmegaCurvature;

	//Iterate using modified Halley's method
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;

		ainv = 1.0 + zCur;
		Dc0 = Dc0_Cosmic(uni, zCur, true, epsabsDl0, epsrelDl0, 128, &dummy);
		Einv_ = Einv(uni, ainv);
		sinhReduced = sinh(Dc0 * rootWk)/rootWk;
		coshReduced = cosh(Dc0 * rootWk);
		tanhReduced = tanh(Dc0 * rootWk)/rootWk;

		f = FMA_m(sinhReduced, ainv, -Dl0);
		fp = FMA_m(Einv_ * coshReduced, ainv, sinhReduced);
		foverfp = f / fp;

		fppoverfp = (FMA_m(FMA_m(c3, ainv, c2), ainv, c1) * \
			ainv * Einv_ * Einv_);
		fppoverfp /= FMA_m(Einv_, ainv, tanhReduced);
		fppoverfp += FMA_m(tanhReduced, ainv * Einv_ / Wkpos, 2.0);
		fppoverfp *= Einv_;

		delta = -foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-4);
		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= -1.0) {
			if (zCur < -1.0) {
				zCur = fabs(zCur + 1.0) - 1.0;
			}
			else {
				zCur = -0.5;
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) < epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double DlInv_Cosmic(UniverseLCDM *uni, double Dl, unsigned int branchNum,
	double epsrel, unsigned int maxiter) {
	return Dl0Inv_Cosmic(uni, Dl / uni->DH, branchNum, epsrel, maxiter);
}

double Vc0Inv_Cosmic(UniverseLCDM *uni, double Vc0, unsigned int branchNum,
	double epsrel, unsigned int maxiter) {

	return (*(uni->F_Vc0Inv))(uni, Vc0, branchNum, epsrel, maxiter);
}

double VcInv_Cosmic(UniverseLCDM *uni, double Vc,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double Vc0 = Vc / uni->VH;

	return Vc0Inv_Cosmic(uni, Vc0, branchNum, epsrel, maxiter);
}

static double Vc0InvFlat(struct UniverseLCDM *uni, double Vc0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {

	return Dc0Inv_Cosmic(uni, pow(Vc0 * 3.0 / (4.0 * Pi), 1.0/3.0),
		branchNum, epsrel, maxiter);
}

static double Vc0InvNegK(struct UniverseLCDM *uni, double Vc0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double OmegaKpos = -uni->OmegaCurvature;
	double nrm = sqrt(OmegaKpos);
	double x, f, fp, delta;
	unsigned int iter;
	bool converged = false;
	double xCur, xLast, y, sx;
	double y0, dy;
	double epsrelx = epsrel * 0.25;

	y = -OmegaKpos * nrm * Vc0 / Pi; //This is sin(x) - x
	y0 = -2.0 * Pi * round(y / (2*Pi));
	dy = y - y0;

	xCur = -y0 - copysign(pow(6.0 * fabs(dy), 1.0/3.0), dy) ;

	//Iterate using Newton's method
	for (iter = 0; iter < maxiter; ++iter) {
		xLast = xCur;
		sx = sin(xCur);
		f = sx - xCur - y;
		fp = sin(0.5*xCur);
		fp *= -2.0 * fp;

		delta = -f/fp;

		xCur += delta;

		if (fabs(delta) <= epsrelx * fabs(xCur)) {
			converged = true;
			break;
		}
	}

	if (converged) {
		return Dc0Inv_Cosmic(uni, xCur * 0.5 * nrm / OmegaKpos, branchNum, epsrel, maxiter);
	}
	return NAN;
}

static double Vc0InvPosK(struct UniverseLCDM *uni, double Vc0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double OmegaKpos = uni->OmegaCurvature;
	double nrm = sqrt(OmegaKpos);
	double x, f, fp, delta;
	unsigned int iter;
	bool converged = false;
	double xCur, xLast, y, ypos, sx;
	double epsrelx = epsrel * 0.25;

	y = OmegaKpos * nrm * Vc0 / Pi; //This is sinh(x) - x
	ypos = fabs(y);

	if (ypos < 4.5) {
		xCur = copysign(pow(6.0 * ypos, 1.0/3.0), y);
	} else {
		xCur = copysign(y, log(0.5*ypos));
	}

	//Iterate using Newton's method
	for (iter = 0; iter < maxiter; ++iter) {
		xLast = xCur;
		sx = sinh(xCur);
		f = sx - xCur - y;
		fp = sinh(0.5*xCur);
		fp *= 2.0 * fp;

		delta = -f/fp;
		xCur += delta;

		if (fabs(delta) <= epsrelx * fabs(xCur)) {
			converged = true;
			break;
		}
	}

	if (converged) {
		return Dc0Inv_Cosmic(uni, xCur * 0.5 * nrm / OmegaKpos, branchNum, epsrel, maxiter);
	}
	return NAN;
}

double tL0Inv_Cosmic(UniverseLCDM *uni, double tL0,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double zCur, zLast, f, fp, dummy;
	double ainv, foverfp, fppoverfp, delta, c4, c3, c2, c0;
	double epsreltL0 = 0.0625 * epsrel;
	double epsabstL0 = tL0 * epsreltL0;
	unsigned int iter;
	bool converged = false;

	if (isnan(tL0) || isnan(epsrel)) {
		return NAN;
	}
	else if (tL0 < uni->age0) {
		//do nothing
	}
	else if (!isnan(uni->tL0LastTurn) && tL0 == uni->tL0LastTurn) {
		return uni->zLastTurn;
	}
	else if (tL0 == uni->age0) {
		return INFINITY;
	}
	else if (isnan(tL0) || isnan(epsrel)) {
		return NAN;
	}
	else {
		fprintf(stderr, "tL0Inv_Cosmic: tL0 outside the invertable region.\n");
		return NAN;
	}

	c4 = -3.0 * uni->OmegaRelativistic;
	c3 = -2.5 * uni->OmegaMatter;
	c2 = -2.0 * uni->OmegaCurvature;
	c0 = -uni->OmegaLambda;

	//Approximation for tL= tAge * ( 1 - 1/pow(1+z, 1.5) )
	zCur = pow(1.0 - tL0 / uni->age0, -2.0/3.0);

	//Iterate using Halley's method
	for (iter = 0; iter < maxiter; ++iter) {
		zLast = zCur;
		f = (tL0_Cosmic(uni, zCur, true, epsabstL0, epsreltL0, MaxGSLlevel, &dummy)
			- tL0);
		ainv = 1.0 + zCur;
		fp = Einv(uni, ainv) / ainv;

		fppoverfp = FMA_m(FMA_m(FMA_m(c4, ainv, c3), ainv, c2), ainv * ainv, c0) * ainv;
		fppoverfp *= fp * fp;
		foverfp = f / fp;
		delta = -foverfp / MAX(FMA_m(-0.5 * fppoverfp, foverfp, 1.0), 0x1.0p-1);

		zCur += delta;

		//Prevent the current estimate from becoming invalid
		if (zCur <= -1.0){
			if (zCur < -1.0) {
				zCur = fabs(zCur + 1.0) - 1.0;
			}
			else {
				zCur = -0.5;
			}
			delta = zCur - zLast;
		}

		if (fabs(delta) <= epsrel * fabs(zCur)) {
			converged = true;
			break;
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}

double tLInv_Cosmic(UniverseLCDM *uni, double tL,
	unsigned int branchNum, double epsrel, unsigned int maxiter) {
	double tL0 = tL / uni->tH;

	return tL0Inv_Cosmic(uni, tL0, branchNum, epsrel, maxiter);
}

//Function to find the redshift at which Da has a maximum
static double FindDa0Max(UniverseLCDM *uni, double epsrel, unsigned int maxiter) {
	double Einv_, ainv, a, zCur, zLast, dummy;
	double epsrelDa0 = 0.0625 * epsrel;
	unsigned int iter;
	bool converged = false;

	if (uni->flat || fabs(uni->OmegaCurvature) <= NumericallyFlat) {
		//Iterate using Newton's method
		double f, fp, delta, Dc0, dummy;
		double c3, c2, c1;

		zCur = 0.0;
		f = (DcT0Flat(uni, 1.0, true, epsrelDa0, epsrel, 128, &dummy) \
				/ ainv);

		c3 = -2.0 * uni->OmegaRelativistic;
		c2 = -1.5 * uni->OmegaMatter;
		c1 = -uni->OmegaCurvature;
		for (iter = 0; iter < maxiter; ++iter) {
			zLast = zCur;

			ainv = 1.0 + zCur;
			a = 1.0 / ainv;
			Dc0 = DcT0Flat(uni, zCur, \
				true, f * epsrelDa0, epsrel, 128, &dummy);
			Einv_ = Einv(uni, ainv);
			f = FMA_m(-Dc0, a, Einv_);
			fp = FMA_m(FMA_m(c3, ainv, c2), ainv, c1);
			fp = FMA_m(fp, ainv * Einv_ * Einv_ * Einv_, \
				FMA_m(Dc0, a, -Einv_) * a * 2.0);

			delta = -f / fp;
			zCur += delta;

			//Prevent the current estimate from becoming invalid
			if (zCur <= -1.0){
				if (zCur < -1.0) {
					zCur = fabs(zCur + 1.0) - 1.0;
				}
				else {
					zCur = -0.5;
				}
			}

			if (fabs(delta) <= epsrel * fabs(zCur) || isinf(zCur)) {
				converged = true;
				break;
			}
		}
	}

	else if (uni->OmegaCurvature < 0.0) {
		//Iterate using Newton's method
		double f, fp, delta, Dc0, dummy;
		double c3, c2, c1;
		double tanReduced, Einv_, ainv, a;
		double Wkpos = -uni->OmegaCurvature;
		double rootWk = sqrt(Wkpos);

		zCur = 0.0;
		Dc0 = (Dc0_Cosmic(uni, 1.0, true, epsrelDa0, epsrel, 128, &dummy) \
				/ ainv);

		c3 = -2.0 * uni->OmegaRelativistic;
		c2 = -1.5 * uni->OmegaMatter;
		c1 = -uni->OmegaCurvature;
		for (iter = 0; iter < maxiter; ++iter) {
			zLast = zCur;

			ainv = 1.0 + zCur;
			a = 1.0 / ainv;
			Dc0 = Dc0_Cosmic(uni, zCur, true,
				epsrelDa0 * Dc0, epsrelDa0, 128, &dummy);
			tanReduced = tan(Dc0 * rootWk) / rootWk;
			Einv_ = Einv(uni, ainv);

			f = FMA_m(-tanReduced, a, Einv_);
			fp = FMA_m(2.0*tanReduced, a, -2.0 * Einv_);
			fp = FMA_m(fp, a, Einv_*Einv_ * \
				FMA_m(FMA_m(FMA_m(c3, ainv, c2), ainv, c1), \
					ainv * Einv_, Wkpos * tanReduced));

			delta = -f / fp;
			zCur += delta;

			//Prevent the current estimate from becoming invalid
			if (zCur <= -1.0){
				if (zCur < -1.0) {
					zCur = fabs(zCur + 1.0) - 1.0;
				}
				else {
					zCur = -0.5;
				}
			}

			if (fabs(delta) <= epsrel * fabs(zCur) || isinf(zCur)) {
				converged = true;
				break;
			}
		}
	}

	else { //uni->OmegaCurvature > 0.0
		//Iterate using Newton's method
		double f, fp, delta, Dc0, dummy;
		double c3, c2, c1;
		double tanhReduced, Einv_, ainv, a;
		double Wkpos = uni->OmegaCurvature;
		double rootWk = sqrt(Wkpos);

		zCur = 0.0;
		Dc0 = (Dc0_Cosmic(uni, 1.0, true, epsrelDa0, epsrel, 128, &dummy) \
				/ ainv);

		c3 = -2.0 * uni->OmegaRelativistic;
		c2 = -1.5 * uni->OmegaMatter;
		c1 = -uni->OmegaCurvature;
		for (iter = 0; iter < maxiter; ++iter) {
			zLast = zCur;

			ainv = 1.0 + zCur;
			a = 1.0 / ainv;
			Dc0 = Dc0_Cosmic(uni, zCur, true,
				epsrelDa0 * Dc0, epsrelDa0, 128, &dummy);
			tanhReduced = tanh(Dc0 * rootWk) / rootWk;
			Einv_ = Einv(uni, ainv);

			f = FMA_m(-tanhReduced, a, Einv_);
			fp = FMA_m(2.0*tanhReduced, a, -2.0 * Einv_);
			fp = FMA_m(fp, a, Einv_*Einv_ * \
				FMA_m(FMA_m(FMA_m(c3, ainv, c2), ainv, c1), \
					ainv * Einv_, Wkpos * tanhReduced));

			delta = -f / fp;
			zCur += delta;

			//Prevent the current estimate from becoming invalid
			if (zCur <= -1.0){
				if (zCur < -1.0) {
					zCur = fabs(zCur + 1.0) - 1.0;
				}
				else {
					zCur = -0.5;
				}
			}

			if (fabs(delta) <= epsrel * fabs(zCur) || isinf(zCur)) {
				converged = true;
				break;
			}
		}
	}

	if (!converged) {
		zCur = NAN;
	}

	return zCur;
}


text:
/* filter/gaussian.c
 *
 * Gaussian smoothing filters
 *
 * Copyright (C) 2018 Patrick Alken
 *
 * 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 3 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.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_poly.h>

/* maximum derivative order allowed for Gaussian filter */
#define GSL_FILTER_GAUSSIAN_MAX_ORDER     10

typedef double gaussian_type_t;
typedef double ringbuf_type_t;
#include "ringbuf.c"

typedef struct
{
  size_t n;        /* window size */
  double * window; /* linear array with current window */
  ringbuf * rbuf;  /* ring buffer storing current window */
} gaussian_state_t;

static size_t gaussian_size(const size_t n);
static int gaussian_init(const size_t n, void * vstate);
static int gaussian_insert(const gaussian_type_t x, void * vstate);
static int gaussian_delete(void * vstate);
static int gaussian_get(void * params, gaussian_type_t * result, const void * vstate);

static const gsl_movstat_accum gaussian_accum_type;

/*
gsl_filter_gaussian_alloc()
  Allocate a workspace for Gaussian filtering.

Inputs: K - number of samples in window; if even, it is rounded up to
            the next odd, to have a symmetric window

Return: pointer to workspace
*/

gsl_filter_gaussian_workspace *
gsl_filter_gaussian_alloc(const size_t K)
{
  const size_t H = K / 2;
  gsl_filter_gaussian_workspace *w;
  size_t state_size;

  w = calloc(1, sizeof(gsl_filter_gaussian_workspace));
  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->K = 2 * H + 1;

  w->kernel = malloc(w->K * sizeof(double));
  if (w->kernel == 0)
    {
      gsl_filter_gaussian_free(w);
      GSL_ERROR_NULL ("failed to allocate space for kernel", GSL_ENOMEM);
      return NULL;
    }

  state_size = gaussian_size(w->K);

  w->movstat_workspace_p = gsl_movstat_alloc_with_size(state_size, H, H);
  if (!w->movstat_workspace_p)
    {
      gsl_filter_gaussian_free(w);
      GSL_ERROR_NULL ("failed to allocate space for movstat workspace", GSL_ENOMEM);
    }

  return w;
}

void
gsl_filter_gaussian_free(gsl_filter_gaussian_workspace * w)
{
  if (w->kernel)
    free(w->kernel);

  if (w->movstat_workspace_p)
    gsl_movstat_free(w->movstat_workspace_p);

  free(w);
}

/*
gsl_filter_gaussian()
  Apply a Gaussian filter to an input vector:

G_{sigma}(x) = exp [ -x^2 / (2 sigma^2) ]

Inputs: alpha - number of standard deviations to include in Gaussian kernel
        order - derivative order of Gaussian
        x     - input vector, size n
        y     - (output) filtered vector, size n
        w     - workspace

Notes:
1) If alpha = 3, then the Gaussian kernel will be a Gaussian of +/- 3 standard deviations
*/

int
gsl_filter_gaussian(const gsl_filter_end_t endtype, const double alpha, const size_t order, const gsl_vector * x,
                    gsl_vector * y, gsl_filter_gaussian_workspace * w)
{
  if (x->size != y->size)
    {
      GSL_ERROR("input and output vectors must have same length", GSL_EBADLEN);
    }
  else if (alpha <= 0.0)
    {
      GSL_ERROR("alpha must be positive", GSL_EDOM);
    }
  else
    {
      int status;
      gsl_vector_view kernel = gsl_vector_view_array(w->kernel, w->K);

      /* construct Gaussian kernel of length K */
      gsl_filter_gaussian_kernel(alpha, order, 1, &kernel.vector);

      status = gsl_movstat_apply_accum(endtype, x, &gaussian_accum_type, (void *) w->kernel, y,
                                       NULL, w->movstat_workspace_p);

      return status;
    }
}

/*
gsl_filter_gaussian_kernel()
  Construct Gaussian kernel with given sigma and order

Inputs: alpha     - number of standard deviations to include in window
        order     - kernel order (0 = gaussian, 1 = first derivative, ...)
        normalize - normalize so sum(G) = 1
        kernel    - (output) Gaussian kernel

Return: success/error

Notes:
1) If alpha = 3, then the output kernel will contain a Gaussian with +/- 3 standard deviations
*/

int
gsl_filter_gaussian_kernel(const double alpha, const size_t order, const int normalize, gsl_vector * kernel)
{
  const size_t N = kernel->size;

  if (alpha <= 0.0)
    {
      GSL_ERROR("alpha must be positive", GSL_EDOM);
    }
  else if (order > GSL_FILTER_GAUSSIAN_MAX_ORDER)
    {
      GSL_ERROR("derivative order is too large", GSL_EDOM);
    }
  else
    {
      const double half = 0.5 * (N - 1.0); /* (N - 1) / 2 */
      double sum = 0.0;
      size_t i;

      /* check for quick return */
      if (N == 1)
        {
          if (order == 0)
            gsl_vector_set(kernel, 0, 1.0);
          else
            gsl_vector_set(kernel, 0, 0.0);

          return GSL_SUCCESS;
        }

      for (i = 0; i < N; ++i)
        {
          double xi = ((double)i - half) / half;
          double yi = alpha * xi;
          double gi = exp(-0.5 * yi * yi);

          gsl_vector_set(kernel, i, gi);
          sum += gi;
        }

      /* normalize so sum(kernel) = 1 */
      if (normalize)
        gsl_vector_scale(kernel, 1.0 / sum);

      if (order > 0)
        {
          const double beta = -0.5 * alpha * alpha;
          double q[GSL_FILTER_GAUSSIAN_MAX_ORDER + 1];
          size_t k;

          /*
           * Need to calculate derivatives of the Gaussian window; define
           *
           * w(n) = C * exp [ p(n) ]
           *
           * p(n) = beta * n^2
           * beta = -1/2 * ( alpha / ((N-1)/2) )^2
           *
           * Then:
           *
           * d^k/dn^k w(n) = q_k(n) * w(n)
           *
           * where q_k(n) is a degree-k polynomial in n, which satisfies:
           *
           * q_k(n) = d/dn q_{k-1}(n) + q_{k-1}(n) * dp(n)/dn
           * q_0(n) = 1 / half^{order}
           */

          /* initialize q_0(n) = 1 / half^{order} */
          q[0] = 1.0 / gsl_pow_uint(half, order);
          for (i = 1; i <= GSL_FILTER_GAUSSIAN_MAX_ORDER; ++i)
            q[i] = 0.0;

          /* loop through derivative orders and calculate q_k(n) for k = 1,...,order */
          for (k = 1; k <= order; ++k)
            {
              double qm1 = q[0];

              q[0] = q[1];
              for (i = 1; i <= k; ++i)
                {
                  double tmp = q[i];
                  q[i] = (i + 1.0) * q[i + 1] + /* d/dn q_{k-1} */
                         2.0 * beta * qm1;      /* q_{k-1}(n) p'(n) */
                  qm1 = tmp;
                }
            }

          /* now set w(n) := q(n) * w(n) */
          for (i = 0; i < N; ++i)
            {
              double xi = ((double)i - half) / half;
              double qn = gsl_poly_eval(q, order + 1, xi);
              double *wn = gsl_vector_ptr(kernel, i);

              *wn *= qn;
            }
        }

      return GSL_SUCCESS;
    }
}

static size_t
gaussian_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(gaussian_state_t);
  size += n * sizeof(gaussian_type_t);
  size += ringbuf_size(n);

  return size;
}

static int
gaussian_init(const size_t n, void * vstate)
{
  gaussian_state_t * state = (gaussian_state_t *) vstate;

  state->n = n;

  state->window = (gaussian_type_t *) ((unsigned char *) vstate + sizeof(gaussian_state_t));
  state->rbuf = (ringbuf *) ((unsigned char *) state->window + n * sizeof(gaussian_type_t));

  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
gaussian_insert(const gaussian_type_t x, void * vstate)
{
  gaussian_state_t * state = (gaussian_state_t *) vstate;

  /* add new element to ring buffer */
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
gaussian_delete(void * vstate)
{
  gaussian_state_t * state = (gaussian_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    ringbuf_pop_back(state->rbuf);

  return GSL_SUCCESS;
}

static int
gaussian_get(void * params, gaussian_type_t * result, const void * vstate)
{
  const gaussian_state_t * state = (const gaussian_state_t *) vstate;
  const double * kernel = (const double *) params;
  size_t n = ringbuf_copy(state->window, state->rbuf);
  double sum = 0.0;
  size_t i;

  for (i = 0; i < n; ++i)
    sum += state->window[i] * kernel[n - i - 1];

  *result = sum;

  return GSL_SUCCESS;
}

static const gsl_movstat_accum gaussian_accum_type =
{
  gaussian_size,
  gaussian_init,
  gaussian_insert,
  gaussian_delete,
  gaussian_get
};


text:
#include <stdlib.h>
#include <limits.h>
#include <stdio.h>
#include <errno.h>
#include <time.h>
#include <sys/resource.h>
#include <string.h>

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

double drand() {
    return (double) rand() / RAND_MAX;
}

void ptime(struct rusage *ru0, struct rusage *ru1) {

    double utime = 0, stime = 0, ttime = 0;

    utime = (double) ru1->ru_utime.tv_sec
            + 1.e-6 * (double) ru1->ru_utime.tv_usec
            - ru0->ru_utime.tv_sec
            - 1.e-6 * (double) ru0->ru_utime.tv_usec;
    stime = (double) ru1->ru_stime.tv_sec
            + 1.e-6 * (double) ru1->ru_stime.tv_usec
            - ru0->ru_stime.tv_sec
            - 1.e-6 * (double) ru0->ru_stime.tv_usec;
    ttime = stime + utime;

    printf("%3f\n", ttime);
}

int main(int argc, char *argv[]) {

    if (argc < 2) {
        errno = EINVAL;
        perror("");
        return errno;
    }
    int n = strtol(argv[1], NULL, 10);
    if (0 != errno) {
        perror("");
        return errno;
    }

    srand(time(NULL));

    int els = n * n;
    double *A_data = (double *) malloc(els * sizeof(double));
    for (int i = 0; i < els; i++) {
        A_data[i] = drand();
    }
    double *b_data = (double *) malloc(n * sizeof(double));
    for (int i = 0; i < n; i++) {
        b_data[i] = drand();
    }

    gsl_matrix_view A
            = gsl_matrix_view_array(A_data, n, n);
    gsl_matrix *A_copy = gsl_matrix_alloc(n, n);
    gsl_matrix_memcpy(A_copy, &A.matrix);
    fprintf(stderr, "A = [\n");
    gsl_matrix_fprintf(stderr, &A.matrix, "%g");
    fprintf(stderr, "]\n");

    gsl_vector_view b
            = gsl_vector_view_array(b_data, n);
    fprintf(stderr, "b = [\n");
    gsl_vector_fprintf(stderr, &b.vector, "%g");
    fprintf(stderr, "]\n");

    gsl_vector *x = gsl_vector_alloc(n);

    struct rusage t0, t1, t2;
    int s;

    gsl_permutation *p = gsl_permutation_alloc(n);

    getrusage(RUSAGE_SELF, &t0);
    gsl_linalg_LU_decomp(&A.matrix, p, &s);
    getrusage(RUSAGE_SELF, &t1);
    gsl_linalg_LU_solve(&A.matrix, p, &b.vector, x);
    getrusage(RUSAGE_SELF, &t2);

    fprintf(stderr, "x = [\n");
    gsl_vector_fprintf(stderr, x, "%g");
    fprintf(stderr, "]\n");

    fprintf(stderr, "%s \n", "decomposition time:");
    ptime(&t0, &t1);
    fprintf(stderr, "%s \n", "solve time:");
    ptime(&t1, &t2);
    fprintf(stderr, "%s \n", "together time:");
    ptime(&t0, &t2);

    gsl_vector *y = gsl_vector_alloc(n);
    gsl_blas_dgemv(CblasNoTrans, 1.0, A_copy, x, 0, y);
    fprintf(stderr, "y = [\n");
    gsl_vector_fprintf(stderr, y, "%g");
    fprintf(stderr, "]\n");

    gsl_permutation_free(p);
    gsl_vector_free(x);
    return EXIT_SUCCESS;
}

max_stars_repo_name: mistyfiky/agh-mownit

max_issues_repo_name: mistyfiky/agh-mownit

max_forks_repo_name: mistyfiky/agh-mownit