text:
/* trav.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_bst.h>
#include <gsl/gsl_errno.h>

int
gsl_bst_trav_init(gsl_bst_trav * trav, const gsl_bst_workspace * w)
{
  int status = (w->type->trav_init)((void *) &trav->trav_data, (const void *) &w->table);
  trav->type = w->type;
  return status;
}

/*
gsl_bst_trav_first()
  Initialize traverser to least-valued item in tree and return
a pointer to it. Return NULL if tree has no nodes.
*/

void *
gsl_bst_trav_first(gsl_bst_trav * trav, const gsl_bst_workspace * w)
{
  trav->type = w->type;
  return (w->type->trav_first)((void *) &trav->trav_data, (const void *) &w->table);
}

/*
gsl_bst_trav_last()
  Initializes |trav| for |tree| and selects and returns a pointer
to its greatest-valued item. Returns NULL if |tree| contains no nodes.
*/

void *
gsl_bst_trav_last (gsl_bst_trav * trav, const gsl_bst_workspace * w)
{
  trav->type = w->type;
  return (w->type->trav_last)((void * ) &trav->trav_data, (const void *) &w->table);
}

/*
gsl_bst_trav_find()
  Searches for |item| in tree. If found, initializes |trav| to
the item found and returns the item as well.  If there is no matching
item, initializes |trav| to the null item and returns |NULL|.
*/

void *
gsl_bst_trav_find (const void * item, gsl_bst_trav * trav, const gsl_bst_workspace * w)
{
  trav->type = w->type;
  return (w->type->trav_find)(item, (void * ) &trav->trav_data, (const void *) &w->table);
}

/*
gsl_bst_trav_insert()
  Attempts to insert |item| into tree.  If |item| is inserted
successfully, it is returned and |trav| is initialized to its location.
If a duplicate is found, it is returned and |trav| is initialized to
its location.  No replacement of the item occurs.
If a memory allocation failure occurs, |NULL| is returned and |trav|
is initialized to the null item.
*/

void *
gsl_bst_trav_insert (void * item, gsl_bst_trav * trav, gsl_bst_workspace * w)
{
  trav->type = w->type;
  return (w->type->trav_insert)(item, (void * ) &trav->trav_data, (void *) &w->table);
}

/*
gsl_bst_trav_copy()
  Copy traverser 'src' into 'dest'
*/

void *
gsl_bst_trav_copy(gsl_bst_trav * dest, const gsl_bst_trav * src)
{
  dest->type = src->type;
  return (src->type->trav_copy)((void * ) &dest->trav_data, (const void *) &src->trav_data);
}

/*
gsl_bst_trav_next()
  Update traverser to point to next sequential node, and
return a pointer to its data
*/

void *
gsl_bst_trav_next(gsl_bst_trav * trav)
{
  return (trav->type->trav_next)((void *) &trav->trav_data);
}

/*
gsl_bst_trav_prev()
  Update traverser to point to previous sequential node, and
return a pointer to its data
*/

void *
gsl_bst_trav_prev(gsl_bst_trav * trav)
{
  return (trav->type->trav_prev)((void *) &trav->trav_data);
}

/*
gsl_bst_trav_cur()
  Return a pointer to data of current traverser node
*/

void *
gsl_bst_trav_cur(const gsl_bst_trav * trav)
{
  return (trav->type->trav_cur)((const void *) &trav->trav_data);
}

/*
gsl_bst_trav_replace()
  Replace current item in trav with new_item and returns the item
replaced. The new item must not change the ordering of the tree.
*/

void *
gsl_bst_trav_replace (gsl_bst_trav * trav, void * new_item)
{
  return (trav->type->trav_replace)((void * ) &trav->trav_data, new_item);
}


text:
/**************************************************************/
/*! \brief This file contains all .h files
 */
/**************************************************************/

#ifndef __SLIM_H__
#define __SLIM_H__

#define _GNU_SOURCE 1

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <time.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <bcls.h>
#include <GKlib.h>
#include <def.h>
#include <struct.h>
#include <proto.h>
#include <mpi.h>
#include <omp.h>

#endif

ext: h

max_stars_repo_name: echristakopoulou/glslim

max_issues_repo_name: echristakopoulou/glslim

max_forks_repo_name: echristakopoulou/glslim


text:
#ifndef __gslInterp_
#define __gslInterp__

#include "include/stdinc.h"
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_integration.h>
#include <vector>

class gslInterp {
protected:
  gsl_interp_accel *acc;
  gsl_interp       *interp;
  vector<double>   x, y;
public:

  gslInterp(vector<double> &x0, vector<double> &y0) {
    acc   = gsl_interp_accel_alloc();
    for (size_t i = 0; i < x0.size(); i++) {
      x.push_back(x0[i]);
      y.push_back(y0[i]);
    }
    interp = gsl_interp_alloc(gsl_interp_linear, x.size());
    gsl_interp_init(interp, &x[0], &y[0], x.size());
  }
  ~gslInterp() {
    gsl_interp_accel_free(acc);
    gsl_interp_free(interp);
  }
  void init(vector<double> &x0, vector<double> &y0) {
    gsl_interp_accel_free(acc);
    gsl_interp_free(interp);
    acc    = gsl_interp_accel_alloc();
    x.clear();
    y.clear();
    for (size_t i = 0; i < x0.size(); i++) {
      x.push_back(x0[i]);
      y.push_back(y0[i]);
    }
    interp = gsl_interp_alloc(gsl_interp_linear, x.size());
    gsl_interp_init(interp, &x[0], &y[0], x.size());
  }

  inline double eval(double xc) {
    return gsl_interp_eval(interp, &x[0], &y[0], xc, acc);
  }
  inline double deriv(double xc) {
    return gsl_interp_eval_deriv(interp, &x[0], &y[0], xc, acc);
  }
  inline double deriv2(double xc) {
    return gsl_interp_eval_deriv2(interp, &x[0], &y[0], xc, acc);
  }
  inline double integ(double x_lo, double x_up) {
    return gsl_interp_eval_integ(interp, &x[0], &y[0], x_lo, x_up, acc);
  }
};

#endif // __gslInterp__

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
/* specfunc/test_gamma.c
 *
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"

int test_gamma(void)
{
  gsl_sf_result r;
  gsl_sf_result r1, r2;
  double sgn;
  int s = 0;

  TEST_SF(s,  gsl_sf_lngamma_e, (-0.1, &r), 2.368961332728788655 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (-1.0/256.0, &r), 5.547444766967471595  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (1.0e-08, &r), 18.420680738180208905 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (0.1, &r), 2.252712651734205 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (1.0 + 1.0/256.0, &r), -0.0022422226599611501448 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (2.0 + 1.0/256.0, &r), 0.0016564177556961728692 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (100.0, &r), 359.1342053695753 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (-1.0-1.0/65536.0, &r), 11.090348438090047844 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (-1.0-1.0/268435456.0, &r), 19.408121054103474300 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (-100.5, &r), -364.9009683094273518 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lngamma_e, (-100-1.0/65536.0, &r), -352.6490910117097874 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (0.7, &r, &sgn), 0.26086724653166651439, TEST_TOL1, 1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (0.1, &r, &sgn), 2.2527126517342059599, TEST_TOL0, 1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (-0.1, &r, &sgn), 2.368961332728788655, TEST_TOL0, -1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (-1.0-1.0/65536.0, &r, &sgn), 11.090348438090047844, TEST_TOL0, 1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (-2.0-1.0/256.0, &r, &sgn), 4.848447725860607213, TEST_TOL0, -1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (-2.0-1.0/65536.0, &r, &sgn), 10.397193628164674967, TEST_TOL0, -1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (-3.0-1.0/8.0, &r, &sgn), 0.15431112768404182427, TEST_TOL2, 1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lngamma_sgn_e, (-100.5, &r, &sgn), -364.9009683094273518, TEST_TOL0, -1.0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gamma_e, (1.0 + 1.0/4096.0, &r), 0.9998591371459403421 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (1.0 + 1.0/32.0, &r), 0.9829010992836269148 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (2.0 + 1.0/256.0, &r), 1.0016577903733583299 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (9.0, &r), 40320.0                   , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (10.0, &r), 362880.0                  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (100.0, &r), 9.332621544394415268e+155 , TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (170.0, &r), 4.269068009004705275e+304 , TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (171.0, &r), 7.257415615307998967e+306 , TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (-10.5, &r), -2.640121820547716316e-07  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gamma_e, (-11.25, &r), 6.027393816261931672e-08  , TEST_TOL0, GSL_SUCCESS); /* exp()... not my fault */
  TEST_SF(s,  gsl_sf_gamma_e, (-1.0+1.0/65536.0, &r), -65536.42280587818970 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gammastar_e, (1.0e-08, &r), 3989.423555759890865  , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (1.0e-05, &r), 126.17168469882690233 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (0.001, &r), 12.708492464364073506 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (1.5, &r), 1.0563442442685598666 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (3.0, &r), 1.0280645179187893045 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (9.0, &r), 1.0092984264218189715 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (11.0, &r), 1.0076024283104962850 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (100.0, &r), 1.0008336778720121418 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (1.0e+05, &r), 1.0000008333336805529 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammastar_e, (1.0e+20, &r), 1.0 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gammainv_e, (1.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (2.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (3.0, &r), 0.5, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (4.0, &r), 1.0/6.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gammainv_e, (10.0, &r), 1.0/362880.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (100.0, &r), 1.0715102881254669232e-156, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gammainv_e, (0.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (-1.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (-2.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (-3.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (-4.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_gammainv_e, (-10.5, &r), -1.0/2.640121820547716316e-07, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (-11.25, &r), 1.0/6.027393816261931672e-08, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_gammainv_e, (-1.0+1.0/65536.0, &r), -1.0/65536.42280587818970 , TEST_TOL1, GSL_SUCCESS);

  TEST_SF_2(s, gsl_sf_lngamma_complex_e, (5.0, 2.0, &r1, &r2),
            2.7487017561338026749, TEST_TOL0,
            3.0738434100497007915, TEST_TOL0,
            GSL_SUCCESS);

  TEST_SF_2(s, gsl_sf_lngamma_complex_e, (100.0, 100.0, &r1, &r2),
            315.07804459949331323, TEST_TOL1,
            2.0821801804113110099, TEST_TOL3,
            GSL_SUCCESS);

  TEST_SF_2(s, gsl_sf_lngamma_complex_e, (100.0, -1000.0, &r1, &r2),
            -882.3920483010362817000, TEST_TOL1,
            -2.1169293725678813270, TEST_TOL3,
            GSL_SUCCESS);

  TEST_SF_2(s, gsl_sf_lngamma_complex_e, (-100.0, -1.0, &r1, &r2),
            -365.0362469529239516000, TEST_TOL1,
            -3.0393820262864361140, TEST_TOL1,
            GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_taylorcoeff_e, (10,   1.0/1048576.0, &r), 1.7148961854776073928e-67  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (10,   1.0/1024.0, &r), 2.1738891788497900281e-37  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (10,   1.0, &r), 2.7557319223985890653e-07  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (10,   5.0, &r), 2.6911444554673721340      , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (10,   500.0, &r), 2.6911444554673721340e+20  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (100,  100.0, &r), 1.0715102881254669232e+42  , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (1000, 200.0, &r), 2.6628790558154746898e-267 , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_taylorcoeff_e, (1000, 500.0, &r), 2.3193170139740855074e+131 , TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s,   gsl_sf_fact_e, (0, &r), 1.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,   gsl_sf_fact_e, (1, &r), 1.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,   gsl_sf_fact_e, (7, &r), 5040.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_fact_e, (33, &r), 8.683317618811886496e+36 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,   gsl_sf_doublefact_e, (0, &r), 1.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,   gsl_sf_doublefact_e, (1, &r), 1.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,   gsl_sf_doublefact_e, (7, &r), 105.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_doublefact_e, (33, &r), 6.332659870762850625e+18 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_lnfact_e, (0, &r), 0.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnfact_e, (1, &r), 0.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnfact_e, (7, &r), 8.525161361065414300 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnfact_e, (33, &r), 85.05446701758151741 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_lndoublefact_e, (0, &r), 0.0  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lndoublefact_e, (7, &r), 4.653960350157523371  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lndoublefact_e, (33, &r), 43.292252022541719660 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lndoublefact_e, (34, &r), 45.288575519655959140 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lndoublefact_e, (1034, &r), 3075.6383796271197707 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lndoublefact_e, (1035, &r), 3078.8839081731809169 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_lnchoose_e, (7,3, &r), 3.555348061489413680 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnchoose_e, (5,2, &r), 2.302585092994045684 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_choose_e, (7,3, &r), 35.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_choose_e, (7,4, &r), 35.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_choose_e, (5,2, &r), 10.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_choose_e, (5,3, &r), 10.0 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_choose_e, (500,495, &r), 255244687600.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_choose_e, (500,5, &r), 255244687600.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_choose_e, (500,200, &r), 5.054949849935532221e+144 , TEST_TOL5, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_choose_e, (500,300, &r), 5.054949849935532221e+144 , TEST_TOL5, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_lnpoch_e, (5, 0.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnpoch_e, (5, 1.0/65536.0, &r), 0.000022981557571259389129, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnpoch_e, (5, 1.0/256.0, &r),   0.005884960217985189004,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnpoch_e, (7,3, &r), 6.222576268071368616, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_lnpoch_e, (5,2, &r), 3.401197381662155375, TEST_TOL0, GSL_SUCCESS);

  TEST_SF_SGN(s, gsl_sf_lnpoch_sgn_e, (5.0, 0.0, &r, &sgn), 0.0, TEST_TOL1, 1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lnpoch_sgn_e, (-4.5, 0.25, &r, &sgn), 0.7430116475119920117, TEST_TOL1, 1.0, GSL_SUCCESS);
  TEST_SF_SGN(s, gsl_sf_lnpoch_sgn_e, (-4.5, 1.25, &r, &sgn), 2.1899306304483174731, TEST_TOL1, -1.0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_poch_e, (5, 0.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (7,3, &r), 504.0 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (5,2, &r), 30.0  , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (5,1.0/256.0, &r), 1.0059023106151364982 , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_pochrel_e, (5,0, &r), 1.506117668431800472, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (7,3, &r), 503.0/3.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (5,2, &r), 29.0/2.0, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (5,0.01, &r), 1.5186393661368275330, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (-5.5,0.01, &r), 1.8584945633829063516, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (-5.5,-1.0/8.0, &r), 1.0883319303552135488, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (-5.5,-1.0/256.0, &r), 1.7678268037726177453, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_pochrel_e, (-5.5,-11.0, &r), 0.09090909090939652475, TEST_TOL0, GSL_SUCCESS);

  /* Add tests for special cases with negative arguments */
  TEST_SF(s,  gsl_sf_poch_e, (-9.0, -4.0, &r), 1.0/17160.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-9.0, -3.0, &r), -1.0/1320.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-9.0, -3.5, &r), 0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_poch_e, (-9.0, 4.0, &r), 3024.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-9.0, 3.0, &r), -504.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-9.0, 3.5, &r), 0.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_poch_e, (-9.0, 0.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_poch_e, (-8.0, -4.0, &r), 1.0/11880.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-8.0, -3.0, &r), -1.0/990.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-8.0, +4.0, &r), 1680.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-8.0, +3.0, &r), -336.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-3.0, +4.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);

  /* FIXME: we should be able to get an exact answer for poch(-a,a) if
     gsl_sf_lngamma functions were fixed to handle integer arguments
     exactly as a special case */
  TEST_SF(s,  gsl_sf_poch_e, (-3.0, +3.0, &r), -6.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_poch_e, (-4.0, +4.0, &r), 24.0, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_poch_e, (-3.0, +100.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1e-100, 0.001, &r), 1.0, TEST_TOL0, GSL_SUCCESS) ;
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (0.001, 0.001, &r), 0.9936876467088602902, TEST_TOL0, GSL_SUCCESS) ;
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (0.001, 1.0, &r), 0.9997803916424144436, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (0.001, 10.0, &r), 0.9999999958306921828, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1.0, 0.001, &r), 0.0009995001666250083319, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1.0, 1.01, &r), 0.6357810204284766802, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1.0, 10.0, &r), 0.9999546000702375151, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (10.0, 10.01, &r), 0.5433207586693410570, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (10.0, 20.0, &r), 0.9950045876916924128, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1000.0, 1000.1, &r), 0.5054666401440661753, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1000.0, 2000.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  /* Test for failure of the Gautschi recurrence (now fixed) for x = a - 2 */
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (34.0, 32.0, &r), 0.3849626436463866776322932129, TEST_TOL2, GSL_SUCCESS);
  /* and the next test is gamma_inc_P(37,35-20*eps) */
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (37.0, 3.499999999999999289e+01, &r), 0.3898035054195570860969333039, TEST_TOL2, GSL_SUCCESS);

  /* Regression test Martin Jansche <jansche@ling.ohio-state.edu> BUG#12 */
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (10, 1e-16, &r), 2.755731922398588814734648067e-167, TEST_TOL2, GSL_SUCCESS);

  /* Regression test for gsl_cdf_chisq_Pinv, (0.05, 1263131.0) */
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1263131.0, 1261282.3637, &r), 0.04994777516935182963821362168, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (1263131.0, 1263131.0, &r), 0.500118321758657770672882362502514254, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (0.0, 0.001, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (0.001, 0.001, &r), 0.006312353291139709793, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (0.001, 1.0, &r), 0.00021960835758555639171, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (0.001, 2.0, &r), 0.00004897691783098147880, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (0.001, 5.0, &r), 1.1509813397308608541e-06, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0, 0.001, &r), 0.9990004998333749917, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0, 1.01, &r), 0.3642189795715233198, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0, 10.0, &r), 0.00004539992976248485154, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (10.0, 10.01, &r), 0.4566792413306589430, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (10.0, 100.0, &r), 1.1253473960842733885e-31, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1000.0, 1000.1, &r), 0.4945333598559338247, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1000.0, 2000.0, &r), 6.847349459614753180e-136, TEST_TOL2, GSL_SUCCESS);

  /* designed to trap the a-x=1 problem */
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (100,  99.0, &r), 0.5266956696005394, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (200, 199.0, &r), 0.5188414119121281, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (100,  99.0, &r), 0.4733043303994607, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (200, 199.0, &r), 0.4811585880878718, TEST_TOL2, GSL_SUCCESS);

  /* Test for x86 cancellation problems */
  TEST_SF(s, gsl_sf_gamma_inc_P_e, (5670, 4574, &r),  3.063972328743934e-55, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (5670, 4574, &r), 1.0000000000000000, TEST_TOL2, GSL_SUCCESS);

  /* test suggested by Michel Lespinasse [gsl-discuss Sat, 13 Nov 2004] */
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0e+06-1.0, 1.0e+06-2.0, &r), 0.50026596175224547004, TEST_TOL3, GSL_SUCCESS);

  /* tests in asymptotic regime related to Lespinasse test */
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0e+06+2.0, 1.0e+06+1.0, &r), 0.50026596135330304336, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0e+06, 1.0e+06-2.0, &r), 0.50066490399940144811, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0e+07, 1.0e+07-2.0, &r), 0.50021026104978614908, TEST_TOL2, GSL_SUCCESS);

  /* non-normalized "Q" function */
  TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0/1048576.0, 1.0/1048576.0, &r), 13.285819596290624271, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 1.0/1048576.0, &r), 13.381275128625328858, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0,   1.0/1048576.0, &r), 1.0485617142715768655e+06, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.00001,0.001, &r), 6.3317681434563592142, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.0001,0.001, &r), 6.3338276439767189385, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 0.001, &r), 6.3544709102510843793, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.5,   0.001, &r), 59.763880515942196981, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0,   0.001, &r), 992.66896046923884234, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-3.5,   0.001, &r), 9.0224404490639003706e+09, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-10.5,  0.001, &r), 3.0083661558184815656e+30, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 0.1,   &r), 1.8249109609418620068, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.5,   0.1,   &r), 3.4017693366916154163, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-10.0,  0.1,   &r), 8.9490757483586989181e+08, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-10.5,  0.1,   &r), 2.6967403834226421766e+09, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 1.0,   &r), 0.21928612679072766340, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.5,   1.0,   &r), 0.17814771178156069019, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0,   1.0,   &r), 0.14849550677592204792, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-2.5,   1.0,   &r), 0.096556648631275160264, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0,   10.0,  &r), 3.8302404656316087616e-07, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 10.0,  &r), 4.1470562324807320961e-06, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-0.5,   10.0,  &r), 1.2609042613241570681e-06, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0,   10.0,  &r), 3.8302404656316087616e-07, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-10.5,  10.0,  &r), 6.8404927328441566785e-17, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-100.0, 10.0,  &r), 4.1238327669858313997e-107, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (-200.0, 10.0,  &r), 2.1614091830529343423e-207, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_gamma_inc_e, (  0.0,     0.001, &r), 6.3315393641361493320, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  0.001,   0.001, &r), 6.3087159394864007261, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  1.0,     0.001, &r), 0.99900049983337499167, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, ( 10.0,     0.001, &r), 362880.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  0.0,     1.0,   &r), 0.21938393439552027368, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  0.001,   1.0,   &r), 0.21948181320730279613, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  1.0,     1.0,   &r), 0.36787944117144232160, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, ( 10.0,     1.0,   &r), 362879.95956592242045, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (100.0,     1.0,   &r), 9.3326215443944152682e+155, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  0.0,   100.0, &r), 3.6835977616820321802e-46, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  0.001, 100.0, &r), 3.7006367674063550631e-46, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (  1.0,   100.0, &r), 3.7200759760208359630e-44, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, ( 10.0,   100.0, &r), 4.0836606309106112723e-26, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_gamma_inc_e, (100.0,   100.0, &r), 4.5421981208626694294e+155, TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_lnbeta_e, (1.0e-8, 1.0e-8, &r),  19.113827924512310617 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0e-8, 0.01, &r),  18.420681743788563403 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0e-8, 1.0, &r),  18.420680743952365472 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0e-8, 10.0, &r),  18.420680715662683009 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0e-8, 1000.0, &r),  18.420680669107656949 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (0.1, 0.1, &r), 2.9813614810376273949 , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (0.1, 1.0, &r),  2.3025850929940456840 , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (0.1, 100.0, &r),  1.7926462324527931217 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (0.1, 1000, &r),  1.5619821298353164928 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0, 1.00025, &r),  -0.0002499687552073570, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0, 1.01, &r),  -0.009950330853168082848 , TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (1.0, 1000.0, &r),  -6.907755278982137052 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (100.0, 100.0, &r),  -139.66525908670663927 , TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (100.0, 1000.0, &r),  -336.4348576477366051 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_lnbeta_e, (100.0, 1.0e+8, &r),  -1482.9339185256447309 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (1.0,   1.0, &r), 1.0                   , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (1.0, 1.001, &r), 0.9990009990009990010 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (1.0,   5.0, &r), 0.2                   , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (1.0,  100.0, &r), 0.01                  , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (10.0, 100.0, &r), 2.3455339739604649879e-15 , TEST_TOL2, GSL_SUCCESS);

  /* Test negative arguments */
  TEST_SF(s,  gsl_sf_beta_e, (2.5, -0.1, &r), -11.43621278354402041480, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (2.5, -1.1, &r), 14.555179906328753255202, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-0.25, -0.1, &r), -13.238937960945229110, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-1.25, -0.1, &r), -14.298052997820847439, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-100.1, -99.1, &r), -1.005181917797644630375787297e60, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-100.1, 99.3, &r), 0.0004474258199579694011200969001, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (100.1, -99.3, &r), 1.328660939628876472028853747, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-100.1, 1.2, &r), 0.00365530364287960795444856281, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (100.1, -1.2, &r), 1203.895236907821059270698160, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-100.1, -1.2, &r), -3236.073671884748847700283841, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (-100.001, 0.0099, &r), -853.946649365611147996495177, TEST_TOL4, GSL_SUCCESS);

  /* Other test cases */
  TEST_SF(s,  gsl_sf_beta_e, (1e-32, 1.5, &r), 1e32, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_e, (1e-6, 0.5, &r), 1000001.386293677092419390336, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_beta_e, (-1.5, 0.5, &r), 0.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_beta_inc_e, (1.0, 1.0, 0.0, &r), 0.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (1.0, 1.0, 1.0, &r), 1.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (0.1, 0.1, 1.0, &r), 1.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0,  1.0, 0.5, &r), 0.5, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 0.1,  1.0, 0.5, &r), 0.9330329915368074160, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (10.0,  1.0, 0.5, &r), 0.0009765625000000000000, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (50.0,  1.0, 0.5, &r), 8.881784197001252323e-16, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0,  0.1, 0.5, &r), 0.06696700846319258402, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0, 10.0, 0.5, &r), 0.99902343750000000000, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0, 50.0, 0.5, &r), 0.99999999999999911180, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0,  1.0, 0.1, &r), 0.10, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0,  2.0, 0.1, &r), 0.19, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 1.0,  2.0, 0.9, &r), 0.99, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (50.0, 60.0, 0.5, &r), 0.8309072939016694143, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (90.0, 90.0, 0.5, &r), 0.5, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, ( 500.0,  500.0, 0.6, &r), 0.9999999999157549630, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (5000.0, 5000.0, 0.4, &r), 4.518543727260666383e-91, TEST_TOL5, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (5000.0, 5000.0, 0.6, &r), 1.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (5000.0, 2000.0, 0.6, &r), 8.445388773903332659e-89, TEST_TOL5, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.1, -0.1, 1.0, &r), 1.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.1, -0.2, 1.0, &r), 1.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.2, -0.1, 1.0, &r), 1.0, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.1, -0.2, 0.5, &r), 0.675252001958389971991335, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.2, -0.1, 0.5, &r), 0.324747998041610028008665, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.1, -0.1, 0.0, &r), 0.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.1, -0.2, 0.0, &r), 0.0, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.2, -0.1, 0.0, &r), 0.0, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.1, -0.2, 0.3, &r), 0.7469186777964287252, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_beta_inc_e, (-0.2, -0.1, 0.3, &r), 0.3995299653262016818, TEST_TOL2, GSL_SUCCESS);

  /* Bug report from Thomas Tanner <tanner@gmx.de> */

  TEST_SF(s,  gsl_sf_beta_inc_e, (0.5, 101.5, 0.999457, &r), 1.0, TEST_TOL2, GSL_SUCCESS);

  return s;
}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
/**
 *
 * @file timing.c
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @author Dulceneia Becker
 * @date 2010-11-15
 *
 **/

 #define _FMULS FMULS_GEQRF( M, N )
#define _FADDS FADDS_GEQRF( M, N )

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef _OPENMP
#include <omp.h>
#endif

#ifdef PLASMA_EZTRACE
#include <eztrace.h>
#endif

#include <unistd.h>
#include <sys/time.h>
#include <sys/resource.h>

#ifdef USE_MKL
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#else
#include <cblas.h>
#include <lapacke.h>
#endif
#include "plasma.h"
#include "core_blas.h"
#include "flops.h"
#include "timing.h"
#include "auxiliary.h"
#include "main.h"
#include "workspace.h"
#include "Utils.h"
#define EPSILON 1.0E-9

static double RunTest(double *t_, struct user_parameters*);

static double
RunTest_dpotrf(real_Double_t *t_, struct user_parameters* params)
{
    double  t;
    int64_t N     = params->matrix_size;
    int64_t NB    = params->blocksize;
    int check     = params->check;
    int uplo = PlasmaUpper;
    double check_res = 0;

    /* Allocate Data */
    PLASMA_desc *descA = NULL;
    double* ptr = malloc(N * N * sizeof(double));
    PLASMA_Desc_Create(&descA, ptr, PlasmaRealDouble, NB, NB, NB*NB, N, N, 0, 0, N, N);

#pragma omp parallel
#pragma omp master
    plasma_pdplgsy_quark( (double)N, *descA, 51 );

    /* Save A for check */
    double *A = NULL;
    if(check) {
        A = (double*)malloc(N * N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descA, (void*)A, N);
    }

    double t_start, t_end;

    //START_TIMING();
    t_start = rtclock();
#pragma omp parallel
#pragma omp master
    plasma_pdpotrf_quark(uplo, *descA);
    //STOP_TIMING();
    t_end = rtclock();

    *t_ = t_end - t_start;

    /* Check the solution */
    if ( check )
    {
        PLASMA_desc *descB = NULL;
        double* ptr = (double*)malloc(N * sizeof(double));
        PLASMA_Desc_Create(&descB, ptr, PlasmaRealDouble, NB, NB, NB*NB, N, 1, 0, 0, N, 1);

        plasma_pdpltmg_seq(* descB, 7672 );
        double* B = (double*)malloc(N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descB, (void*)B, N);

        PLASMA_dpotrs_Tile( uplo, descA, descB );

        double* X = (double*)malloc(N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descB, (void*)X, N);

        check_res = d_check_solution(N, N, 1, A, N, B, X, N);

        PASTE_CODE_FREE_MATRIX( descB );
        free( A );
        free( B );
        free( X );
    }

    PASTE_CODE_FREE_MATRIX( descA );

    return check_res;
}

static double
RunTest_dgetrf(real_Double_t *t_, struct user_parameters* params)
{
    double  t;
    int64_t N     = params->matrix_size;
    int64_t NB    = params->blocksize;
    int check     = params->check;
    double check_res = 0;

    /* Allocate Data */
    PLASMA_desc *descA = NULL;
    double* ptr = malloc(N * N * sizeof(double));
    PLASMA_Desc_Create(&descA, ptr, PlasmaRealDouble, NB, NB, NB*NB, N, N, 0, 0, N, N);

    int* piv = (int*)malloc(N * sizeof(double));

#pragma omp parallel
#pragma omp master
    plasma_pdpltmg_quark(*descA, 3456);

    /* Save AT in lapack layout for check */
    double *A = NULL;
    if(check) {
        A = (double*)malloc(N * N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descA, (void*)A, N);
    }

    double t_start, t_end;

    //START_TIMING();
    t_start = rtclock();
#pragma omp parallel
#pragma omp master
    plasma_pdgetrf_rectil_quark(*descA, piv);
    //STOP_TIMING();
    t_end = rtclock();

    *t_ = t_end - t_start;

    /* Check the solution */
    if ( check )
    {
        PLASMA_desc *descB = NULL;
        double* ptr = (double*)malloc(N * sizeof(double));
        PLASMA_Desc_Create(&descB, ptr, PlasmaRealDouble, NB, NB, NB*NB, N, 1, 0, 0, N, 1);

        plasma_pdpltmg_seq(*descB, 7732 );
        double* b = (double*)malloc(N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descB, (void*)b, N);

        PLASMA_dgetrs_Tile( PlasmaNoTrans, descA, piv, descB );

        double* x = (double*)malloc(N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descB, (void*)x, N);

        check_res = d_check_solution(N, N, 1, A, N, b, x, N);

        PASTE_CODE_FREE_MATRIX( descB );
        free(A); free(b); free(x);
    }

    PASTE_CODE_FREE_MATRIX( descA );
    free( piv );

    return check_res;
}

static double
RunTest_dgeqrf(real_Double_t *t_, struct user_parameters* params)
{
    double t;
    PLASMA_desc *descT;
    int64_t N     = params->matrix_size;
    int64_t IB    = params->iblocksize;
    int64_t NB    = params->blocksize;
    int check     = params->check;
    double check_res = 0;

    /* Allocate Data */
    PLASMA_desc *descA = NULL;
    double *ptr = (double*)malloc(N * N * sizeof(double));
    PLASMA_Desc_Create(&descA, ptr, PlasmaRealDouble, NB, NB, NB*NB, N, N, 0, 0, N, N);

#pragma omp parallel
#pragma omp master
    plasma_pdpltmg_quark(*descA, 5373 );

    /* Save A for check */
    double *A = NULL;
    if ( check ) {
        A = (double*)malloc(N * N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descA, (void*)A, N);
    }

    /* Allocate Workspace */
    plasma_alloc_ibnb_tile(N, N, PlasmaRealDouble, &descT, IB, NB);

    double t_start, t_end;

    //START_TIMING();
    t_start = rtclock();
#pragma omp parallel
#pragma omp master
    plasma_pdgeqrf_quark( *descA, *descT , IB);
    //STOP_TIMING();
    t_end = rtclock();

    *t_ = t_end - t_start;

    /* Check the solution */
    if ( check )
    {
        /* Allocate B for check */
        PLASMA_desc *descB = NULL;
        double* ptr = (double*)malloc(N * sizeof(double));
        PLASMA_Desc_Create(&descB, ptr, PlasmaRealDouble, NB, NB, NB*NB, N, 1, 0, 0, N, 1);

        /* Initialize and save B */
        plasma_pdpltmg_seq(*descB, 2264 );
        double *B = (double*)malloc(N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descB, (void*)B, N);

        /* Compute the solution */
        PLASMA_dgeqrs_Tile( descA, descT, descB , IB);

        /* Copy solution to X */
        double *X = (double*)malloc(N * sizeof(double));
        plasma_pdtile_to_lapack_quark(*descB, (void*)X, N);

        check_res = d_check_solution(N, N, 1, A, N, B, X, N);

        /* Free checking structures */
        PASTE_CODE_FREE_MATRIX( descB );
        free( A );
        free( B );
        free( X );
    }

    /* Free data */
    PLASMA_Dealloc_Handle_Tile(&descT);
    PASTE_CODE_FREE_MATRIX( descA );

    return check_res;
}

int ISEED[4] = {0,0,0,1};   /* initial seed for zlarnv() */

double run(struct user_parameters* params)
{
    double   t;
    double   fmuls, fadds;
    double   flops;

    params->succeed = 1;
    int type = params->type;
    if (params->matrix_size <= 0) {
        params->matrix_size = 2048;
    }

    if (params->blocksize <= 0) {
        params->blocksize = 128;
    }
//ifdef
    if (params->iblocksize <= 0) {
        params->iblocksize = params->blocksize;
    }

    if (params->type <= 0) {
        params->type = 1;
	type = params->type;
    }

    int64_t N    = params->matrix_size;
    int64_t M    = params->matrix_size;
    fadds = (double)(_FADDS);
    fmuls = (double)(_FMULS);
    flops = 1e-9 * (fmuls + fadds);
    //mod aqui

    if (type == 1) {
      if (RunTest_dgeqrf(&t, params) > EPSILON && params->check)
          params->succeed = 0;
    } else if (type == 2) {
      if (RunTest_dgetrf(&t, params) > EPSILON && params->check)
          params->succeed = 0;
    } else if (type == 3) {
      if (RunTest_dpotrf(&t, params) > EPSILON && params->check)
          params->succeed = 0;
    }

    // return gflops
    return flops / t;
}

max_stars_repo_name: fabiodrg/Unibench

max_issues_repo_name: fabiodrg/Unibench

max_forks_repo_name: fabiodrg/Unibench


text:
#pragma once

#include <algorithm>
#include <gsl/gsl_assert>
#include "family_type.h"
#include "family_mask.h"
#include "entity_id.h"
#include "halley/data_structures/nullable_reference.h"
#include "halley/support/exception.h"
#include "halley/support/debug.h"
#include "halley/utils/utils.h"

namespace Halley {
	class Entity;
	class FamilyBindingBase;

	class Family {
		friend class World;

	public:
		explicit Family(FamilyMaskType inclusionMask, FamilyMaskType optionalMask);
		virtual ~Family() {}

		size_t count() const
		{
			return elemCount;
		}

		void* getElement(size_t n) const
		{
			return static_cast<char*>(elems) + (n * elemSize);
		}

		void addOnEntitiesAdded(FamilyBindingBase* bind);
		void removeOnEntityAdded(FamilyBindingBase* bind);
		void addOnEntitiesRemoved(FamilyBindingBase* bind);
		void removeOnEntityRemoved(FamilyBindingBase* bind);
		void addOnEntitiesReloaded(FamilyBindingBase* bind);
		void removeOnEntitiesReloaded(FamilyBindingBase* bind);

		void notifyAdd(void* entities, size_t count);
		void notifyRemove(void* entities, size_t count);
		void notifyReload(void* entities, size_t count);

	protected:
		virtual void addEntity(Entity& entity) = 0;
		virtual void refreshEntity(Entity& entity) = 0;
		void removeEntity(Entity& entity);
		void reloadEntity(Entity& entity);
		virtual void updateEntities() = 0;
		virtual void clearEntities() = 0;

		void* elems = nullptr;
		size_t elemCount = 0;
		size_t elemSize = 0;
		Vector<EntityId> toRemove;
		Vector<EntityId> toReload;

		Vector<FamilyBindingBase*> addEntityCallbacks;
		Vector<FamilyBindingBase*> removeEntityCallbacks;
		Vector<FamilyBindingBase*> modifiedEntityCallbacks;

	private:
		FamilyMaskType inclusionMask;
		FamilyMaskType optionalMask;
	};

	class FamilyBase {
	protected:
		NullableReferenceAnchor anchor;

	public:
		EntityId entityId;
	};

	template <typename T>
	class FamilyBaseOf : public FamilyBase {
	public:
		NullableReferenceOf<T> getReference()
		{
			return anchor.getReferenceOf<T>();
		}

		NullableReferenceOf<const T> getReference() const
		{
			return anchor.getReferenceOf<const T>();
		}
	};

	// Apple's Clang 3.5 does not seem to have constexpr std::max...
	constexpr size_t maxSize(size_t a, size_t b)
	{
		return a > b ? a : b;
	}

	template <typename T>
	class FamilyImpl : public Family
	{
		constexpr static size_t storageSize = sizeof(T) - alignUp(sizeof(FamilyBase), alignof(void*));
		static_assert(std::is_base_of<FamilyBase, T>::value, "Family type does not derive from FamilyBase");

		// I don't know why this needs to be aligned up to 8 on Win32. :|
		static_assert(alignUp(T::Type::getNumComponents() * sizeof(void*), size_t(8)) == storageSize, "Family type has unexpected storage size");

		struct StorageType : public FamilyBase
		{
			alignas(alignof(void*)) std::array<char, storageSize> data;
		};

	public:
		explicit FamilyImpl(MaskStorage& storage)
			: Family(T::Type::inclusionMask(storage), T::Type::optionalMask(storage))
		{
		}

	protected:
		void addEntity(Entity& entity) override
		{
			auto& e = entities.emplace_back();
			e.entityId = entity.getEntityId();
			T::Type::loadComponents(entity, &e.data[0]);

			dirty = true;
		}

		void refreshEntity(Entity& entity) override
		{
			for (auto& e: entities) {
				if (e.entityId == entity.getEntityId()) {
					T::Type::loadComponents(entity, &e.data[0]);
					break;
				}
			}
		}

		void updateEntities() override
		{
			if (dirty) {
				// Notify additions
				HALLEY_DEBUG_TRACE();
				size_t prevSize = elemCount;
				size_t curSize = entities.size();
				updateElems();
				Expects(curSize >= prevSize);
				if (curSize > prevSize) {
					notifyAdd(entities.data() + prevSize, curSize - prevSize);
				}

				dirty = false;
			}

			if (!toReload.empty()) {
				// Notify reloads
				HALLEY_DEBUG_TRACE();
				std::vector<StorageType*> reloadedEntities;
				for (auto& entity : entities) {
					if (std::find(toReload.begin(), toReload.end(), entity.entityId) != toReload.end()) {
						reloadedEntities.push_back(&entity);
					}
				}
				notifyReload(reloadedEntities.data(), reloadedEntities.size());
				toReload.clear();
			}

			// Remove
			removeDeadEntities();
		}

		void clearEntities() override
		{
			notifyRemove(entities.data(), entities.size());
			entities.clear();
			updateElems();
		}

	private:
		Vector<StorageType> entities;
		bool dirty = false;

		void updateElems()
		{
			elems = entities.empty() ? nullptr : entities.data();
			elemCount = entities.size();
			elemSize = sizeof(StorageType);
		}

		void removeDeadEntities()
		{
			// Performance-critical code
			// Benchmarks suggest that using a Vector is faster than std::set and std::unordered_set
			if (!toRemove.empty()) {
				HALLEY_DEBUG_TRACE();
				size_t removeCount = toRemove.size();
				Expects(removeCount > 0);
				Expects(removeCount <= entities.size());
				std::sort(toRemove.begin(), toRemove.end());

				for (size_t i = 1; i < toRemove.size(); ++i) {
					Expects(toRemove[i - 1] != toRemove[i]);
				}

				// Move all entities to be removed to the back of the vector
				{
					int n = int(entities.size());
					// Note: it's important to scan it forward. Scanning backwards would improve performance for short-lived entities,
					// but it causes an issue where an entity is removed and added to the same family in one frame.
					for (int i = 0; i < n; i++) {
						EntityId id = entities[i].entityId;
						auto iter = std::lower_bound(toRemove.begin(), toRemove.end(), id);
						if (iter != toRemove.end() && id == *iter) {
							toRemove.erase(iter);
							if (i != n - 1) {
								std::swap(entities[i], entities[n - 1]);
								i--;
							}
							n--;
							if (toRemove.empty()) {
								break;
							}
						}
					}
					Ensures(size_t(n) + removeCount == entities.size());
				}

				Expects(toRemove.empty());

				// Notify removal
				size_t newSize = entities.size() - removeCount;
				Ensures(newSize < entities.size());
				notifyRemove(entities.data() + newSize, removeCount);

				// Remove them
				entities.resize(newSize);
				updateElems();
			}
			Ensures(toRemove.empty());
		}
	};
}

max_stars_repo_name: code-disaster/halley

max_issues_repo_name: code-disaster/halley

max_forks_repo_name: code-disaster/halley


text:
/**
 * This file is part of the "libterminal" project
 *   Copyright (c) 2019-2020 Christian Parpart <christian@parpart.family>
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#pragma once

#include <text_shaper/font.h>

#include <fmt/format.h>

#include <gsl/span>
#include <gsl/span_ext>

#include <optional>
#include <variant>
#include <vector>

namespace text
{

/// Holds the system path to a font file.
struct font_path
{
    std::string value;
};

/// Holds a view into the contents of a font file.
struct font_memory_ref
{
    std::string identifier;  //!< a unique identifier for this font
    gsl::span<uint8_t> data; //!< font file contents (non-owned)
};

/// Represents a font source (such as file path or memory).
using font_source = std::variant<font_path, font_memory_ref>;

/// Holds a list of fonts.
using font_source_list = std::vector<font_source>;

/**
 * Font location API.
 *
 * Used for locating fonts and fallback fonts to be used
 * for text shaping and glyph rendering.
 */
class font_locator
{
  public:
    virtual ~font_locator() = default;

    /**
     * Enumerates all available fonts.
     */
    virtual font_source_list all() = 0;

    /**
     * Locates the font matching the given description the best
     * and an ordered list of fallback fonts.
     */
    virtual font_source_list locate(font_description const& description) = 0;

    /**
     * Resolves the given codepoint sequence into an ordered list of
     * possible fonts that can be used for text shaping the given
     * codepoint sequence.
     */
    virtual font_source_list resolve(gsl::span<const char32_t> codepoints) = 0;
};

} // namespace text

namespace fmt // {{{
{
template <>
struct formatter<text::font_path>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(text::font_path path, FormatContext& ctx)
    {
        return fmt::format_to(ctx.out(), "path {}", path.value);
    }
};

template <>
struct formatter<text::font_memory_ref>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(text::font_memory_ref ref, FormatContext& ctx)
    {
        return fmt::format_to(ctx.out(), "in-memory: {}", ref.identifier);
    }
};

template <>
struct formatter<text::font_source>
{
    template <typename ParseContext>
    constexpr auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(text::font_source source, FormatContext& ctx)
    {
        if (std::holds_alternative<text::font_path>(source))
            return fmt::format_to(ctx.out(), "{}", std::get<text::font_path>(source));
        if (std::holds_alternative<text::font_memory_ref>(source))
            return fmt::format_to(ctx.out(), "{}", std::get<text::font_memory_ref>(source));
        return fmt::format_to(ctx.out(), "UNKNOWN SOURCE");
    }
};
} // namespace fmt

max_stars_repo_name: christianparpart/libterminal

max_issues_repo_name: christianparpart/libterminal

max_forks_repo_name: christianparpart/libterminal


text:
/***************************************************************************
 *  *
 *  *   Copyright (C) 2011 by Stefan Oslowski
 *  *   Licensed under the Academic Free License version 2.1
 *  *
 ****************************************************************************/

#ifndef __Pulsar__TimeDomainCovariance_h
#define __Pulsar__TimeDomainCovariance_h

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

#include "Pulsar/ProfileCovariance.h"

namespace Pulsar
{
  class TimeDomainCovariance : public ProfileCovariance
  {
  public:
    //! Default constructor
    TimeDomainCovariance ();

    //! Destructor
    ~TimeDomainCovariance ();

    //! Choose which bins are to be used in the calculation of covariance matrix
    void choose_bins ( unsigned, unsigned );

    //! Get the name of class used
    std::string get_name () const;

    //! Adds the given Profile to the covariance matrix estimate
    void add_Profile (const Profile* );
    void add_Profile (const Profile*, float wt );

    //! Get the whole covariance matrix
    void get_covariance_matrix ( double * );

    //! Get the whole covariance matrix - gsl
    void get_covariance_matrix_gsl ( gsl_matrix * );

    //! Set the whole covariance matrix
    void set_covariance_matrix ( double * );

    //! Get the i, j covariance matrix element, where i is the row and j is the column
    double get_covariance_matrix_value ( unsigned, unsigned );

    //! This normalizes the covariance matrix after all profiles were added
    void finalize ();

    //! Return a text interface that can be used to configure this instance
    //virtual TextInterface::Parser* get_interface ();

    //! Get the count of profiles used to estimate the covariance matrix
    unsigned get_count ();

    //! Set the count of profiles used to estimate the covariance matrix
    void set_count ( unsigned );

    //! Get the rank of covariance matrix
    unsigned get_rank ();

    //! Set the rank of covariance matrix
    void set_rank ( unsigned );

    //! Reset the covariance matrix calculations
    void reset ();

  protected:
    //! first bin to be used
    unsigned first_bin;

    //! last bin to be used
    unsigned last_bin;

    //! Auxiliary variable for covariance matrix_calculation
    gsl_matrix *covariance_gsl;

    //! Auxiliary variable for storing the profile
    gsl_vector_view p_gsl;
  };
}

#endif


text:
/* specfunc/debye.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */
/* augmented to n=5 and 6 2005-11-08 by R. J. Mathar, http://www.strw.leidenuniv.nl/~mathar */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_debye.h>

#include "error.h"
#include "check.h"

#include "chebyshev.h"
#include "cheb_eval.c"

static double adeb1_data[17] = {
   2.4006597190381410194,
   0.1937213042189360089,
  -0.62329124554895770e-02,
   0.3511174770206480e-03,
  -0.228222466701231e-04,
   0.15805467875030e-05,
  -0.1135378197072e-06,
   0.83583361188e-08,
  -0.6264424787e-09,
   0.476033489e-10,
  -0.36574154e-11,
   0.2835431e-12,
  -0.221473e-13,
   0.17409e-14,
  -0.1376e-15,
   0.109e-16,
  -0.9e-18
};
static cheb_series adeb1_cs = {
  adeb1_data,
  16,
  -1.0, 1.0,
  9
};

static double adeb2_data[18] = {
   2.5943810232570770282,
   0.2863357204530719834,
  -0.102062656158046713e-01,
   0.6049109775346844e-03,
  -0.405257658950210e-04,
   0.28633826328811e-05,
  -0.2086394303065e-06,
   0.155237875826e-07,
  -0.11731280087e-08,
   0.897358589e-10,
  -0.69317614e-11,
   0.5398057e-12,
  -0.423241e-13,
   0.33378e-14,
  -0.2645e-15,
   0.211e-16,
  -0.17e-17,
   0.1e-18
};
static cheb_series adeb2_cs = {
  adeb2_data,
  17,
  -1.0, 1.0,
  10
};

static double adeb3_data[17] = {
   2.707737068327440945,
   0.340068135211091751,
  -0.12945150184440869e-01,
   0.7963755380173816e-03,
  -0.546360009590824e-04,
   0.39243019598805e-05,
  -0.2894032823539e-06,
   0.217317613962e-07,
  -0.16542099950e-08,
   0.1272796189e-09,
  -0.987963460e-11,
   0.7725074e-12,
  -0.607797e-13,
   0.48076e-14,
  -0.3820e-15,
   0.305e-16,
  -0.24e-17
};
static cheb_series adeb3_cs = {
  adeb3_data,
  16,
  -1.0, 1.0,
  10
};

static double adeb4_data[17] = {
   2.781869415020523460,
   0.374976783526892863,
  -0.14940907399031583e-01,
   0.945679811437042e-03,
  -0.66132916138933e-04,
   0.4815632982144e-05,
  -0.3588083958759e-06,
   0.271601187416e-07,
  -0.20807099122e-08,
   0.1609383869e-09,
  -0.125470979e-10,
   0.9847265e-12,
  -0.777237e-13,
   0.61648e-14,
  -0.4911e-15,
   0.393e-16,
  -0.32e-17
};
static cheb_series adeb4_cs = {
  adeb4_data,
  16,
  -1.0, 1.0,
  10
};

static double adeb5_data[17] = {
   2.8340269546834530149,
   0.3994098857106266445,
  -0.164566764773099646e-1,
   0.10652138340664541e-2,
  -0.756730374875418e-4,
   0.55745985240273e-5,
  -0.4190692330918e-6,
   0.319456143678e-7,
  -0.24613318171e-8,
   0.1912801633e-9,
  -0.149720049e-10,
   0.11790312e-11,
  -0.933329e-13,
   0.74218e-14,
  -0.5925e-15,
   0.475e-16,
  -0.39e-17
};
static cheb_series adeb5_cs = {
  adeb5_data,
  16,
  -1.0, 1.0,
  10
};

static double adeb6_data[17] = {
  2.8726727134130122113,
  0.4174375352339027746,
 -0.176453849354067873e-1,
  0.11629852733494556e-2,
 -0.837118027357117e-4,
  0.62283611596189e-5,
 -0.4718644465636e-6,
  0.361950397806e-7,
 -0.28030368010e-8,
  0.2187681983e-9,
 -0.171857387e-10,
  0.13575809e-11,
 -0.1077580e-12,
  0.85893e-14,
 -0.6872e-15,
  0.552e-16,
 -0.44e-17
};
static cheb_series adeb6_cs = {
  adeb6_data,
  16,
  -1.0, 1.0,
  10
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_debye_1_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 1.64493406684822644;
  const double xcut = -GSL_LOG_DBL_MIN;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - 0.25*x + x*x/36.0;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double t = x*x/8.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&adeb1_cs, t, &c);
    result->val = c.val - 0.25 * x;
    result->err = c.err + 0.25 * x * GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
    const int nexp = floor(xcut/x);
    const double ex  = exp(-x);
    double sum = 0.0;
    double xk  = nexp * x;
    double rk  = nexp;
    int i;
    for(i=nexp; i>=1; i--) {
      sum *= ex;
      sum += (1.0 + 1.0/xk)/rk;
      rk -= 1.0;
      xk -= x;
    }
    result->val = val_infinity/x - sum*ex;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < xcut) {
    result->val = (val_infinity - exp(-x)*(x+1.0)) / x;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = val_infinity/x;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_debye_2_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 4.80822761263837714;
  const double xcut = -GSL_LOG_DBL_MIN;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - x/3.0 + x*x/24.0;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double t = x*x/8.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&adeb2_cs, t, &c);
    result->val = c.val - x/3.0;
    result->err = c.err + GSL_DBL_EPSILON * x/3.0;
    return GSL_SUCCESS;
  }
  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
    const int nexp = floor(xcut/x);
    const double ex  = exp(-x);
    double xk  = nexp * x;
    double rk  = nexp;
    double sum = 0.0;
    int i;
    for(i=nexp; i>=1; i--) {
      sum *= ex;
      sum += (1.0 + 2.0/xk + 2.0/(xk*xk)) / rk;
      rk -= 1.0;
      xk -= x;
    }
    result->val = val_infinity/(x*x) - 2.0 * sum * ex;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < xcut) {
    const double x2  = x*x;
    const double sum = 2.0 + 2.0*x + x2;
    result->val = (val_infinity - 2.0 * sum * exp(-x)) / x2;
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    result->val = (val_infinity/x)/x;
    result->err = GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_debye_3_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 19.4818182068004875;
  const double xcut = -GSL_LOG_DBL_MIN;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - 3.0*x/8.0 + x*x/20.0;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double t = x*x/8.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&adeb3_cs, t, &c);
    result->val = c.val - 0.375*x;
    result->err = c.err + GSL_DBL_EPSILON * 0.375*x;
    return GSL_SUCCESS;
  }
  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
    const int nexp = floor(xcut/x);
    const double ex  = exp(-x);
    double xk  = nexp * x;
    double rk  = nexp;
    double sum = 0.0;
    int i;
    for(i=nexp; i>=1; i--) {
      double xk_inv = 1.0/xk;
      sum *= ex;
      sum += (((6.0*xk_inv + 6.0)*xk_inv + 3.0)*xk_inv + 1.0) / rk;
      rk -= 1.0;
      xk -= x;
    }
    result->val = val_infinity/(x*x*x) - 3.0 * sum * ex;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x < xcut) {
    const double x3 = x*x*x;
    const double sum = 6.0 + 6.0*x + 3.0*x*x + x3;
    result->val = (val_infinity - 3.0 * sum * exp(-x)) / x3;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    result->val = ((val_infinity/x)/x)/x;
    result->err = GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_debye_4_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 99.5450644937635129;
  const double xcut = -GSL_LOG_DBL_MIN;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - 2.0*x/5.0 + x*x/18.0;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double t = x*x/8.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&adeb4_cs, t, &c);
    result->val = c.val - 2.0*x/5.0;
    result->err = c.err + GSL_DBL_EPSILON * 2.0*x/5.0;
    return GSL_SUCCESS;
  }
  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
    const int nexp = floor(xcut/x);
    const double ex  = exp(-x);
    double xk  = nexp * x;
    double rk  = nexp;
    double sum = 0.0;
    int i;
    for(i=nexp; i>=1; i--) {
      double xk_inv = 1.0/xk;
      sum *= ex;
      sum += ((((24.0*xk_inv + 24.0)*xk_inv + 12.0)*xk_inv + 4.0)*xk_inv + 1.0) / rk;
      rk -= 1.0;
      xk -= x;
    }
    result->val = val_infinity/(x*x*x*x) - 4.0 * sum * ex;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x < xcut) {
    const double x2 = x*x;
    const double x4 = x2*x2;
    const double sum = 24.0 + 24.0*x + 12.0*x2 + 4.0*x2*x + x4;
    result->val = (val_infinity - 4.0 * sum * exp(-x)) / x4;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    result->val = (((val_infinity/x)/x)/x)/x;
    result->err = GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_debye_5_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 610.405837190669483828710757875 ;
  const double xcut = -GSL_LOG_DBL_MIN;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - 5.0*x/12.0 + 5.0*x*x/84.0;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double t = x*x/8.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&adeb5_cs, t, &c);
    result->val = c.val - 5.0*x/12.0;
    result->err = c.err + GSL_DBL_EPSILON * 5.0*x/12.0;
    return GSL_SUCCESS;
  }
  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
    const int nexp = floor(xcut/x);
    const double ex  = exp(-x);
    double xk  = nexp * x;
    double rk  = nexp;
    double sum = 0.0;
    int i;
    for(i=nexp; i>=1; i--) {
      double xk_inv = 1.0/xk;
      sum *= ex;
      sum += (((((120.0*xk_inv + 120.0)*xk_inv + 60.0)*xk_inv + 20.0)*xk_inv + 5.0)*xk_inv+ 1.0) / rk;
      rk -= 1.0;
      xk -= x;
    }
    result->val = val_infinity/(x*x*x*x*x) - 5.0 * sum * ex;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x < xcut) {
    const double x2 = x*x;
    const double x4 = x2*x2;
    const double x5 = x4*x;
    const double sum = 120.0 + 120.0*x + 60.0*x2 + 20.0*x2*x + 5.0*x4 + x5;
    result->val = (val_infinity - 5.0 * sum * exp(-x)) / x5;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    result->val = ((((val_infinity/x)/x)/x)/x)/x;
    result->err = GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

int gsl_sf_debye_6_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 4356.06887828990661194792541535 ;
  const double xcut = -GSL_LOG_DBL_MIN;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0 - 3.0*x/7.0 + x*x/16.0;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double t = x*x/8.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&adeb6_cs, t, &c);
    result->val = c.val - 3.0*x/7.0;
    result->err = c.err + GSL_DBL_EPSILON * 3.0*x/7.0;
    return GSL_SUCCESS;
  }
  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
    const int nexp = floor(xcut/x);
    const double ex  = exp(-x);
    double xk  = nexp * x;
    double rk  = nexp;
    double sum = 0.0;
    int i;
    for(i=nexp; i>=1; i--) {
      double xk_inv = 1.0/xk;
      sum *= ex;
      sum += ((((((720.0*xk_inv + 720.0)*xk_inv + 360.0)*xk_inv + 120.0)*xk_inv + 30.0)*xk_inv+ 6.0)*xk_inv+ 1.0) / rk;
      rk -= 1.0;
      xk -= x;
    }
    result->val = val_infinity/(x*x*x*x*x*x) - 6.0 * sum * ex;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else if(x < xcut) {
    const double x2 = x*x;
    const double x4 = x2*x2;
    const double x6 = x4*x2;
    const double sum = 720.0 + 720.0*x + 360.0*x2 + 120.0*x2*x + 30.0*x4 + 6.0*x4*x +x6 ;
    result->val = (val_infinity - 6.0 * sum * exp(-x)) / x6;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    result->val = (((((val_infinity/x)/x)/x)/x)/x)/x ;
    result->err = GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_debye_1(const double x)
{
  EVAL_RESULT(gsl_sf_debye_1_e(x, &result));
}

double gsl_sf_debye_2(const double x)
{
  EVAL_RESULT(gsl_sf_debye_2_e(x, &result));
}

double gsl_sf_debye_3(const double x)
{
  EVAL_RESULT(gsl_sf_debye_3_e(x, &result));
}

double gsl_sf_debye_4(const double x)
{
  EVAL_RESULT(gsl_sf_debye_4_e(x, &result));
}

double gsl_sf_debye_5(const double x)
{
  EVAL_RESULT(gsl_sf_debye_5_e(x, &result));
}

double gsl_sf_debye_6(const double x)
{
  EVAL_RESULT(gsl_sf_debye_6_e(x, &result));
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once

#include "BoundingVolume.h"
#include "Library.h"
#include "TileContentLoadResult.h"
#include "TileContentLoader.h"
#include "TileID.h"
#include "TileRefine.h"

#include <CesiumAsync/AsyncSystem.h>
#include <CesiumAsync/Future.h>
#include <CesiumAsync/HttpHeaders.h>
#include <CesiumAsync/IAssetAccessor.h>
#include <CesiumGltfReader/GltfReader.h>

#include <glm/mat4x4.hpp>
#include <gsl/span>
#include <spdlog/fwd.h>

#include <cstddef>
#include <optional>

namespace CesiumGeospatial {
class GlobeRectangle;
}

namespace Cesium3DTilesSelection {

class Tileset;

/**
 * @brief Creates {@link TileContentLoadResult} from glTF data.
 */
class CESIUM3DTILESSELECTION_API GltfContent final : public TileContentLoader {
public:
  /**
   * @copydoc TileContentLoader::load
   *
   * The result will only contain the `model`. Other fields will be
   * empty or have default values.
   */
  CesiumAsync::Future<std::unique_ptr<TileContentLoadResult>>
  load(const TileContentLoadInput& input) override;

  /**
   * @brief Create a {@link TileContentLoadResult} from the given data.
   *
   * (Only public to be called from `Batched3DModelContent`)
   *
   * @param asyncSystem The async system to use for requesting any external
   * content.
   * @param pLogger Only used for logging
   * @param url The URL, only used for logging
   * @param headers The http headers to use for resolving any external content.
   * @param pAssetAccessor The asset accessor to use to resolve external
   * content.
   * @param data The actual glTF data
   * @param contentOptions The content options for loading this glTF.
   * @return The {@link TileContentLoadResult}
   */
  static CesiumAsync::Future<std::unique_ptr<TileContentLoadResult>> load(
      const CesiumAsync::AsyncSystem& asyncSystem,
      const std::shared_ptr<spdlog::logger>& pLogger,
      const std::string& url,
      const CesiumAsync::HttpHeaders& headers,
      const std::shared_ptr<CesiumAsync::IAssetAccessor>& pAssetAccessor,
      const gsl::span<const std::byte>& data,
      const TilesetContentOptions& contentOptions);

  /**
   * @brief Creates texture coordinates for mapping {@link RasterOverlay} tiles
   * to {@link Tileset} tiles.
   *
   * Generates new texture coordinates for the `gltf` using the given
   * `projections`. The first new texture coordinate (`u` or `s`) will be 0.0 at
   * the `minimumX` of the given `rectangle` and 1.0 at the `maximumX`. The
   * second texture coordinate (`v` or `t`) will be 0.0 at the `minimumY` of
   * the given `rectangle` and 1.0 at the `maximumY`.
   *
   * Coordinate values for vertices in between these extremes are determined by
   * projecting the vertex position with the `projection` and then computing the
   * fractional distance of that projected position between the minimum and
   * maximum.
   *
   * Projected positions that fall outside the `globeRectangle` will be clamped
   * to the edges, so the coordinate values will never be less then 0.0 or
   * greater than 1.0.
   *
   * These texture coordinates are stored in the provided glTF, and a new
   * primitive attribute named `_CESIUMOVERLAY_n` is added to each primitive,
   * where `n` starts with the `firstTextureCoordinateID` passed to this
   * function and increases with each projection.
   *
   * @param gltf The glTF model.
   * @param modelToEcefTransform The transformation of this glTF to ECEF
   * coordinates.
   * @param firstTextureCoordinateID The texture coordinate ID of the first
   * projection.
   * @param globeRectangle The rectangle that all projected vertex positions are
   * expected to lie within. If this parameter is std::nullopt, it is computed
   * from the vertices.
   * @param projections The projections for which to generate texture
   * coordinates. There is a linear relationship between the coordinates of this
   * projection and the generated texture coordinates.
   * @return The detailed of the generated texture coordinates.
   */
  static std::optional<TileContentDetailsForOverlays>
  createRasterOverlayTextureCoordinates(
      CesiumGltf::Model& gltf,
      const glm::dmat4& modelToEcefTransform,
      int32_t firstTextureCoordinateID,
      const std::optional<CesiumGeospatial::GlobeRectangle>& globeRectangle,
      std::vector<CesiumGeospatial::Projection>&& projections);

  /**
   * @brief Computes a bounding region from the vertex positions in a glTF
   * model.
   *
   * If the glTF model spans the anti-meridian, the west and east longitude
   * values will be in the usual -PI to PI range, but east will have a smaller
   * value than west.
   *
   * @param gltf The model.
   * @param transform The transform from model coordinates to ECEF coordinates.
   * @return The computed bounding region.
   */
  static CesiumGeospatial::BoundingRegion computeBoundingRegion(
      const CesiumGltf::Model& gltf,
      const glm::dmat4& transform);

  /**
   * @brief Applies the glTF's RTC_CENTER, if any, to the given rootTransform.
   *
   * @param gltf
   * @param rootTransform
   * @return glm::dmat4x4
   */

  /**
   * @brief Applies the glTF's RTC_CENTER, if any, to the given transform.
   *
   * If the glTF has a `CESIUM_RTC` extension, this function will multiply the
   * given matrix with the (translation) matrix that is created from the
   * `RTC_CENTER` in the. If the given model does not have this extension, then
   * this function will return the `rootTransform` unchanged.
   *
   * @param model The glTF model
   * @param rootTransform The matrix that will be multiplied with the transform
   * @return The result of multiplying the `RTC_CENTER` with the
   * `rootTransform`.
   */
  static glm::dmat4x4 applyRtcCenter(
      const CesiumGltf::Model& gltf,
      const glm::dmat4x4& rootTransform);

  /**
   * @brief Applies the glTF's `gltfUpAxis`, if any, to the given transform.
   *
   * By default, the up-axis of a glTF model will the the Y-axis.
   *
   * If the tileset that contained the model had the `asset.gltfUpAxis` string
   * property, then the information about the up-axis has been stored in as a
   * number property called `gltfUpAxis` in the `extras` of the given model.
   *
   * Depending on whether this value is `CesiumGeometry::Axis::X`, `Y`, or `Z`,
   * the given matrix will be multiplied with a matrix that converts the
   * respective axis to be the Z-axis, as required by the 3D Tiles standard.
   *
   * @param model The glTF model
   * @param rootTransform The matrix that will be multiplied with the transform
   * @return The result of multiplying the `rootTransform` with the
   * `gltfUpAxis`.
   */
  static glm::dmat4x4 applyGltfUpAxisTransform(
      const CesiumGltf::Model& model,
      const glm::dmat4x4& rootTransform);

private:
  static CesiumGltfReader::GltfReader _gltfReader;
};

} // namespace Cesium3DTilesSelection

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_roots.h>

#include "cosmocalc.h"

static double gaussiannorm_linear_powspec_exact_lnk_integ_funct(double lnk, void *p);
static double onederiv_gaussiannorm_linear_powspec_exact_lnk_integ_funct(double lnk, void *p);
static double twoderiv_gaussiannorm_linear_powspec_exact_lnk_integ_funct(double lnk, void *p);
static double nonlinear_gaussnorm_scale_funct(double gaussR, void *p);

static double gaussiannorm_linear_powspec_exact_lnk_integ_funct(double lnk, void *p)
{
  double gaussRad = (*(double*)p);
  double k = exp(lnk);
  return linear_powspec(k,1.0)*k*k*k/2.0/M_PI/M_PI*exp(-1.0*k*k*gaussRad*gaussRad);
}

static double onederiv_gaussiannorm_linear_powspec_exact_lnk_integ_funct(double lnk, void *p)
{
  double gaussRad = (*(double*)p);
  double k = exp(lnk);
  return linear_powspec(k,1.0)*k*k*k/2.0/M_PI/M_PI*exp(-1.0*k*k*gaussRad*gaussRad)*(-1.0*k*k*2.0*gaussRad);
}

static double twoderiv_gaussiannorm_linear_powspec_exact_lnk_integ_funct(double lnk, void *p)
{
  double gaussRad = (*(double*)p);
  double k = exp(lnk);
  return  linear_powspec(k,1.0)*k*k*k/2.0/M_PI/M_PI*exp(-1.0*k*k*gaussRad*gaussRad)*(-2.0*k*k + 4.0*k*k*k*k*gaussRad*gaussRad);
}

//uses Takahashi et al. (2012) arXiv:1208.2701 unless macro below is set to use true Smith+03
//#define SMITH03

double nonlinear_powspec(double k, double a)
{
  static int initFlag = 1;
  static int currCosmoNum;
  static gsl_spline *spline[4];
  static gsl_interp_accel *accel[4];
  int i;
  double xtab[COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH],ytab[COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH];
  double Rsigma,C,neff,ksigma,sigma2;
  double an,bn,cn,alphan,gamman,betan,mun,nun;
  double f1,f2,f3;
  double DeltakNL,dsigma2dR,d2sigma2d2R,PkNL,PkL;
  double y,DeltakL,fy,DeltakQ,DeltakHprime,DeltakH;
  //double t;
  double I0,I1;
  double abserr;
  gsl_integration_workspace *workspace;
  gsl_function F;
  double gaussRad;

#define WORKSPACE_NUM 10000000
#define ABSERR 1e-6
#define RELERR 0.0

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      currCosmoNum = cosmoData.cosmoNum;

      if(initFlag)
	{
	  for(i=0;i<4;++i)
	    spline[i] = gsl_spline_alloc(gsl_interp_akima,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));
	  for(i=0;i<4;++i)
	    accel[i] = gsl_interp_accel_alloc();

	  initFlag = 0;
	}
      else
	{
	  for(i=0;i<4;++i)
	    gsl_spline_free(spline[i]);
	  for(i=0;i<4;++i)
	    spline[i] = gsl_spline_alloc(gsl_interp_akima,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));
	  for(i=0;i<4;++i)
	    gsl_interp_accel_reset(accel[i]);
	}

      //t = -wtime();

      for(i=0;i<COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH;++i)
	{
	  xtab[i] = i*(1.0-AEXPN_MIN_NONLINPK)/(COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH-1.0) + AEXPN_MIN_NONLINPK;
	  ytab[i] = get_nonlinear_gaussnorm_scale(xtab[i]);
	}
      gsl_spline_init(spline[0],xtab,ytab,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));

      gsl_sort(ytab,(size_t) 1,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));
      for(i=0;i<COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH;++i)
	{
	  xtab[i] = ytab[i];
	  ytab[i] = gaussiannorm_linear_powspec(xtab[i]);
	}
      gsl_spline_init(spline[1],xtab,ytab,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));

      workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);
      F.params = &gaussRad;

      F.function = &onederiv_gaussiannorm_linear_powspec_exact_lnk_integ_funct;
      for(i=0;i<COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH;++i)
	{
	  gaussRad = xtab[i];
	  gsl_integration_qags(&F,log(1e-4),log(2.0*M_PI/gaussRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
	  gsl_integration_qags(&F,log(2.0*M_PI/gaussRad),log(1e3),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);
	  ytab[i] = I0+I1;
	}
      gsl_spline_init(spline[2],xtab,ytab,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));

      F.function = &twoderiv_gaussiannorm_linear_powspec_exact_lnk_integ_funct;
      for(i=0;i<COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH;++i)
	{
	  gaussRad = xtab[i];
	  gsl_integration_qags(&F,log(1e-4),log(2.0*M_PI/gaussRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
	  gsl_integration_qags(&F,log(2.0*M_PI/gaussRad),log(1e3),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);
	  ytab[i] = I0+I1;
	}
      gsl_spline_init(spline[3],xtab,ytab,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));

      gsl_integration_workspace_free(workspace);

#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM

      //t += wtime();
      //fprintf(stderr,"comp of non-linear Pk took %f seconds.\n",t);
    }

  Rsigma = gsl_spline_eval(spline[0],a,accel[0]);
  sigma2 = gsl_spline_eval(spline[1],Rsigma,accel[1]);
  dsigma2dR = gsl_spline_eval(spline[2],Rsigma,accel[2]);
  d2sigma2d2R = gsl_spline_eval(spline[3],Rsigma,accel[3]);

  ksigma = 1.0/Rsigma;
  neff = -1.0*Rsigma/sigma2*dsigma2dR - 3.0;
  C = -1.0*(d2sigma2d2R*Rsigma*Rsigma/sigma2 + dsigma2dR*Rsigma/sigma2 - dsigma2dR*dsigma2dR*Rsigma*Rsigma/sigma2/sigma2);

  double ha,weffa,omegaMz,omegaDEwz;

  ha = hubble_noscale(a);
  weffa = weff(a);
  omegaMz = cosmoData.OmegaM/a/a/a/ha/ha;
  omegaDEwz = (1.0-cosmoData.OmegaM)/ha/ha/pow(a,3.0*(1.0 + weffa));

#ifdef SMITH03
  an = pow(10.0,1.4861 + 1.8369*neff + 1.6762*neff*neff + 0.7940*neff*neff*neff + 0.1670*neff*neff*neff*neff - 0.6206*C);
  bn = pow(10.0,0.9463 + 0.9466*neff + 0.3084*neff*neff - 0.9400*C);
  cn = pow(10.0,-0.2807 + 0.6669*neff + 0.3214*neff*neff - 0.0793*C);
  gamman = 0.8649 + 0.2989*neff + 0.1631*C;
  alphan = 1.3884 + 0.3700*neff - 0.1452*neff*neff;
  betan = 0.8291 + 0.9854*neff + 0.3401*neff*neff;
  mun = pow(10.0,-3.5442 + 0.1908*neff);
  nun = pow(10.0,0.9589 + 1.2857*neff);

  f1 = pow(OmegaMz,-0.0307);
  f2 = pow(OmegaMz,-0.0585);
  f3 = pow(OmegaMz,0.0743);
#else
  an = pow(10.0,1.5222 + 2.8553*neff + 2.3706*neff*neff + 0.9903*neff*neff*neff + 0.2250*neff*neff*neff*neff - 0.6038*C + 0.1749*omegaDEwz*(1.0 + weffa));
  bn = pow(10.0,-0.5642 + 0.5864*neff + 0.5716*neff*neff - 1.5474*C + 0.2279*omegaDEwz*(1.0 + weffa));
  cn = pow(10.0,0.3698 + 2.0404*neff + 0.8161*neff*neff + 0.5869*C);
  gamman = 0.1971 - 0.0843*neff + 0.8460*C;
  alphan = fabs(6.0835 + 1.3373*neff - 0.1959*neff*neff - 5.5274*C);
  betan = 2.0379 - 0.7354*neff + 0.3157*neff*neff + 1.2490*neff*neff*neff + 0.3980*neff*neff*neff*neff - 0.1682*C;
  mun = 0.0;
  nun = pow(10.0,5.2105 + 3.6902*neff);

  f1 = pow(omegaMz,-0.0307);
  f2 = pow(omegaMz,-0.0585);
  f3 = pow(omegaMz,0.0743);
#endif

  PkL = linear_powspec(k,a);
  y = k/ksigma;
  fy = y/4.0 + y*y/8.0;
  DeltakL = PkL*k*k*k/2.0/M_PI/M_PI;

  DeltakQ = DeltakL*pow(1.0 + DeltakL,betan)/(1.0 + alphan*DeltakL)*exp(-1.0*fy);

  DeltakHprime = an*pow(y,3.0*f1)/(1.0 + bn*pow(y,f2) + pow(cn*f3*y,3.0 - gamman));
  DeltakH = DeltakHprime/(1.0 + mun/y + nun/y/y);

  DeltakNL = DeltakQ + DeltakH;
  PkNL = DeltakNL/(k*k*k/2.0/M_PI/M_PI);

  return PkNL;
}

double gaussiannorm_linear_powspec(double gaussRad)
{
#define PNL_RGAUSS_MIN 0.0001
#define PNL_RGAUSS_MAX 100.0
#define WORKSPACE_NUM 10000000
#define ABSERR 1e-6
#define RELERR 0.0

  static int initFlag = 1;
  static int currCosmoNum;
  static gsl_spline *spline = NULL;
  static gsl_interp_accel *accel = NULL;

  double I0,I1;
  double abserr;
  gsl_integration_workspace *workspace;
  gsl_function F;
  int i;
  double xtab[COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH];
  double ytab[COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH];
  double dlnr;
  double lnrmin;
  double lnr;

  if(initFlag == 1 || currCosmoNum != cosmoData.cosmoNum)
    {
      initFlag = 0;
      currCosmoNum = cosmoData.cosmoNum;

      workspace = gsl_integration_workspace_alloc((size_t) WORKSPACE_NUM);
      F.function = &gaussiannorm_linear_powspec_exact_lnk_integ_funct;

      dlnr = log(PNL_RGAUSS_MAX/PNL_RGAUSS_MIN)/(COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH-1.0);
      lnrmin = log(PNL_RGAUSS_MIN);

      for(i=0;i<COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH;++i)
	{
	  lnr = dlnr*i + lnrmin;
	  xtab[i] = exp(lnr);
	  F.params = &(xtab[i]);

	  gsl_integration_qags(&F,log(1e-4),log(2.0*M_PI/gaussRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
	  gsl_integration_qags(&F,log(2.0*M_PI/gaussRad),log(1e3),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);

	  //gsl_integration_qagil(&F,log(2.0*M_PI/gaussRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I0,&abserr);
	  //gsl_integration_qagiu(&F,log(2.0*M_PI/gaussRad),ABSERR,RELERR,(size_t) WORKSPACE_NUM,workspace,&I1,&abserr);

	  xtab[i] = lnr;
	  ytab[i] = log(I0+I1);
	}

      gsl_integration_workspace_free(workspace);

      if(spline != NULL)
	gsl_spline_free(spline);
      spline = gsl_spline_alloc(GSL_SPLINE_TYPE,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));
      gsl_spline_init(spline,xtab,ytab,(size_t) (COSMOCALC_NONLINEAR_POWSPEC_TABLE_LENGTH));
      if(accel != NULL)
	gsl_interp_accel_reset(accel);
      else
	accel = gsl_interp_accel_alloc();

#undef ABSERR
#undef RELERR
#undef WORKSPACE_NUM
#undef PNL_RGAUSS_MIN
#undef PNL_RGAUSS_MAX
    }

  return exp(gsl_spline_eval(spline,log(gaussRad),accel));
}

static double nonlinear_gaussnorm_scale_funct(double gaussR, void *p)
{
  double gf = ((double*)p)[0];

  return gaussiannorm_linear_powspec(gaussR)*gf*gf-1.0;
}

double get_nonlinear_gaussnorm_scale(double a)
{
  double gf = growth_function(a);
  double Rsigma,Rlow=0.001,Rhigh=10.0;
  int itr,maxItr=1000,status;

#define ABSERR 1e-6
#define RELERR 1e-6
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  gsl_function F;

  F.function = &nonlinear_gaussnorm_scale_funct;
  F.params = &gf;

  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc(T);
  gsl_root_fsolver_set(s,&F,Rlow,Rhigh);
  itr = 0;

  do
    {
      itr++;
      status = gsl_root_fsolver_iterate(s);
      Rsigma = gsl_root_fsolver_root(s);
      Rlow = gsl_root_fsolver_x_lower(s);
      Rhigh = gsl_root_fsolver_x_upper(s);
      status = gsl_root_test_interval(Rlow,Rhigh,ABSERR,RELERR);
    }
  while(status == GSL_CONTINUE && itr < maxItr);

#undef ABSERR
#undef RELERR

  gsl_root_fsolver_free(s);

  return Rsigma;
}

max_stars_repo_name: beckermr/cosmocalc

max_issues_repo_name: beckermr/cosmocalc

max_forks_repo_name: beckermr/cosmocalc


text: #include <cblas.h>


text:
#pragma once
#include <fstream>
#include <memory>
#include <optional>
#include <sstream>
#include <utility>

#include <gsl/pointers>

#include <fmt/core.h>
#include <fmt/ranges.h>

#include <nop/serializer.h>
#include <nop/status.h>
#include <nop/structure.h>
#include <nop/utility/die.h>
#include <nop/utility/stream_reader.h>
#include <nop/utility/stream_writer.h>

#include <core/Util.h>
#include <core/storage/Page.h>
#include <core/storage/PageCache.h>
#include <core/storage/btree/Config.h>
#include <core/storage/btree/Node.h>

namespace internal::storage::btree {

namespace util {
template<BtreeConfig Config = DefaultConfig>
class BtreePrinter;
}

template<BtreeConfig Config = DefaultConfig>
class Btree final {
	using Self = Btree<Config>;

	using Key = typename Config::Key;
	using Val = typename Config::Val;
	using Ref = typename Config::Ref;
	using Nod = Node<Config>;

	friend util::BtreePrinter<Config>;

public:
	//! Number of entries in branch and leaf nodes may differ
	//! Directly unwrap with `.value()` since we _want to fail at compile time_ in case their is no value which
	//! satisfies the predicates
	int NUM_LINKS_BRANCH = ::internal::binsearch_primitive(2l, Page::size(), [](long current, long, long)
			{ return nop::Encoding<Nod>::Size({typename Nod::Metadata(typename Nod::Branch(std::vector<Ref>(current), std::vector<Position>(current))), {}, true}) - Page::size(); }).value();
	int NUM_RECORDS_BRANCH = NUM_LINKS_BRANCH - 1;

	//! Equivalent to `m` in Knuth's definition
	//! Make sure that when a leaf is split, its contents could be distributed among the two branch nodes.
	//! Directly unwrap with `.value()` since we _want to fail at compile time_ in case their is no value which
	//! satisfies the predicates
	int _NUM_RECORDS_LEAF = ::internal::binsearch_primitive(2l, Page::size(), [](long current, long, long)
			{ return nop::Encoding<Nod>::Size({typename Nod::Metadata(typename Nod::Leaf(std::vector<Key>(current), std::vector<Val>(current))), {}, true}) - Page::size(); }).value();
	int NUM_RECORDS_LEAF = _NUM_RECORDS_LEAF - 1 >= NUM_RECORDS_BRANCH * 2
	        ? NUM_RECORDS_BRANCH * 2 - 1
	        : _NUM_RECORDS_LEAF;

private:
	static inline constexpr bool APPLY_COMPRESSION = Config::APPLY_COMPRESSION;
	static inline constexpr int PAGE_CACHE_SIZE = Config::PAGE_CACHE_SIZE;

	static inline constexpr uint32_t MAGIC = 0xB75EEA41;

public:
	struct Header {
		Position m_rootpos;
		std::size_t m_size{};
		std::size_t m_depth{};
		uint32_t m_magic{MAGIC};
		uint32_t m_pgcache_size{PAGE_CACHE_SIZE};
		uint8_t m_apply_compression{APPLY_COMPRESSION};

		/*
		 * The storage locations of the tree header and the tree contents differ.
		 * This one stores the name of the file which contains the header of the tree,
		 * whereas the other name passed to the Btree constructor is the one where
		 * the actual nodes of the tree are stored.
		 */
		std::string m_content_file;

		/*
		 * We don't want to serialize this. It is used only to check whether the header
		 * we are currently possessing is valid.
		 */
		bool m_dirty{false};

		Header() = default;
		explicit Header(Position rootpos, std::size_t size, std::size_t depth, std::string_view content_file)
		    : m_rootpos{rootpos},
		      m_size{size},
		      m_depth{depth},
		      m_content_file{std::string{content_file}} {}

		[[nodiscard]] auto &rootpos() noexcept { return m_rootpos; }

		[[nodiscard]] auto &size() noexcept { return m_size; }

		[[nodiscard]] auto &depth() noexcept { return m_depth; }

		[[nodiscard]] auto &dirty() noexcept { return m_dirty; }

		auto operator<=>(const Header &) const noexcept = default;

		friend std::ostream &operator<<(std::ostream &os, const Header &h) {
			os << "Header { .rootpos = " << h.m_rootpos << ", .size =" << h.m_size << ", .depth =" << h.m_depth << " }";
			return os;
		}

		NOP_STRUCTURE(Header, m_rootpos, m_size, m_depth, m_magic, m_pgcache_size, m_apply_compression, m_content_file);
	};

private:
	[[nodiscard]] bool is_node_full(const Self::Nod &node) {
		if (node.is_branch())
			return node.is_full(NUM_RECORDS_BRANCH);
		return node.is_full(NUM_RECORDS_LEAF);
	}

	[[nodiscard]] auto node_split(Self::Nod &node) {
		if (node.is_branch())
			return node.split(NUM_RECORDS_BRANCH);
		return node.split(NUM_RECORDS_LEAF);
	}

	[[nodiscard]] std::optional<Val> search_subtree(const Self::Nod &node, const Self::Key &target_key) const noexcept {
		if (node.is_branch()) {
			const auto &refs = node.branch().m_refs;
			const std::size_t index = std::lower_bound(refs.cbegin(), refs.cend(), target_key) - refs.cbegin();
			const Position pos = node.branch().m_links[index];
			const auto other = Nod::from_page(m_pgcache.get_page(pos));
			return search_subtree(other, target_key);
		}
		assert(node.is_leaf());

		const auto &keys = node.leaf().m_keys;
		const auto &vals = node.leaf().m_vals;
		const auto it = std::lower_bound(keys.cbegin(), keys.cend(), target_key);
		if (it == keys.cend() || *it != target_key)
			return {};
		return vals[it - keys.cbegin()];
	}

	Nod make_new_root() {
		auto old_root = root();
		auto old_pos = m_rootpos;

		auto new_pos = m_pgcache.get_new_pos();

		old_root.set_parent(new_pos);
		old_root.set_root(false);

		auto [midkey, sibling] = node_split(old_root);
		auto sibling_pos = m_pgcache.get_new_pos();

		Nod new_root{typename Nod::Metadata(typename Nod::Branch({midkey}, {old_pos, sibling_pos})), new_pos, true};

		m_pgcache.put_page(new_pos, new_root.make_page());
		m_pgcache.put_page(old_pos, old_root.make_page());
		m_pgcache.put_page(sibling_pos, sibling.make_page());

		m_rootpos = new_pos;
		++m_depth;
		m_header.m_dirty = true;

		return new_root;
	}

public:
	[[nodiscard]] auto root() noexcept { return std::as_const(*this).root(); }

	[[nodiscard]] auto root() const noexcept { return Nod::from_page(m_pgcache.get_page(rootpos())); }

	[[nodiscard]] auto &header() noexcept { return m_header; }

	[[nodiscard]] const auto &rootpos() const noexcept { return m_rootpos; }

	[[nodiscard]] std::size_t size() noexcept { return m_size; }

	[[nodiscard]] std::size_t size() const noexcept { return m_size; }

	[[nodiscard]] bool empty() const noexcept { return size() == 0; }

	[[nodiscard]] bool empty() noexcept { return size() == 0; }

	[[nodiscard]] std::size_t depth() { return m_depth; }

	[[nodiscard]] auto num_records_leaf() const noexcept { return NUM_RECORDS_LEAF; }

	[[nodiscard]] auto num_records_branch() const noexcept { return NUM_RECORDS_BRANCH; }

public:
	/*
	 *  Operations API
	 */

	void put(const Self::Key &key, const Self::Val &val) {
		Position currpos{rootpos()};
		Nod curr{root()};

		if (is_node_full(curr))
			curr = make_new_root();

		while (true) {
			if (curr.is_leaf()) {
				/* fmt::print(" -- Putting kv pair in leaf ... \n"); */
				auto &keys = curr.leaf().m_keys;
				auto &vals = curr.leaf().m_vals;

				const std::size_t index = std::lower_bound(keys.cbegin(), keys.cend(), key) - keys.cbegin();
				if (!keys.empty() && index < keys.size() && keys[index] == key)
					return;

				keys.insert(keys.begin() + index, key);
				vals.insert(vals.begin() + index, val);
				m_pgcache.put_page(currpos, curr.make_page());
				++m_size;
				break;
			}

			auto &refs = curr.branch().m_refs;
			auto &links = curr.branch().m_links;

			const std::size_t index = std::lower_bound(refs.cbegin(), refs.cend(), key) - refs.cbegin();
			const Position child_pos = links[index];
			assert(child_pos.is_set());
			auto child = Nod::from_page(m_pgcache.get_page(child_pos));

			if (!is_node_full(child)) {
				currpos = child_pos;
				curr = std::move(child);
				continue;
			}

			auto [midkey, sibling] = node_split(child);
			auto sibling_pos = m_pgcache.get_new_pos();

			assert(std::find(refs.cbegin(), refs.cend(), midkey) == refs.cend());
			assert(std::find(links.cbegin(), links.cend(), sibling_pos) == links.cend());

			refs.insert(refs.begin() + index, midkey);
			links.insert(links.begin() + index + 1, sibling_pos);

			m_pgcache.put_page(sibling_pos, std::move(sibling.make_page()));
			m_pgcache.put_page(child_pos, std::move(child.make_page()));
			m_pgcache.put_page(currpos, std::move(curr.make_page()));

			if (key < midkey) {
				currpos = child_pos;
				curr = std::move(child);
			} else if (key > midkey) {
				currpos = sibling_pos;
				curr = std::move(sibling);
			}
		}
	}

	[[nodiscard]] std::optional<Val> get(const Self::Key &key) const noexcept {
		return search_subtree(root(), key);
	}

	[[nodiscard]] bool contains(const Self::Key &key) const noexcept {
		return search_subtree(root(), key).has_value();
	}

private:
	void save_header() const noexcept {
		if (m_header.dirty()) {
			m_header.rootpos() = m_rootpos;
			m_header.size() = m_size;
			m_header.dirty() = false;
		}

		nop::Serializer<nop::StreamWriter<std::ofstream>> serializer{m_header_name.data(), std::ios::trunc};
		serializer.Write(m_header) || nop::Die(std::cerr);
	}

	bool load_header() {
		nop::Deserializer<nop::StreamReader<std::ifstream>> deserializer{m_header_name.data()};
		deserializer.Read(&m_header) || nop::Die(std::cerr);
		std::cout << m_header << "\n";
		return true;
	}

	void bare_init() noexcept {
		m_rootpos = m_pgcache.get_new_pos();
		Nod root_initial{typename Nod::Metadata(typename Nod::Leaf({}, {})), m_rootpos, true};
		m_pgcache.put_page(m_rootpos, std::move(root_initial.make_page()));

		// Fill in initial header
		m_header.rootpos() = m_rootpos;
		m_header.size() = m_size;
		m_header.depth() = m_depth;
	}

public:
	/*
	 *  Persistence API
	 */

	void load() {
		auto ok = load_header();
		assert(ok);

		m_rootpos = m_header.m_rootpos;
		m_size = m_header.m_size;
		m_depth = m_header.m_depth;
	}

	void save() const {
		save_header();
		m_pgcache.flush_all();
	}

public:
	explicit Btree(std::string_view pgcache_name,
	               std::string_view btree_header_name = "/tmp/eu-btree-header",
	               bool should_load = false)
	    : m_pgcache{pgcache_name, PAGE_CACHE_SIZE},
	      m_header_name{btree_header_name} {

		assert(NUM_RECORDS_BRANCH > 1);
		assert(NUM_RECORDS_LEAF > 1);

		if (should_load)
			load();
		else
			bare_init();
	}

private:
	mutable PageCache m_pgcache;

	mutable Header m_header;
	const std::string_view m_header_name;

	Position m_rootpos;

	std::size_t m_size{0};

	std::size_t m_depth{0};
};

}// namespace internal::storage::btree