text:
/*
 *     Copyright 2021-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */

#pragma once

#include <mcbp/protocol/status.h>
#include <cstdint>
#include <functional>
#include <gsl/gsl>

struct EngineIface;

/**
 * Callback for any function producing stats.
 *
 * @param key the stat's key
 * @param value the stat's value in an ascii form (e.g. text form of a number)
 * @param cookie magic callback cookie
 */
using AddStatFn = std::function<void(std::string_view key,
                                     std::string_view value,
                                     gsl::not_null<const void*> cookie)>;

/**
 * Callback for adding a response backet
 * @param key The key to put in the response
 * @param extras The data to put in the extended field in the response
 * @param body The data body
 * @param datatype This is currently not used and should be set to 0
 * @param status The status code of the return packet (see in protocol_binary
 *               for the legal values)
 * @param cas The cas to put in the return packet
 * @param cookie The cookie provided by the frontend
 * @return true if return message was successfully created, false if an
 *              error occured that prevented the message from being sent
 */
using AddResponseFn = std::function<bool(std::string_view key,
                                         std::string_view extras,
                                         std::string_view body,
                                         uint8_t datatype,
                                         cb::mcbp::Status status,
                                         uint64_t cas,
                                         const void* cookie)>;

ext: h


text:
/*This program evaluates the differential equations for a photon's geodesic in a perturbed spacetime in SPHERICAL coordinates and restricted to RADIAL MOTION with a Runge-Kutta-Fehlberg 45 method.*/

/*This program solves the particular case for flat FRW perturbed spacetime with metric: $g_{ab} = {[g]}_{ab} + h_{ab}$. Where ${[g]}_{ab}$ corresponds to the flat FRW metric and $h_{ab}$ corresponds to the perturbation in the conformal Newtonian gauge. A Plummer potential with adequate parameters have been used to simulate the perturbation.
The equations are written in the form $\frac{d(x or p)^{\alpha}}{d\lambda}=f(x^{\alpha},p^{\beta})$ and the indices $\alpha$ and $\beta$ run from 0 to 1 since the motion is only radial.
Where $p^{\alpha}={\dot{x}}^{\alpha}$.
The coordinates for the photon's geodesics are then: (ct,r) = (x0,x1).*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>          //GSL error management module
#include <gsl/gsl_spline.h>         //GSL interpolation module

/* PHYSICAL CONSTANTS AND PARAMETERS */
#define A 10.0     //Distance parameter of the perturbations
#define G 43007.01     //Gravitational constant
#define M 0.0     //Mass of the perturbation
#define C 299792.458  //Speed of light

/* PROGRAM PARAMETERS (ONLY USING NLINES FRW) */
#define NLINES 100000 //Number of lines in geodesic_solution.dat file
#define NSTEPS 75000000 //Number of steps for solving geodesics
#define NLINESFRW 10000 //Number of lines in frw.dat file
#define DLAMBDA 0.01   //Geodesics parameter step

//We define mydbl as a generic variable type which can be changed when needed to: float, double, long double
typedef long double mydbl;

/* RUNGE-KUTTA-FEHLBERG METHOD CONSTANTS */
const mydbl c2 = 0.25, c3 = 3.0/8.0, c4 = 12.0/13.0, c5 = 1.0, c6 = 0.5; //Multiplying increment in t
const mydbl a21 = 0.25, a31 = 3.0/32.0, a32 = 9.0/32.0, a41 = 1932.0/2197.0, a42 = -7200.0/2197.0, a43 = 7296.0/2197.0, a51 = 439.0/216.0, a52 = -8.0, a53 = 3680.0/513.0, a54 = -845.0/4104.0, a61 = -8.0/27.0, a62 = 2.0, a63 = -3544.0/2565.0, a64 = 1859.0/4104.0, a65 = -11.0/40.0; //This adds to the approximated solution in step i. First script means K value and second which K is multiplying
const mydbl r1 = 1.0/360.0, r3 = -128.0/4275.0, r4 = -2197.0/75240.0, r5 = 1.0/50.0, r6 = 2.0/55.0; //Coefficients to calculate R
const mydbl b1 = 25.0/216.0, b3 = 1408.0/2565.0, b4 = 2197.0/4104.0, b5 = -1.0/5.0;

/* FUNCTIONS OF THE PROGRAM */

/*Find and return the maximum value for objects in a vector*/
mydbl maxValue(mydbl arr[], size_t size)
{
  size_t i;
  mydbl maxValue = arr[0];

  for(i = 1; i < size; i++)
    {
      if(arr[i] > maxValue)
	{
	  maxValue = arr[i];
	}
    }
  return maxValue;
}

/*Interpolation of function inside *spline object evaluated at an abscisa 'x'.
Argument *spline is a pointer to a spline object which stores the type of interpolation to be made. x is the independent variable where the function is evaluated. *acc is a pointer to a lookup object for interpolations.*/
double interpolator(gsl_spline *spline, double x, gsl_interp_accel *acc)
{
  double a = gsl_spline_eval(spline, x, acc);  //Interpolates data to abcisa x using method in spline and acceleration object acc

  return a;  //Return value of interpolated function at x
}

/*Function for the gravitational potential to be used. Potential for Plummer model.*/
mydbl potential(mydbl r)
{
  return -G*M/(sqrtl(A*A + r*r));
}

/*Derivative of potential respecto to radial coordinate.*/
mydbl der_potential(mydbl r)
{
  return G*M*r/(powl(A*A+r*r, 1.5));
}

/*Function of the 0th momentum component differential equation for the geodesics.
${p0}^{dot} = f0(x^{\alpha},p^{\alpha})$.*/
mydbl geodesic_equation_0(gsl_spline *spline1, gsl_interp_accel *acc1, gsl_spline *spline2, gsl_interp_accel *acc2, mydbl p0, mydbl pr, mydbl x0, mydbl r)
{
  double t = (double)(1.0*x0/C);
  mydbl a = (mydbl) 1.0*interpolator(spline1, t, acc1);
  mydbl adot = (mydbl) 1.0*interpolator(spline2, t, acc2);
  mydbl f = -2.0*der_potential(r)*p0*pr/(C*C + 2.0*potential(r)) - (1.0 - 2.0*potential(r)/(C*C))*(a*adot)*(pr*pr)/(C + 2.0*potential(r)/C);
  return f;
}

/*Function of the 1th (radial) momentum component differential equation for the geodesics.
${p1}^{dot} = f1(x^{\alpha},p^{\alpha})$.*/
mydbl geodesic_equation_r(gsl_spline *spline1, gsl_interp_accel *acc1, gsl_spline *spline2, gsl_interp_accel *acc2, mydbl p0, mydbl pr, mydbl x0, mydbl r)
{
  double t = (double)(1.0*x0/C);
  mydbl a = (mydbl) 1.0*interpolator(spline1, t, acc1);
  mydbl adot = (mydbl) 1.0*interpolator(spline2, t, acc2);
  mydbl f = - (der_potential(r)*p0*p0)/(a*a*(C*C - 2.0*potential(r))) - (2.0*adot*p0*pr)/(C*a) + der_potential(r)*(pr*pr)/(C*C - 2.0*potential(r));
  return f;
}

/*To set the initial value of pr, it must hold $g_{\mu\nu}p^{\mu}p^{\nu} = 0$.
This factor multiplies p0 to guarantee that p1 fulfill the null geodesic condition.*/
mydbl condition_factor(mydbl r, double a)
{
  return (mydbl)(1.0/a)*sqrtl((1.0+2.0*potential(r)/(C*C))/(1.0 - 2.0*potential(r)/(C*C)));
}

/*$cp^{0}$ multiplied by this factor allows to obtain the energy for a local inertial observer in this spacetime.*/
mydbl energy_factor(mydbl r)
{
  mydbl g = sqrtl(1.0 + 2.0*potential(r)/(C*C));
  return g;
}

/*Violation of null geodesics condition $g_{\mu\nu}p^{\mu}p^{\nu} = 0$.*/
mydbl violation(mydbl r, mydbl p0, mydbl pr, double a)
{
  mydbl f = -(1.0+2.0*potential(r)/(C*C))*p0*p0 + (mydbl)(1.0*a*a)*(1.0-2.0*potential(r)/(C*C))*pr*pr;
  return f;
}

/*Function for solving the geodesic's differential equations using a Runge-Kutta-Fehlber method with adaptative stepsize.
Arguments are pointer so variables in that memory addresses are changed every time this function is called. There is no need to define auxiliary variables to store the results of computations.*/
void runge_kutta_fehlberg(gsl_spline *spline1, gsl_interp_accel *acc1, gsl_spline *spline2, gsl_interp_accel *acc2, mydbl *x0, mydbl *x1, mydbl *p0, mydbl *p1, mydbl *lambda, mydbl *dlambda, mydbl lambda_endpoint, mydbl tol, mydbl dlambda_max, mydbl dlambda_min, FILE *fp, double aem, mydbl energy1)
{
  /*Define the k quantities for all the variables*/
  mydbl k1x0, k1x1, k1p0, k1p1;
  mydbl k2x0, k2x1, k2p0, k2p1;
  mydbl k3x0, k3x1, k3p0, k3p1;
  mydbl k4x0, k4x1, k4p0, k4p1;
  mydbl k5x0, k5x1, k5p0, k5p1;
  mydbl k6x0, k6x1, k6p0, k6p1;

  /*Define a vector that carries the errors for each differential equation and a variable for storing the biggest error*/
  mydbl r[4], rmax;

  /*Define the increments for the values to be solved in the differential equations*/
  mydbl dp0, dp1, dx0, dx1;

  /*Define some relevant quantities to be printed in a file*/
  mydbl energy, difft, t, v, difference;
  double aobs, difftfrw;

  /*Relevant quantity to determine the new stepsize*/
  mydbl delta;

  while(*lambda < lambda_endpoint)
    {
      /*Check if the independent variable plus the actual stepsize is still in the desired interval*/
      if((*lambda + *dlambda) > lambda_endpoint)
	{
	  *dlambda = lambda_endpoint - *lambda;
	}

      /*Calculate the k quantities for all the variables*/
      /*This section calculates the k1 quantities*/
      k1x0 = (*dlambda)*(*p0); k1x1 = (*dlambda)*(*p1);
      k1p0 = *dlambda*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0, *p1, *x0, *x1);
      k1p1 = *dlambda*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0, *p1, *x0, *x1);

      /*This section calculates the k2 quantities*/
      k2x0 = *dlambda*(*p0 + a21*k1p0); k2x1 = *dlambda*(*p1 + a21*k1p1);
      k2p0 = *dlambda*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + a21*k1p0, *p1 + a21*k1p1, *x0 + a21*k1x0, *x1 + a21*k1x1);
      k2p1 = *dlambda*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + a21*k1p0, *p1 + a21*k1p1, *x0 + a21*k1x0, *x1 + a21*k1x1);

      /*This section calculates the k3 quantities*/
      k3x0 = *dlambda*(*p0 + a31*k1p0 + a32*k2p0); k3x1 = *dlambda*(*p1 + a31*k1p1 + a32*k2p1);
      k3p0 = *dlambda*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + a31*k1p0 + a32*k2p0, *p1 + a31*k1p1 + a32*k2p1, *x0 + a31*k1x0 + a32*k2x0, *x1 + a31*k1x1 + a32*k2x1);
      k3p1 = *dlambda*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + a31*k1p0 + a32*k2p0, *p1 + a31*k1p1 + a32*k2p1, *x0 + a31*k1x0 + a32*k2x0, *x1 + a31*k1x1 + a32*k2x1);

      /*This section calculates the k4 quantities*/
      k4x0 = *dlambda*(*p0 + a41*k1p0 + a42*k2p0 + a43*k3p0); k4x1 = *dlambda*(*p1 + a41*k1p1 + a42*k2p1 + a43*k3p1);
      k4p0 = *dlambda*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + a41*k1p0 + a42*k2p0 + a43*k3p0, *p1 + a41*k1p1 + a42*k2p1 + a43*k3p1, *x0 + a41*k1x0 + a42*k2x0 + a43*k3x0, *x1 + a41*k1x1 + a42*k2x1 + a43*k3x1);
      k4p1 = *dlambda*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + a41*k1p0 + a42*k2p0 + a43*k3p0, *p1 + a41*k1p1 + a42*k2p1 + a43*k3p1, *x0 + a41*k1x0 + a42*k2x0 + a43*k3x0, *x1 + a41*k1x1 + a42*k2x1 + a43*k3x1);

      /*This section calculates the k5 quantities*/
      k5x0 = *dlambda*(*p0 + a51*k1p0 + a52*k2p0 + a53*k3p0 + a54*k4p0); k5x1 = *dlambda*(*p1 + a51*k1p1 + a52*k2p1 + a53*k3p1 + a54*k4p1);
      k5p0 = *dlambda*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + a51*k1p0 + a52*k2p0 + a53*k3p0 + a54*k4p0, *p1 + a51*k1p1 + a52*k2p1 + a53*k3p1 + a54*k4p1, *x0 + a51*k1x0 + a52*k2x0 + a53*k3x0 + a54*k4x0, *x1 + a51*k1x1 + a52*k2x1 + a53*k3x1 + a54*k4x1);
      k5p1 = *dlambda*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + a51*k1p0 + a52*k2p0 + a53*k3p0 + a54*k4p0, *p1 + a51*k1p1 + a52*k2p1 + a53*k3p1 + a54*k4p1, *x0 + a51*k1x0 + a52*k2x0 + a53*k3x0 + a54*k4x0, *x1 + a51*k1x1 + a52*k2x1 + a53*k3x1 + a54*k4x1);

      /*This section calculates the k6 quantities*/
      k6x0 = *dlambda*(*p0 + a61*k1p0 + a62*k2p0 + a63*k3p0 + a64*k4p0 + a65*k5p0); k6x1 = *dlambda*(*p1 + a61*k1p1 + a62*k2p1 + a63*k3p1 + a64*k4p1 + a65*k5p1);
      k6p0 = *dlambda*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + a61*k1p0 + a62*k2p0 + a63*k3p0 + a64*k4p0 + a65*k5p0, *p1 + a61*k1p1 + a62*k2p1 + a63*k3p1 + a64*k4p1 + a65*k5p1, *x0 + a61*k1x0 + a62*k2x0 + a63*k3x0 + a64*k4x0 + a65*k5x0, *x1 + a61*k1x1 + a62*k2x1 + a63*k3x1 + a64*k4x1 + a65*k5x1);
      k6p1 = *dlambda*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + a61*k1p0 + a62*k2p0 + a63*k3p0 + a64*k4p0 + a65*k5p0, *p1 + a61*k1p1 + a62*k2p1 + a63*k3p1 + a64*k4p1 + a65*k5p1, *x0 + a61*k1x0 + a62*k2x0 + a63*k3x0 + a64*k4x0 + a65*k5x0, *x1 + a61*k1x1 + a62*k2x1 + a63*k3x1 + a64*k4x1 + a65*k5x1);

      /*Determination of the biggest error between the errors for each differential equation*/
      r[0] = fabsl(r1*k1p0 + r3*k3p0 + r4*k4p0 + r5*k5p0 + r6*k6p0)/(*dlambda);
      r[1] = fabsl(r1*k1p1 + r3*k3p1 + r4*k4p1 + r5*k5p1 + r6*k6p1)/(*dlambda);
      r[2] = fabsl(r1*k1x0 + r3*k3x0 + r4*k4x0 + r5*k5x0 + r6*k6x0)/(*dlambda);
      r[3] = fabsl(r1*k1x1 + r3*k3x1 + r4*k4x1 + r5*k5x1 + r6*k6x1)/(*dlambda);
      rmax = maxValue(r, 4);

      /*Check whether the worst approximation is less than the introduced tolerance*/
      if(rmax <= tol)
	{
	  *lambda = *lambda + *dlambda;

	  dp0 = b1*k1p0 + b3*k3p0 + b4*k4p0 + b5*k5p0;
	  dp1 = b1*k1p1 + b3*k3p1 + b4*k4p1 + b5*k5p1;
	  dx0 = b1*k1x0 + b3*k3x0 + b4*k4x0 + b5*k5x0;
	  dx1 = b1*k1x1 + b3*k3x1 + b4*k4x1 + b5*k5x1;
	  *p0 = *p0 + dp0; *p1 = *p1 + dp1; *x0 = *x0 + dx0; *x1 = *x1 + dx1;

	  /*Relevant quantities to save in the .dat file. These only are calculated if the approximation is accepted*/
	  energy = C*energy_factor(*x1)*(*p0);
  	  difft = (energy - energy1)/energy1;
  	  t = *x0/C;
  	  aobs = interpolator(spline1, (double)(1.0*t), acc1);
  	  v = violation(*x1, *p0, *p1, aobs);
  	  difftfrw = (aem/aobs) - 1.0;
  	  difference = difft - (mydbl)(1.0*difftfrw);
  	  fprintf(fp, "%16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8e %16.8Le %16.8Le\n", *lambda, *dlambda, *x0, *x1, *p0, *p1, energy, v, difft, difftfrw, difference, rmax);
	}

      /*Calculate a new stepsize*/
      delta = 0.84*powl(tol/rmax, 0.25);

      if(delta <= 0.1)
	{
	  *dlambda = *dlambda*0.1;
	}
      else if(delta >= 4.0)
	{
	  *dlambda = *dlambda*4.0;
	}
      else
	{
	  *dlambda = *dlambda*delta;
	}

      /*Check whether the stepsize is too big or too small*/
      if(*dlambda > dlambda_max)
	{
	  *dlambda = dlambda_max;
	}
      else if(*dlambda < dlambda_min)
	{
	  printf("Minimum value for stepsize dlambda exceeded.\n");
	  break;
	}
    }
}

int main(void)
{
  /***READ SCALE FACTOR DATA AND PREPARE OBJECTS FOR INTERPOLATION ***/
  int i;            //For array manipulation

  /*Pointer to scale_factor.data file*/
  FILE *frw;
  frw = fopen("../scale_factor.dat","r");

  /*Variables and arrays to read the data*/
  double cosmictime[NLINESFRW], conftime, scale_factor[NLINESFRW], der_scale_factor[NLINESFRW];

  /*Reading the data*/
  for(i=0; i<NLINESFRW; i++)
    {
      fscanf(frw,"%lf %lf %lf %lf", &cosmictime[i], &conftime, &scale_factor[i], &der_scale_factor[i]);
    }

  /*Free space in memory*/
  fclose(frw);

  /*** Initializes objects for interpolation. 1 is for interpolation of scale factor, 2 is for interpolation of derivative of scale factor ***/

  /*Allocate space in memory*/
  gsl_interp_accel *acc1 = gsl_interp_accel_alloc();  //Acceleration type object (for index lookup)
  gsl_interp_accel *acc2 = gsl_interp_accel_alloc();
  gsl_spline *spline1 = gsl_spline_alloc(gsl_interp_cspline, NLINESFRW);  //Spline type object (define interpolation type and space in memory, works for both)
  gsl_spline *spline2 = gsl_spline_alloc(gsl_interp_cspline, NLINESFRW);

  /*Initializes objects for interpolation*/
  gsl_spline_init(spline1, cosmictime, scale_factor, NLINESFRW);  //Initializes spline object for data cosmictime, scale_factor of size NLINES
  gsl_spline_init(spline2, cosmictime, der_scale_factor, NLINESFRW);  //Initializes spline object for data cosmictime, der_scale_factor of size NLINES

  /************************************************************************************/

  /***SOLVES GEODESIC EQUATIONS FOR PERTURBED FRW UNIVERSE WITH STATIC PLUMMER POTENTIAL ***/

  /*Parameters for RKF method*/
  mydbl dlambda_max = 1.0e-01, dlambda_min = 1.0e-05, tol = 1.0e-10, lambdaEndPoint = 1.0e+05;

  /*Initial conditions*/
  mydbl ti = 7.0 ,x0, r = -500.0, p0 = 1.0e-3, pr, lambda = 0.0, energy1, energy, v, difft, difference, dlambda = dlambda_max;
  double difftfrw, aem, aobs;
  x0 = C*ti;
  aem = interpolator(spline1, (double)(1.0*ti), acc1);
  pr = condition_factor(r, aem)*p0;
  energy1 = C*energy_factor(r)*p0;
  v = violation(r, p0, pr, aem);
  difft = (energy1 - energy1)/energy1;
  difftfrw = (aem/aem) - 1.0;
  difference = difft - (mydbl)(1.0*difftfrw);

  /*Pointer to file where solution of differential equation will be saved.*/
  FILE *geodesic;
  geodesic = fopen("geodesic_solution.dat","w");

  /*Write line of initial values in file*/
  fprintf(geodesic, "%16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8e %16.8Le %16.8Le\n", lambda, dlambda, x0, r, p0, pr, energy1, v, difft, difftfrw, difference, (mydbl)0.0);

  /*Solution of the differential equation*/
  runge_kutta_fehlberg(spline1, acc1, spline2, acc2, &x0, &r, &p0, &pr, &lambda, &dlambda, lambdaEndPoint, tol, dlambda_max, dlambda_min, geodesic, aem, energy1);

  /************************************************************************************/

  /***RELEASING ALL SPACE USED IN MEMORY***/
  fclose(geodesic); //Close file storing the results
  gsl_spline_free(spline1);  //Free memory of spline object
  gsl_spline_free(spline2);
  gsl_interp_accel_free(acc1);  //Free memory of accel object
  gsl_interp_accel_free(acc2);
}

ext: c


text:
/**
* Copyright 2016 BitTorrent Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*    http://www.apache.org/licenses/LICENSE-2.0
*
* 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 <scraps/config.h>

#include <scraps/base64.h>
#include <scraps/Byte.h>
#include <scraps/hex.h>
#include <scraps/Temp.h>
#include <scraps/random.h>
#include <scraps/hash.h>
#include <stdts/optional.h>

#include <gsl.h>

#include <cerrno>
#include <vector>
#include <iterator>
#include <random>
#include <chrono>

namespace scraps {

constexpr double kPi = 3.1415926535897932385;

/**
* Returns a string escaped according to RFC 4627.
*
* @param the string to be escaped
* @return the string escaped according to RFC 4627
*/
std::string JSONEscape(const char* str);

/**
* Returns a string in which all non-alphanumeric characters except dashes, underscores,
* spaces, and periods are replaced with a percent sign followed by their hexadecimal
* value. Spaces are replaced with plus signs.
*
* @return the url-encoded string
*/
std::string URLEncode(const char* str);

inline std::string URLEncode(const std::string& str) { return URLEncode(str.c_str()); }

/**
* Returns a string in which the effects of URLEncode have been reversed.
*
* @return the url-decoded string
*/
std::string URLDecode(const char* str);

inline std::string URLDecode(const std::string& str) { return URLDecode(str.c_str()); }

/**
 * Clamp a value between min and max, inclusive.
 */
template <typename T, typename MinT, typename MaxT>
constexpr auto Clamp(const T& value, const MinT& min, const MaxT& max) {
    return std::max<std::common_type_t<T, MinT, MaxT>>(min, std::min<std::common_type_t<T, MinT, MaxT>>(value, max));
}

/**
* Parse out an address and port from the given host string. If a port is not found, the
* defaultPort is returned.
*/
std::tuple<std::string, uint16_t> ParseAddressAndPort(const std::string& host, uint16_t defaultPort);

/**
* @return the amount of physical memory that the system has, or 0 on error
*/
size_t PhysicalMemory();

/**
* Iterates over a container, allowing for safe modification of the container at any point.
*/
template <typename T, typename F>
void NonatomicIteration(T&& iterable, F&& function) {
    auto copy = iterable;
    for (auto& element : copy) {
        auto it = std::find(iterable.begin(), iterable.end(), element);
        if (it == iterable.end()) {
            continue;
        }
        function(element);
    }
}

template <typename T, std::ptrdiff_t N>
gsl::basic_string_span<T> TrimLeft(gsl::basic_string_span<T, N> str) {
    decltype(str.size()) i = 0;
    while(i < str.size() && isspace(str[i])) { ++i; };
    return {str.data() + i, static_cast<std::ptrdiff_t>(str.size() - i)};
}

template <typename T, std::ptrdiff_t N>
gsl::basic_string_span<T> TrimRight(gsl::basic_string_span<T, N> str) {
    decltype(str.size()) i = str.size();
    while(i > 0 && isspace(str[i - 1])) { --i; };
    return {str.data(), static_cast<std::ptrdiff_t>(i)};
}

template <typename T, std::ptrdiff_t N>
auto Trim(gsl::basic_string_span<T, N> str) {
    return TrimLeft(TrimRight(str));
}

stdts::optional<std::vector<Byte>> BytesFromFile(const std::string& path);

/**
* Sets the given file descriptor to blocking or non-blocking.
*
* @param fd the file descriptor
* @param blocking true if the file descriptor should be set to blocking. false if it should be set to non-blocking
* @return true on success
*/
bool SetBlocking(int fd, bool blocking = true);

/**
* Returns a demangled symbol name. If demangling is not supported, returns original mangled name.
*/
std::string Demangle(const char* name);

/**
* Returns whether two strings are equal disregarding case
*/
inline bool CaseInsensitiveEquals(stdts::string_view l, stdts::string_view r) {
    return l.size() == r.size() && std::equal(l.begin(), l.end(), r.begin(), [](char a, char b){ return tolower(a) == tolower(b); });
}

/**
* Split a sequence on a delimiter. The SequenceType template argument should be
* specified to the type of container OutputIt should accept. By default, it is
* set OutputIt::container_type::value_type.
*
* @param begin beginning of input sequence
* @param end end of input sequence
* @param out output iterator to a container of Sequence types
* @param delimiter a single element element in the sequence
*/
template <typename InputIt, typename OutputIt, typename Delimiter, typename SequenceType = typename OutputIt::container_type::value_type>
void Split(InputIt&& begin, InputIt&& end, OutputIt&& out, Delimiter&& delimiter) {
    SequenceType current;
    auto inserter = std::back_inserter(current);
    for (auto it = begin; it != end; ++it) {
        if (*it == delimiter) {
            *(out++) = std::move(current);
            current = {};
            inserter = std::back_inserter(current);
        } else {
            inserter = *it;
        }
    }
    *out = std::move(current);
}

} // namespace scraps

max_stars_repo_name: carlbrown/scraps

max_issues_repo_name: carlbrown/scraps

max_forks_repo_name: carlbrown/scraps


text:
/* integration/gsl_integration.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_INTEGRATION_H__
#define __GSL_INTEGRATION_H__
#include <stdlib.h>
#include <gsl/gsl_math.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* Workspace for adaptive integrators */

typedef struct
  {
    size_t limit;
    size_t size;
    size_t nrmax;
    size_t i;
    size_t maximum_level;
    double *alist;
    double *blist;
    double *rlist;
    double *elist;
    size_t *order;
    size_t *level;
  }
gsl_integration_workspace;

gsl_integration_workspace *
  gsl_integration_workspace_alloc (const size_t n);

void
  gsl_integration_workspace_free (gsl_integration_workspace * w);

/* Workspace for QAWS integrator */

typedef struct
{
  double alpha;
  double beta;
  int mu;
  int nu;
  double ri[25];
  double rj[25];
  double rg[25];
  double rh[25];
}
gsl_integration_qaws_table;

gsl_integration_qaws_table *
gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu);

int
gsl_integration_qaws_table_set (gsl_integration_qaws_table * t,
                                double alpha, double beta, int mu, int nu);

void
gsl_integration_qaws_table_free (gsl_integration_qaws_table * t);

/* Workspace for QAWO integrator */

enum gsl_integration_qawo_enum { GSL_INTEG_COSINE, GSL_INTEG_SINE };

typedef struct
{
  size_t n;
  double omega;
  double L;
  double par;
  enum gsl_integration_qawo_enum sine;
  double *chebmo;
}
gsl_integration_qawo_table;

gsl_integration_qawo_table *
gsl_integration_qawo_table_alloc (double omega, double L,
                                  enum gsl_integration_qawo_enum sine,
                                  size_t n);

int
gsl_integration_qawo_table_set (gsl_integration_qawo_table * t,
                                double omega, double L,
                                enum gsl_integration_qawo_enum sine);

int
gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t,
                                       double L);

void
gsl_integration_qawo_table_free (gsl_integration_qawo_table * t);

/* Definition of an integration rule */

typedef void gsl_integration_rule (const gsl_function * f,
                                   double a, double b,
                                   double *result, double *abserr,
                                   double *defabs, double *resabs);

void gsl_integration_qk15 (const gsl_function * f, double a, double b,
                           double *result, double *abserr,
                           double *resabs, double *resasc);

void gsl_integration_qk21 (const gsl_function * f, double a, double b,
                           double *result, double *abserr,
                           double *resabs, double *resasc);

void gsl_integration_qk31 (const gsl_function * f, double a, double b,
                           double *result, double *abserr,
                           double *resabs, double *resasc);

void gsl_integration_qk41 (const gsl_function * f, double a, double b,
                           double *result, double *abserr,
                           double *resabs, double *resasc);

void gsl_integration_qk51 (const gsl_function * f, double a, double b,
                           double *result, double *abserr,
                           double *resabs, double *resasc);

void gsl_integration_qk61 (const gsl_function * f, double a, double b,
                           double *result, double *abserr,
                           double *resabs, double *resasc);

void gsl_integration_qcheb (gsl_function * f, double a, double b,
                            double *cheb12, double *cheb24);

/* The low-level integration rules in QUADPACK are identified by small
   integers (1-6). We'll use symbolic constants to refer to them.  */

enum
  {
    GSL_INTEG_GAUSS15 = 1,      /* 15 point Gauss-Kronrod rule */
    GSL_INTEG_GAUSS21 = 2,      /* 21 point Gauss-Kronrod rule */
    GSL_INTEG_GAUSS31 = 3,      /* 31 point Gauss-Kronrod rule */
    GSL_INTEG_GAUSS41 = 4,      /* 41 point Gauss-Kronrod rule */
    GSL_INTEG_GAUSS51 = 5,      /* 51 point Gauss-Kronrod rule */
    GSL_INTEG_GAUSS61 = 6       /* 61 point Gauss-Kronrod rule */
  };

void
gsl_integration_qk (const int n, const double xgk[],
                    const double wg[], const double wgk[],
                    double fv1[], double fv2[],
                    const gsl_function *f, double a, double b,
                    double * result, double * abserr,
                    double * resabs, double * resasc);

int gsl_integration_qng (const gsl_function * f,
                         double a, double b,
                         double epsabs, double epsrel,
                         double *result, double *abserr,
                         size_t * neval);

int gsl_integration_qag (const gsl_function * f,
                         double a, double b,
                         double epsabs, double epsrel, size_t limit,
                         int key,
                         gsl_integration_workspace * workspace,
                         double *result, double *abserr);

int gsl_integration_qagi (gsl_function * f,
                          double epsabs, double epsrel, size_t limit,
                          gsl_integration_workspace * workspace,
                          double *result, double *abserr);

int gsl_integration_qagiu (gsl_function * f,
                           double a,
                           double epsabs, double epsrel, size_t limit,
                           gsl_integration_workspace * workspace,
                           double *result, double *abserr);

int gsl_integration_qagil (gsl_function * f,
                           double b,
                           double epsabs, double epsrel, size_t limit,
                           gsl_integration_workspace * workspace,
                           double *result, double *abserr);

int gsl_integration_qags (const gsl_function * f,
                          double a, double b,
                          double epsabs, double epsrel, size_t limit,
                          gsl_integration_workspace * workspace,
                          double *result, double *abserr);

int gsl_integration_qagp (const gsl_function * f,
                          double *pts, size_t npts,
                          double epsabs, double epsrel, size_t limit,
                          gsl_integration_workspace * workspace,
                          double *result, double *abserr);

int gsl_integration_qawc (gsl_function *f,
                          const double a, const double b, const double c,
                          const double epsabs, const double epsrel, const size_t limit,
                          gsl_integration_workspace * workspace,
                          double * result, double * abserr);

int gsl_integration_qaws (gsl_function * f,
                          const double a, const double b,
                          gsl_integration_qaws_table * t,
                          const double epsabs, const double epsrel,
                          const size_t limit,
                          gsl_integration_workspace * workspace,
                          double *result, double *abserr);

int gsl_integration_qawo (gsl_function * f,
                          const double a,
                          const double epsabs, const double epsrel,
                          const size_t limit,
                          gsl_integration_workspace * workspace,
                          gsl_integration_qawo_table * wf,
                          double *result, double *abserr);

int gsl_integration_qawf (gsl_function * f,
                          const double a,
                          const double epsabs,
                          const size_t limit,
                          gsl_integration_workspace * workspace,
                          gsl_integration_workspace * cycle_workspace,
                          gsl_integration_qawo_table * wf,
                          double *result, double *abserr);

/* Workspace for fixed-order Gauss-Legendre integration */

typedef struct
  {
    size_t n;         /* number of points */
    double *x;        /* Gauss abscissae/points */
    double *w;        /* Gauss weights for each abscissae */
    int precomputed;  /* high precision abscissae/weights precomputed? */
  }
gsl_integration_glfixed_table;

gsl_integration_glfixed_table * gsl_integration_glfixed_table_alloc (size_t n);

void gsl_integration_glfixed_table_free (gsl_integration_glfixed_table * t);

/* Routine for fixed-order Gauss-Legendre integration */

double gsl_integration_glfixed (const gsl_function *f,
                                double a,
                                double b,
                                const gsl_integration_glfixed_table * t);

/* Routine to retrieve the i-th Gauss-Legendre point and weight from t */

int gsl_integration_glfixed_point (double a,
                                   double b,
                                   size_t i,
                                   double *xi,
                                   double *wi,
                                   const gsl_integration_glfixed_table * t);

/* Cquad integration - Pedro Gonnet */

/* Data of a single interval */
typedef struct
{
  double a, b;
  double c[64];
  double fx[33];
  double igral, err;
  int depth, rdepth, ndiv;
} gsl_integration_cquad_ival;

/* The workspace is just a collection of intervals */
typedef struct
{
  size_t size;
  gsl_integration_cquad_ival *ivals;
  size_t *heap;
} gsl_integration_cquad_workspace;

gsl_integration_cquad_workspace *
gsl_integration_cquad_workspace_alloc (const size_t n);

void
gsl_integration_cquad_workspace_free (gsl_integration_cquad_workspace * w);

int
gsl_integration_cquad (const gsl_function * f, double a, double b,
		                   double epsabs, double epsrel,
		                   gsl_integration_cquad_workspace * ws,
		                   double *result, double *abserr, size_t * nevals);

/* Romberg integration workspace and routines */

typedef struct
{
  size_t n;       /* maximum number of steps */
  double *work1;  /* workspace for a row of R matrix, size n */
  double *work2;  /* workspace for a row of R matrix, size n */
} gsl_integration_romberg_workspace;

gsl_integration_romberg_workspace *gsl_integration_romberg_alloc(const size_t n);
void gsl_integration_romberg_free(gsl_integration_romberg_workspace * w);
int gsl_integration_romberg(const gsl_function * f, const double a, const double b,
                            const double epsabs, const double epsrel, double * result,
                            size_t * neval, gsl_integration_romberg_workspace * w);

/* IQPACK related structures and routines */

typedef struct
{
  double alpha;
  double beta;
  double a;
  double b;
  double zemu;
  double shft;
  double slp;
  double al;
  double be;
} gsl_integration_fixed_params;

typedef struct
{
  int (*check)(const size_t n, const gsl_integration_fixed_params * params);
  int (*init)(const size_t n, double * diag, double * subdiag, gsl_integration_fixed_params * params);
} gsl_integration_fixed_type;

typedef struct
{
  size_t n;        /* number of nodes/weights */
  double *weights; /* quadrature weights */
  double *x;       /* quadrature nodes */
  double *diag;    /* diagonal of Jacobi matrix */
  double *subdiag; /* subdiagonal of Jacobi matrix */
  const gsl_integration_fixed_type * type;
} gsl_integration_fixed_workspace;

/* IQPACK integral types */
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_legendre;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_chebyshev;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_gegenbauer;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_jacobi;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_laguerre;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_hermite;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_exponential;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_rational;
GSL_VAR const gsl_integration_fixed_type * gsl_integration_fixed_chebyshev2;

gsl_integration_fixed_workspace *
gsl_integration_fixed_alloc(const gsl_integration_fixed_type * type, const size_t n,
                            const double a, const double b, const double alpha, const double beta);

void gsl_integration_fixed_free(gsl_integration_fixed_workspace * w);

size_t gsl_integration_fixed_n(const gsl_integration_fixed_workspace * w);

double *gsl_integration_fixed_nodes(const gsl_integration_fixed_workspace * w);

double *gsl_integration_fixed_weights(const gsl_integration_fixed_workspace * w);

int gsl_integration_fixed(const gsl_function * func, double * result,
                          const gsl_integration_fixed_workspace * w);

__END_DECLS

#endif /* __GSL_INTEGRATION_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
/* specfunc/test_sincos_pi.c
 *
 * Copyright (C) 2017 Konrad Griessinger
 *
 * 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: Konrad Griessinger */

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"

/* Any double precision number bigger than this is automatically an even integer. */
#define BIGDBL (2.0 / GSL_DBL_EPSILON)

int
test_sincos_pi(void)
{
  gsl_sf_result r;
  int s = 0;
  int k = 0, kmax = 12;
  double x = 0.0, ix = 0.0, fx = 0.0, exact = 0.0;

  /* sin_pi tests */

  fx = 0.5;
  exact = 1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
  fx = -0.5;
  exact = -1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 1.5;
  exact = -1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
  fx = -1.5;
  exact = 1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 2.5;
  exact = 1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
  fx = -2.5;
  exact = -1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 3.5;
  exact = -1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
  fx = -3.5;
  exact = 1.0;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 0.375;
  exact = 0.923879532511286756128183189397;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = -0.375;
  exact = -0.923879532511286756128183189397;

  TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 0.0;
  exact = 0.0;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.5;
  exact = 1.0;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.03125;
  exact = 0.0980171403295606019941955638886;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.0625;
  exact = 0.195090322016128267848284868477;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.75;
  exact = 0.707106781186547524400844362105;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.0078125;
  exact = 0.0245412285229122880317345294593;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  /* sin_pi tests for very large arguments */

  fx = 0.0625;
  exact = 0.195090322016128267848284868477;
  ix = LONG_MAX + 1.0;
  ix += fabs(fmod(ix,2.0)); /* make sure of even number */

  for (k=0; k<kmax; k++) {
    x = ix + fx;
    x -= ix; /* careful with compiler optimization */
    if ( ( x != fx ) || ( fabs(ix+fx) >= BIGDBL ) ) break;
    printf("ix+fx= %.18e\n", ix+fx);
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix += 101.0;
    exact = -exact;
  }

  fx = -0.0625;
  exact = -0.195090322016128267848284868477;
  ix = LONG_MIN - 1.0;
  ix -= fabs(fmod(ix,2.0)); /* make sure of even number */

  for (k=0; k<kmax; k++) {
    x = ix + fx;
    x -= ix; /* careful with compiler optimization */
    if ( ( x != fx ) || ( fabs(ix+fx) >= BIGDBL ) ) break;
    printf("ix+fx= %.18e\n", ix+fx);
    TEST_SF(s, gsl_sf_sin_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix -= 101.0;
    exact = -exact;
  }

  /* cos_pi tests */

  ix = 0.0;
  fx = 0.0;
  exact = 1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 1.0;
  exact = -1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = -1.0;
  exact = -1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 2.0;
  exact = 1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = -2.0;
  exact = 1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 3.0;
  exact = -1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = -3.0;
  exact = -1.0;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 0.375;
  exact = 0.382683432365089771728459984030;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = -0.375;
  exact = 0.382683432365089771728459984030;

  TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);

  fx = 0.0;
  exact = 1.0;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.5;
  exact = 0.0;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.0625;
  exact = 0.980785280403230449126182236134;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.4375;
  exact = 0.195090322016128267848284868477;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  fx = 0.4921875;
  exact = 0.0245412285229122880317345294593;

  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(3.0,k+1);
    if (k==0) exact = -exact;
  }

  exact = fabs(exact);
  ix = 0.0;
  for (k=0; k<kmax; k++) {
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix = pow(10.0,k+1);
  }

  /* cos_pi tests for very large arguments */

  fx = 0.0625;
  exact = 0.980785280403230449126182236134;
  ix = LONG_MAX + 1.0;
  ix += fabs(fmod(ix,2.0)); /* make sure of even number */

  for (k=0; k<kmax; k++) {
    x = ix + fx;
    x -= ix; /* careful with compiler optimization */
    if ( ( x != fx ) || ( fabs(ix+fx) >= BIGDBL ) ) break;
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix += 101.0;
    exact = -exact;
  }

  fx = -0.0625;
  exact = 0.980785280403230449126182236134;
  ix = LONG_MIN - 1.0;
  ix -= fabs(fmod(ix,2.0)); /* make sure of even number */

  for (k=0; k<kmax; k++) {
    x = ix + fx;
    x -= ix; /* careful with compiler optimization */
    if ( ( x != fx ) || ( fabs(ix+fx) >= BIGDBL ) ) break;
    TEST_SF(s, gsl_sf_cos_pi_e, (ix+fx, &r), exact, TEST_TOL0, GSL_SUCCESS);
    ix -= 101.0;
    exact = -exact;
  }

  return s;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once

#include <cstring>
#include <string>
#include <sstream>
#include <iomanip>
#include <array>
#include <string_view>
#include <iterator>
#include <random>
#include <memory>
#include <functional>
#include <type_traits>
#include <assert.h>
#include <gsl/span>

#ifdef _WIN32
#include <objbase.h>
#elif defined(__linux__) || defined(__unix__)
#include <uuid/uuid.h>
#elif defined(__APPLE__)
#include <CoreFoundation/CFUUID.h>
#endif

namespace uuids
{
   namespace detail
   {
      template <typename TChar>
      constexpr inline unsigned char hex2char(TChar const ch)
      {
         if (ch >= static_cast<TChar>('0') && ch <= static_cast<TChar>('9'))
            return static_cast<unsigned char>(ch - static_cast<TChar>('0'));
         if (ch >= static_cast<TChar>('a') && ch <= static_cast<TChar>('f'))
            return static_cast<unsigned char>(10 + ch - static_cast<TChar>('a'));
         if (ch >= static_cast<TChar>('A') && ch <= static_cast<TChar>('F'))
            return static_cast<unsigned char>(10 + ch - static_cast<TChar>('A'));
         return 0;
      }

      template <typename TChar>
      constexpr inline bool is_hex(TChar const ch)
      {
         return
            (ch >= static_cast<TChar>('0') && ch <= static_cast<TChar>('9')) ||
            (ch >= static_cast<TChar>('a') && ch <= static_cast<TChar>('f')) ||
            (ch >= static_cast<TChar>('A') && ch <= static_cast<TChar>('F'));
      }

      template <typename TChar>
      constexpr inline unsigned char hexpair2char(TChar const a, TChar const b)
      {
         return static_cast<unsigned char>((hex2char(a) << 4) | hex2char(b));
      }

      class sha1
      {
      public:
         using digest32_t = uint32_t[5];
         using digest8_t = uint8_t[20];

         static constexpr unsigned int block_bytes = 64;

         inline static uint32_t left_rotate(uint32_t value, size_t const count)
         {
            return (value << count) ^ (value >> (32 - count));
         }

         sha1() { reset(); }

         void reset()
         {
            m_digest[0] = 0x67452301;
            m_digest[1] = 0xEFCDAB89;
            m_digest[2] = 0x98BADCFE;
            m_digest[3] = 0x10325476;
            m_digest[4] = 0xC3D2E1F0;
            m_blockByteIndex = 0;
            m_byteCount = 0;
         }

         void process_byte(uint8_t octet)
         {
            this->m_block[this->m_blockByteIndex++] = octet;
            ++this->m_byteCount;
            if (m_blockByteIndex == block_bytes)
            {
               this->m_blockByteIndex = 0;
               process_block();
            }
         }

         void process_block(void const * const start, void const * const end)
         {
            const uint8_t* begin = static_cast<const uint8_t*>(start);
            const uint8_t* finish = static_cast<const uint8_t*>(end);
            while (begin != finish)
            {
               process_byte(*begin);
               begin++;
            }
         }

         void process_bytes(void const * const data, size_t const len)
         {
            const uint8_t* block = static_cast<const uint8_t*>(data);
            process_block(block, block + len);
         }

         uint32_t const * get_digest(digest32_t digest)
         {
            size_t const bitCount = this->m_byteCount * 8;
            process_byte(0x80);
            if (this->m_blockByteIndex > 56) {
               while (m_blockByteIndex != 0) {
                  process_byte(0);
               }
               while (m_blockByteIndex < 56) {
                  process_byte(0);
               }
            }
            else {
               while (m_blockByteIndex < 56) {
                  process_byte(0);
               }
            }
            process_byte(0);
            process_byte(0);
            process_byte(0);
            process_byte(0);
            process_byte(static_cast<unsigned char>((bitCount >> 24) & 0xFF));
            process_byte(static_cast<unsigned char>((bitCount >> 16) & 0xFF));
            process_byte(static_cast<unsigned char>((bitCount >> 8) & 0xFF));
            process_byte(static_cast<unsigned char>((bitCount) & 0xFF));

            memcpy(digest, m_digest, 5 * sizeof(uint32_t));
            return digest;
         }

         uint8_t const * get_digest_bytes(digest8_t digest)
         {
            digest32_t d32;
            get_digest(d32);
            size_t di = 0;
            digest[di++] = (uint8_t)((d32[0] >> 24) & 0xFF);
            digest[di++] = (uint8_t)((d32[0] >> 16) & 0xFF);
            digest[di++] = (uint8_t)((d32[0] >> 8) & 0xFF);
            digest[di++] = (uint8_t)((d32[0]) & 0xFF);

            digest[di++] = (uint8_t)((d32[1] >> 24) & 0xFF);
            digest[di++] = (uint8_t)((d32[1] >> 16) & 0xFF);
            digest[di++] = (uint8_t)((d32[1] >> 8) & 0xFF);
            digest[di++] = (uint8_t)((d32[1]) & 0xFF);

            digest[di++] = (uint8_t)((d32[2] >> 24) & 0xFF);
            digest[di++] = (uint8_t)((d32[2] >> 16) & 0xFF);
            digest[di++] = (uint8_t)((d32[2] >> 8) & 0xFF);
            digest[di++] = (uint8_t)((d32[2]) & 0xFF);

            digest[di++] = (uint8_t)((d32[3] >> 24) & 0xFF);
            digest[di++] = (uint8_t)((d32[3] >> 16) & 0xFF);
            digest[di++] = (uint8_t)((d32[3] >> 8) & 0xFF);
            digest[di++] = ((d32[3]) & 0xFF);

            digest[di++] = (uint8_t)((d32[4] >> 24) & 0xFF);
            digest[di++] = (uint8_t)((d32[4] >> 16) & 0xFF);
            digest[di++] = (uint8_t)((d32[4] >> 8) & 0xFF);
            digest[di++] = (uint8_t)((d32[4]) & 0xFF);

            return digest;
         }

      private:
         void process_block()
         {
            uint32_t w[80];
            for (size_t i = 0; i < 16; i++) {
               w[i] = (m_block[i * 4 + 0] << 24);
               w[i] |= (m_block[i * 4 + 1] << 16);
               w[i] |= (m_block[i * 4 + 2] << 8);
               w[i] |= (m_block[i * 4 + 3]);
            }
            for (size_t i = 16; i < 80; i++) {
               w[i] = left_rotate((w[i - 3] ^ w[i - 8] ^ w[i - 14] ^ w[i - 16]), 1);
            }

            uint32_t a = m_digest[0];
            uint32_t b = m_digest[1];
            uint32_t c = m_digest[2];
            uint32_t d = m_digest[3];
            uint32_t e = m_digest[4];

            for (std::size_t i = 0; i < 80; ++i)
            {
               uint32_t f = 0;
               uint32_t k = 0;

               if (i < 20) {
                  f = (b & c) | (~b & d);
                  k = 0x5A827999;
               }
               else if (i < 40) {
                  f = b ^ c ^ d;
                  k = 0x6ED9EBA1;
               }
               else if (i < 60) {
                  f = (b & c) | (b & d) | (c & d);
                  k = 0x8F1BBCDC;
               }
               else {
                  f = b ^ c ^ d;
                  k = 0xCA62C1D6;
               }
               uint32_t temp = left_rotate(a, 5) + f + e + k + w[i];
               e = d;
               d = c;
               c = left_rotate(b, 30);
               b = a;
               a = temp;
            }

            m_digest[0] += a;
            m_digest[1] += b;
            m_digest[2] += c;
            m_digest[3] += d;
            m_digest[4] += e;
         }

      private:
         digest32_t m_digest;
         uint8_t m_block[64];
         size_t m_blockByteIndex;
         size_t m_byteCount;
      };
   }

   // UUID format https://tools.ietf.org/html/rfc4122
   // Field	                     NDR Data Type	   Octet #	Note
   // --------------------------------------------------------------------------------------------------------------------------
   // time_low	                  unsigned long	   0 - 3	   The low field of the timestamp.
   // time_mid	                  unsigned short	   4 - 5	   The middle field of the timestamp.
   // time_hi_and_version	      unsigned short	   6 - 7	   The high field of the timestamp multiplexed with the version number.
   // clock_seq_hi_and_reserved	unsigned small	   8	      The high field of the clock sequence multiplexed with the variant.
   // clock_seq_low	            unsigned small	   9	      The low field of the clock sequence.
   // node	                     character	      10 - 15	The spatially unique node identifier.
   // --------------------------------------------------------------------------------------------------------------------------
   // 0                   1                   2                   3
   //  0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
   // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
   // |                          time_low                             |
   // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
   // |       time_mid                |         time_hi_and_version   |
   // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
   // |clk_seq_hi_res |  clk_seq_low  |         node (0-1)            |
   // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
   // |                         node (2-5)                            |
   // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

   // indicated by a bit pattern in octet 8, marked with N in xxxxxxxx-xxxx-xxxx-Nxxx-xxxxxxxxxxxx
   enum class uuid_variant
   {
      // NCS backward compatibility (with the obsolete Apollo Network Computing System 1.5 UUID format)
      // N bit pattern: 0xxx
      // > the first 6 octets of the UUID are a 48-bit timestamp (the number of 4 microsecond units of time since 1 Jan 1980 UTC);
      // > the next 2 octets are reserved;
      // > the next octet is the "address family";
      // > the final 7 octets are a 56-bit host ID in the form specified by the address family
      ncs,

      // RFC 4122/DCE 1.1
      // N bit pattern: 10xx
      // > big-endian byte order
      rfc,

      // Microsoft Corporation backward compatibility
      // N bit pattern: 110x
      // > little endian byte order
      // > formely used in the Component Object Model (COM) library
      microsoft,

      // reserved for possible future definition
      // N bit pattern: 111x
      reserved
   };

   struct uuid_error : public std::runtime_error
   {
      explicit uuid_error(std::string_view message)
         : std::runtime_error(message.data())
      {
      }

      explicit uuid_error(char const * message)
         : std::runtime_error(message)
      {
      }
   };

   // indicated by a bit pattern in octet 6, marked with M in xxxxxxxx-xxxx-Mxxx-xxxx-xxxxxxxxxxxx
   enum class uuid_version
   {
      none = 0, // only possible for nil or invalid uuids
      time_based = 1,  // The time-based version specified in RFC 4122
      dce_security = 2,  // DCE Security version, with embedded POSIX UIDs.
      name_based_md5 = 3,  // The name-based version specified in RFS 4122 with MD5 hashing
      random_number_based = 4,  // The randomly or pseudo-randomly generated version specified in RFS 4122
      name_based_sha1 = 5   // The name-based version specified in RFS 4122 with SHA1 hashing
   };

   struct uuid
   {
      struct uuid_const_iterator
      {
         using self_type         = uuid_const_iterator;
         using value_type        = uint8_t;
         using reference         = uint8_t const &;
         using pointer           = uint8_t const *;
         using iterator_category = std::random_access_iterator_tag;
         using difference_type   = ptrdiff_t;

      protected:
         pointer ptr = nullptr;
         size_t  index = 0;

         bool compatible(self_type const & other) const noexcept
         {
            return ptr == other.ptr;
         }

      public:
         constexpr explicit uuid_const_iterator(pointer ptr, size_t const index) :
            ptr(ptr), index(index)
         {
         }

         uuid_const_iterator(uuid_const_iterator const & o) = default;
         uuid_const_iterator& operator=(uuid_const_iterator const & o) = default;
         ~uuid_const_iterator() = default;

         self_type & operator++ ()
         {
            if (index >= 16)
               throw std::out_of_range("Iterator cannot be incremented past the end of the data.");
            ++index;
            return *this;
         }

         self_type operator++ (int)
         {
            self_type tmp = *this;
            ++*this;
            return tmp;
         }

         bool operator== (self_type const & other) const
         {
            assert(compatible(other));
            return index == other.index;
         }

         bool operator!= (self_type const & other) const
         {
            return !(*this == other);
         }

         reference operator* () const
         {
            if (ptr == nullptr)
               throw std::bad_function_call();
            return *(ptr + index);
         }

         reference operator-> () const
         {
            if (ptr == nullptr)
               throw std::bad_function_call();
            return *(ptr + index);
         }

         uuid_const_iterator() = default;

         self_type & operator--()
         {
            if (index <= 0)
               throw std::out_of_range("Iterator cannot be decremented past the beginning of the data.");
            --index;
            return *this;
         }

         self_type operator--(int)
         {
            self_type tmp = *this;
            --*this;
            return tmp;
         }

         self_type operator+(difference_type offset) const
         {
            self_type tmp = *this;
            return tmp += offset;
         }

         self_type operator-(difference_type offset) const
         {
            self_type tmp = *this;
            return tmp -= offset;
         }

         difference_type operator-(self_type const & other) const
         {
            assert(compatible(other));
            return (index - other.index);
         }

         bool operator<(self_type const & other) const
         {
            assert(compatible(other));
            return index < other.index;
         }

         bool operator>(self_type const & other) const
         {
            return other < *this;
         }

         bool operator<=(self_type const & other) const
         {
            return !(other < *this);
         }

         bool operator>=(self_type const & other) const
         {
            return !(*this < other);
         }

         self_type & operator+=(difference_type const offset)
         {
            if (static_cast<difference_type>(index) + offset < 0 ||
               static_cast<difference_type>(index) + offset > 16)
               throw std::out_of_range("Iterator cannot be incremented outside data bounds.");

            index += offset;
            return *this;
         }

         self_type & operator-=(difference_type const offset)
         {
            return *this += -offset;
         }

         value_type const & operator[](difference_type const offset) const
         {
            return (*(*this + offset));
         }
      };

      using value_type = uint8_t;

   public:
      constexpr uuid() noexcept : data({}) {};

      explicit uuid(gsl::span<value_type, 16> bytes)
      {
         std::copy(std::cbegin(bytes), std::cend(bytes), std::begin(data));
      }

      template<typename ForwardIterator>
      explicit uuid(ForwardIterator first, ForwardIterator last)
      {
         if (std::distance(first, last) == 16)
            std::copy(first, last, std::begin(data));
      }

      constexpr uuid_variant variant() const noexcept
      {
         if ((data[8] & 0x80) == 0x00)
            return uuid_variant::ncs;
         else if ((data[8] & 0xC0) == 0x80)
            return uuid_variant::rfc;
         else if ((data[8] & 0xE0) == 0xC0)
            return uuid_variant::microsoft;
         else
            return uuid_variant::reserved;
      }

      constexpr uuid_version version() const noexcept
      {
         if ((data[6] & 0xF0) == 0x10)
            return uuid_version::time_based;
         else if ((data[6] & 0xF0) == 0x20)
            return uuid_version::dce_security;
         else if ((data[6] & 0xF0) == 0x30)
            return uuid_version::name_based_md5;
         else if ((data[6] & 0xF0) == 0x40)
            return uuid_version::random_number_based;
         else if ((data[6] & 0xF0) == 0x50)
            return uuid_version::name_based_sha1;
         else
            return uuid_version::none;
      }

      constexpr std::size_t size() const noexcept { return 16; }

      constexpr bool is_nil() const noexcept
      {
         for (size_t i = 0; i < data.size(); ++i) if (data[i] != 0) return false;
         return true;
      }

      void swap(uuid & other) noexcept
      {
         data.swap(other.data);
      }

      constexpr uuid_const_iterator begin() const noexcept { return uuid_const_iterator(&data[0], 0); }
      constexpr uuid_const_iterator end() const noexcept { return uuid_const_iterator(&data[0], 16); }

      inline gsl::span<std::byte const, 16> as_bytes() const
      {
         return gsl::span<std::byte const, 16>(reinterpret_cast<std::byte const*>(data.data()), 16);
      }

      template <typename TChar>
      static uuid from_string(TChar const * const str, size_t const size)
      {
         TChar digit = 0;
         bool firstDigit = true;
         int hasBraces = 0;
         size_t index = 0;
         std::array<uint8_t, 16> data{ { 0 } };

         if (str == nullptr || size == 0)
            throw uuid_error{ "Wrong uuid format" };

         if (str[0] == static_cast<TChar>('{'))
            hasBraces = 1;
         if (hasBraces && str[size - 1] != static_cast<TChar>('}'))
            throw uuid_error{ "Wrong uuid format" };

         for (size_t i = hasBraces; i < size - hasBraces; ++i)
         {
            if (str[i] == static_cast<TChar>('-')) continue;

            if (index >= 16 || !detail::is_hex(str[i]))
            {
               throw uuid_error{ "Wrong uuid format" };
            }

            if (firstDigit)
            {
               digit = str[i];
               firstDigit = false;
            }
            else
            {
               data[index++] = detail::hexpair2char(digit, str[i]);
               firstDigit = true;
            }
         }

         if (index < 16)
         {
            throw uuid_error{ "Wrong uuid format" };
         }

         return uuid{ std::cbegin(data), std::cend(data) };
      }

      static uuid from_string(std::string_view str)
      {
         return from_string(str.data(), str.size());
      }

      static uuid from_string(std::wstring_view str)
      {
         return from_string(str.data(), str.size());
      }

   private:
      std::array<value_type, 16> data{ { 0 } };

      friend bool operator==(uuid const & lhs, uuid const & rhs) noexcept;
      friend bool operator<(uuid const & lhs, uuid const & rhs) noexcept;

      template <class Elem, class Traits>
      friend std::basic_ostream<Elem, Traits> & operator<<(std::basic_ostream<Elem, Traits> &s, uuid const & id);
   };

   inline bool operator== (uuid const& lhs, uuid const& rhs) noexcept
   {
      return lhs.data == rhs.data;
   }

   inline bool operator!= (uuid const& lhs, uuid const& rhs) noexcept
   {
      return !(lhs == rhs);
   }

   inline bool operator< (uuid const& lhs, uuid const& rhs) noexcept
   {
      return lhs.data < rhs.data;
   }

   template <class Elem, class Traits>
   std::basic_ostream<Elem, Traits> & operator<<(std::basic_ostream<Elem, Traits> &s, uuid const & id)
   {
      return s << std::hex << std::setfill(static_cast<Elem>('0'))
         << std::setw(2) << (int)id.data[0]
         << std::setw(2) << (int)id.data[1]
         << std::setw(2) << (int)id.data[2]
         << std::setw(2) << (int)id.data[3]
         << '-'
         << std::setw(2) << (int)id.data[4]
         << std::setw(2) << (int)id.data[5]
         << '-'
         << std::setw(2) << (int)id.data[6]
         << std::setw(2) << (int)id.data[7]
         << '-'
         << std::setw(2) << (int)id.data[8]
         << std::setw(2) << (int)id.data[9]
         << '-'
         << std::setw(2) << (int)id.data[10]
         << std::setw(2) << (int)id.data[11]
         << std::setw(2) << (int)id.data[12]
         << std::setw(2) << (int)id.data[13]
         << std::setw(2) << (int)id.data[14]
         << std::setw(2) << (int)id.data[15];
   }

   inline std::string to_string(uuid const & id)
   {
      std::stringstream sstr;
      sstr << id;
      return sstr.str();
   }

   inline std::wstring to_wstring(uuid const & id)
   {
      std::wstringstream sstr;
      sstr << id;
      return sstr.str();
   }

   inline void swap(uuids::uuid & lhs, uuids::uuid & rhs)
   {
      lhs.swap(rhs);
   }

   class uuid_system_generator
   {
   public:
      using result_type = uuid;

      uuid operator()()
      {
#ifdef _WIN32

         GUID newId;
         ::CoCreateGuid(&newId);

         std::array<uint8_t, 16> bytes =
         { {
               (unsigned char)((newId.Data1 >> 24) & 0xFF),
               (unsigned char)((newId.Data1 >> 16) & 0xFF),
               (unsigned char)((newId.Data1 >> 8) & 0xFF),
               (unsigned char)((newId.Data1) & 0xFF),

               (unsigned char)((newId.Data2 >> 8) & 0xFF),
               (unsigned char)((newId.Data2) & 0xFF),

               (unsigned char)((newId.Data3 >> 8) & 0xFF),
               (unsigned char)((newId.Data3) & 0xFF),

               newId.Data4[0],
               newId.Data4[1],
               newId.Data4[2],
               newId.Data4[3],
               newId.Data4[4],
               newId.Data4[5],
               newId.Data4[6],
               newId.Data4[7]
            } };

         return uuid{ std::begin(bytes), std::end(bytes) };

#elif defined(__linux__) || defined(__unix__)

         uuid_t id;
         uuid_generate(id);

         std::array<uint8_t, 16> bytes =
         { {
               id[0],
               id[1],
               id[2],
               id[3],
               id[4],
               id[5],
               id[6],
               id[7],
               id[8],
               id[9],
               id[10],
               id[11],
               id[12],
               id[13],
               id[14],
               id[15]
            } };

         return uuid{ std::begin(bytes), std::end(bytes) };

#elif defined(__APPLE__)
         auto newId = CFUUIDCreate(NULL);
         auto bytes = CFUUIDGetUUIDBytes(newId);
         CFRelease(newId);

         std::array<uint8_t, 16> arrbytes =
         { {
               bytes.byte0,
               bytes.byte1,
               bytes.byte2,
               bytes.byte3,
               bytes.byte4,
               bytes.byte5,
               bytes.byte6,
               bytes.byte7,
               bytes.byte8,
               bytes.byte9,
               bytes.byte10,
               bytes.byte11,
               bytes.byte12,
               bytes.byte13,
               bytes.byte14,
               bytes.byte15
            } };
         return uuid{ std::begin(arrbytes), std::end(arrbytes) };
#elif
         return uuid{};
#endif
      }
   };

   template <typename UniformRandomNumberGenerator>
   class basic_uuid_random_generator
   {
   public:
      using result_type = uuid;

      basic_uuid_random_generator()
         :generator(new UniformRandomNumberGenerator)
      {
         std::random_device rd;
         generator->seed(rd());
      }

      explicit basic_uuid_random_generator(UniformRandomNumberGenerator& gen) :
         generator(&gen, [](auto) {}) {}
      explicit basic_uuid_random_generator(UniformRandomNumberGenerator* gen) :
         generator(gen, [](auto) {}) {}

      uuid operator()()
      {
         uint8_t bytes[16];
         for (int i = 0; i < 16; i += 4)
            *reinterpret_cast<uint32_t*>(bytes + i) = distribution(*generator);

         // variant must be 10xxxxxx
         bytes[8] &= 0xBF;
         bytes[8] |= 0x80;

         // version must be 0100xxxx
         bytes[6] &= 0x4F;
         bytes[6] |= 0x40;

         return uuid{std::begin(bytes), std::end(bytes)};
      }

   private:
      std::uniform_int_distribution<uint32_t>  distribution;
      std::shared_ptr<UniformRandomNumberGenerator> generator;
   };

   using uuid_random_generator = basic_uuid_random_generator<std::mt19937>;

   class uuid_name_generator
   {
   public:
      using result_type = uuid;

      explicit uuid_name_generator(uuid const& namespace_uuid) noexcept
         : nsuuid(namespace_uuid)
      {}

      uuid operator()(std::string_view name)
      {
         reset();
         process_characters(name.data(), name.size());
         return make_uuid();
      }

      uuid operator()(std::wstring_view name)
      {
         reset();
         process_characters(name.data(), name.size());
         return make_uuid();
      }

   private:
      void reset()
      {
         hasher.reset();
         uint8_t bytes[16];
         std::copy(std::begin(nsuuid), std::end(nsuuid), bytes);
         hasher.process_bytes(bytes, 16);
      }

      template <typename char_type,
                typename = std::enable_if_t<std::is_integral<char_type>::value>>
      void process_characters(char_type const * const characters, size_t const count)
      {
         for (size_t i = 0; i < count; i++)
         {
            uint32_t c = characters[i];
            hasher.process_byte(static_cast<unsigned char>((c >> 0) & 0xFF));
            hasher.process_byte(static_cast<unsigned char>((c >> 8) & 0xFF));
            hasher.process_byte(static_cast<unsigned char>((c >> 16) & 0xFF));
            hasher.process_byte(static_cast<unsigned char>((c >> 24) & 0xFF));
         }
      }

      void process_characters(const char * const characters, size_t const count)
      {
         hasher.process_bytes(characters, count);
      }

      uuid make_uuid()
      {
         detail::sha1::digest8_t digest;
         hasher.get_digest_bytes(digest);

         // variant must be 0b10xxxxxx
         digest[8] &= 0xBF;
         digest[8] |= 0x80;

         // version must be 0b0101xxxx
         digest[6] &= 0x5F;
         digest[6] |= 0x50;

         return uuid{ digest, digest + 16 };
      }

   private:
      uuid nsuuid;
      detail::sha1 hasher;
   };
}

namespace std
{
   template <>
   struct hash<uuids::uuid>
   {
      using argument_type = uuids::uuid;
      using result_type   = std::size_t;

      result_type operator()(argument_type const &uuid) const
      {
         std::hash<std::string> hasher;
         return static_cast<result_type>(hasher(uuids::to_string(uuid)));
      }
   };
}

max_stars_repo_name: cogment/cogment-orchestrator

max_issues_repo_name: cogment/cogment-orchestrator

max_forks_repo_name: cogment/cogment-orchestrator


text:
#pragma once

#include <boost/variant.hpp>
//#include <gsl/span>
#include <gsl/gsl>
#include <boost/system/system_error.hpp>
#include <boost/asio/read.hpp>
#include <boost/asio/write.hpp>

namespace mothbus
{
	template<class ...T>
	using variant = boost::variant<T...>;

	template<class T>//, std::ptrdiff_t Extent=gsl::dynamic_extent>
	using span = gsl::span<T>;//, Extent>;

	using byte = gsl::byte;

	template <typename SyncWriteStream, typename ConstBufferSequence>
	inline std::size_t write(SyncWriteStream& s, const ConstBufferSequence& buffers)
	{
		return boost::asio::write(s, buffers);
	}

	template <typename SyncWriteStream, typename ConstBufferSequence, typename ErrorCode>
	inline std::size_t write(SyncWriteStream& s, const ConstBufferSequence& buffers, ErrorCode& ec)
	{
		return boost::asio::write(s, buffers, ec);
	}

	template <typename SyncReadStream, typename MutableBufferSequence>
	inline std::size_t read(SyncReadStream& s, const MutableBufferSequence& buffers)
	{
		return boost::asio::read(s, buffers);
	}

	template <typename SyncReadStream, typename MutableBufferSequence, typename ErrorCode>
	inline std::size_t read(SyncReadStream& s, const MutableBufferSequence& buffers,  ErrorCode& ec)
	{
		return boost::asio::read(s, buffers, ec);
	}

	template <typename SyncReadStream, typename MutableBufferSequence, typename ReadHandler>
	inline void async_read(SyncReadStream& s, const MutableBufferSequence& buffers, ReadHandler&& handler)
	{
		boost::asio::async_read(s, buffers, std::forward<ReadHandler>(handler));
	}

	class modbus_exception// : public std::runtime_error
	{
	public:
		modbus_exception(int error_code):
		error_code(error_code)
		{}

		int error_code;
	};
}

max_stars_repo_name: ChrisBFX/mothbus

max_issues_repo_name: ChrisBFX/mothbus

max_forks_repo_name: ChrisBFX/mothbus


text:
#pragma once

#include <gsl/gsl>
#include "halley/utils/utils.h"
#include "halley/data_structures/vector.h"
#include "halley/file/path.h"
#include "halley/data_structures/maybe.h"

namespace Halley {
	class String;
	class Path;

	class FileSystem
	{
	public:
		static bool exists(const Path& p);
		static bool createDir(const Path& p);
		static bool createParentDir(const Path& p);

		static int64_t getLastWriteTime(const Path& p);
		static bool isFile(const Path& p);
		static bool isDirectory(const Path& p);

		static void copyFile(const Path& src, const Path& dst);
		static bool remove(const Path& path);

		static void writeFile(const Path& path, gsl::span<const gsl::byte> data);
		static void writeFile(const Path& path, const Bytes& data);
		static void writeFile(const Path& path, const String& data);
		static Bytes readFile(const Path& path);

		static Vector<Path> enumerateDirectory(const Path& path);

		static Path getRelative(const Path& path, const Path& parentPath);
		static Path getAbsolute(const Path& path);

		static size_t fileSize(const Path& path);

		static Path getTemporaryPath();

		static int runCommand(const String& command);
	};

	class ScopedTemporaryFile final {
	public:
		ScopedTemporaryFile();
		ScopedTemporaryFile(const String& extension);
		~ScopedTemporaryFile();

		ScopedTemporaryFile(const ScopedTemporaryFile& other) = delete;
		ScopedTemporaryFile(ScopedTemporaryFile&& other) noexcept;
		ScopedTemporaryFile& operator=(const ScopedTemporaryFile& other) = delete;
		ScopedTemporaryFile& operator=(ScopedTemporaryFile&& other) noexcept;

		const Path& getPath() const;

	private:
		std::optional<Path> path;
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
/**
 */
#ifndef _INOUT_APER_H

#define _INOUT_APER_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "aXe_errors.h"
#include "spc_cfg.h"
#include "spc_trace_functions.h"

#define APER_MAXLINE 14

extern int
nbeams_from_char_array2 (char **apers, int num);

extern gsl_vector_int *
nbeams_from_char_array (char **apers, int num);

extern int
object_list_to_file (object * const *oblist, char *filename,
                     int leaveout_ignored);

extern int
get_beam_from_aper_file (char *filename, int aperID, int beamID,beam * b);

extern gsl_vector_int *
aper_file_aperlist (char *filename);

extern int
aper_file_apernum (char *filename);

extern object *
get_aperture_from_aper_file (char *filename, int aperID);

extern object **
file_to_object_list (char filename[], observation * obs);

extern char **
return_next_aperture(FILE *input);

extern object **
file_to_object_list_seq (char filename[], observation * obs);

extern int
find_object_in_object_list(object **oblist, const int ID);

extern beam
find_beam_in_object_list(object **oblist, const int objID,
                         const int beamID);

extern beam *
find_beamptr_in_object_list(object **oblist, const int objID,
                            const int beamID);

extern void
refurbish_object_list(object **oblist, const int new_default,
                      const int old_value, const int new_value);

extern int
object_list_size(object **oblist);

extern int
get_beamspec_size(object **oblist);
#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
//
// Copyright (c) 2017 Michele Segata <msegata@disi.unitn.it>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 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 Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program.  If not, see http://www.gnu.org/licenses/.
//

#ifndef MATRIX_UTILS_H
#define MATRIX_UTILS_H

#include "Array.hh"
#include <vector>
#ifdef ENABLE_DISCRETIZATION
#include <gsl/gsl_matrix.h>
#endif

/**
 * Copies the values of a matrix over a portion of another
 * @param dst destination matrix
 * @param src source matrix
 * @param row row index of the destination matrix
 * @param col column index of the destination matrix
 */
void set_submatrix(Matrix<double> &dst, const Matrix<double> &src, int row,
                   int col);

/**
 * Copies the values of a vector over a portion of another
 * @param dst destination vector
 * @param src source vector
 * @param pos index of the destination vector
 */
void set_subvector(Vector<double> &dst, const Vector<double> &src, int pos);

/**
 * Generates a diagonal matrix by using a vector as the diagonal
 * @param src vector to use as the diagonal
 * @return the diagonal matrix
 */
Matrix<double> diag(const Vector<double> &src);

/**
 * Returns a diagonal square matrix
 * @param n square matrix size
 * @param v value to put on the diagonal
 * @return the diagonal matrix
 */
Matrix<double> diag(int n, double v = 1);

/**
 * Returns an identity matrix
 * @param n square matrix size
 * @return the identity matrix
 */
Matrix<double> identity(int n);

/**
 * Given a list of comma separated values, returns a vector
 * @param v string of comma separated values
 * @param length number of elements in the string
 * @return the vector including the parsed values
 */
Vector<double> get_vector(const std::string &v, int length);

/**
 * Returns a vector of all zeros
 * @param length length of the vector
 * @return a zero vector of the specified length
 */
Vector<double> zero_vector(int length);

/**
 * Given a string of comma separated values, returns a matrix of the parsed
 * values
 * @param m string of comma separated values. the number of elements must be
 * nrow * ncol
 * @param nrow number of rows
 * @param ncol number of columns
 * @return the matrix of parsed values
 */
Matrix<double> get_matrix(const std::string &m, int nrow, int ncol);

/**
 * Matrix transpose
 * @param m the matrix to be transposed
 * @return the transposed matrix
 */
Matrix<double> transpose(const Matrix<double> &m);

/**
 * Given a matrix of a particular size, changes the size of the matrix
 * without touching the existing values. New matrix values are initialized
 * with zeros
 * @param m matrix to be enlarged in place
 * @param r new number of rows
 * @param c new number of columns
 */
void extend_matrix(Matrix<double> &m, int r, int c);

/**
 * Given a vector of a particular size, changes the size of the vector
 * without touching the existing values. New vector values are initialized
 * with zeros
 * @param v vector to be enlarged in place
 * @param n new size
 */
void extend_vector(Vector<double> &v, int n);

/**
 * Merges two matrices together by appending the rows of the second matrix to
 * the rows of the first
 * @param a first matrix
 * @param b second matrix
 * @return a matrix where a and b are merged by rows
 */
Matrix<double> bind_matrices(const Matrix<double> &a, const Matrix<double> &b);

/**
 * Merges two vectors together
 * @param a first vector
 * @param b second vector
 * @return a vector composed by a and b
 */
Vector<double> bind_vectors(const Vector<double> &a, const Vector<double> &b);

/**
 * Performs the product of two matrices
 * @param a first matrix
 * @param b second matrix
 * @return the matrix product between a and b if the matrices are compatible,
 * an empty matrix otherwise
 */
Matrix<double> multiply(const Matrix<double> &a, const Matrix<double> &b);

/**
 * Performs the product of a matrix and a vector
 * @param a matrix
 * @param b vector
 * @return the product between a and b if the matrices are compatible,
 * an empty matrix otherwise
 */
Matrix<double> multiply(const Matrix<double> &a, const Vector<double> &b);

/**
 * Performs the product of a matrix and a vector
 * @param a matrix
 * @param b vector
 * @return the product between a and b if the matrices are compatible,
 * an empty matrix otherwise
 */
Matrix<double> multiply(const Matrix<double> &a, const double *b);

/**
 * Given a matrix, the function computes the powers of such a matrix and
 * stores it in a vector of matrices, where the index is the exponent. The
 * results is thus a vector v[i] = A^i for i = 0 ... n
 * @param a input matrix
 * @param n maximum exponent
 * @return the vector of matrices A^i
 */
std::vector<Matrix<double> > get_powers(const Matrix<double> &a, int n);

/**
 * Pretty prints a matrix to stdout
 * @param a the matrix to print
 * @param name optional name to print
 */
void pretty_print_matrix(const Matrix<double_t> &a, const std::string &name="");

/**
 * Copies a vector inside another. If the sizes don't match, then the source
 * vector is repeated multiple time in the destination vector. For example,
 * if the source vector is [1, 2, 3] and the destination vector has length 8,
 * then the destination will be filled with [1, 2, 3, 1, 2, 3, 1, 2]
 * @param dst destination vector
 * @param src source vector
 */
void copy_vector(Vector<double>& dst, const Vector<double> &src);

/**
 * Copies a vector inside another. The size of the destination vector must be
 * larger or equal than the one of the source vector
 * @param dst destination vector
 * @param src source vector
 */
void copy_vector(double *dst, const Vector<double> &src);

/**
 * Copies a vector inside another. The size of the destination vector must be
 * larger or equal than the one of the source vector
 * @param dst destination vector
 * @param src source vector
 */
void copy_vector(double *dst, const Matrix<double> &src);

/**
 * Returns a portion of a vector
 * @param src source vector
 * @param from start index
 * @param length elements to copy
 * @return the subvector of elements going from src[from] to src[from+length-1]
 */
Vector<double> subvector(const Vector<double> &src, int from, int length);

#ifdef ENABLE_DISCRETIZATION

/**
 * Converts a gsl matrix into a Matrix<double> object
 * @param m pointer to gsl_matrix
 * @return converted matrix
 */
Matrix<double> from_gsl(const gsl_matrix *m);

/**
 * Converts a Matrix<double> object into a gsl matrix. This function allocates
 * the memory for the gsl matrix and it is a duty of the user to free the
 * memory via the gsl_matrix_free() method
 * @param m matrix object
 * @return converted matrix
 */
gsl_matrix *to_gsl(const Matrix<double> &m);

#endif

#endif //MATRIX_UTILS_H

max_stars_repo_name: michele-segata/mpclib

max_issues_repo_name: michele-segata/mpclib

max_forks_repo_name: michele-segata/mpclib


text:
/*
   Copyright [2017-2021] [IBM Corporation]
   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at
       http://www.apache.org/licenses/LICENSE-2.0
   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.
*/

#ifndef MCAS_HSTORE_HEAP_MM_EPHEMERAL_H
#define MCAS_HSTORE_HEAP_MM_EPHEMERAL_H

#include "heap_ephemeral.h"

#include "hstore_config.h"
#include "histogram_log2.h"
#include "hop_hash_log.h"
#include <mm_plugin_itf.h>
#include "persistent.h"

#include <ccpm/interfaces.h> /* ownership_callback, (IHeap_expandable, region_vector_t) */
#include <common/byte_span.h>
#include <common/string_view.h>
#include <nupm/region_descriptor.h>
#include <gsl/span>

#include <algorithm> /* min, swap */
#include <cstddef> /* size_t */
#include <vector>

struct heap_mm_ephemeral
	: private heap_ephemeral
{
private:
	using byte_span = common::byte_span;

	/* Although the mm heap tracks managed regions, its definition of a managed
	 * region probably differs from what hstore needs:
	 *
	 * The hstore managed region must include the hstore metadata area, as that
	 * must be include in the "maanged regions" returned through the kvstore
	 * interface.
	 *
	 * The mm heap managed region cannot include the hstore metadata area, as
	 * that would allow the allocator to use the area for memory allocation.
	 *
	 * Therefore, heap_mr_ephemeral keeps its own copy of the managed regions.
	 */
	nupm::region_descriptor _managed_regions;
protected:
	using heap_ephemeral::_alloc_mutex;

protected:
	using hist_type = util::histogram_log2<std::size_t>;
	hist_type _hist_alloc;
	hist_type _hist_inject;
	hist_type _hist_free;
private:
	/* Rca_LB seems not to allocate at or above about 2GiB. Limit reporting to 16 GiB. */
	static constexpr unsigned hist_report_upper_bound = 24U;
	static constexpr unsigned log_min_alignment = 3U; /* log (sizeof(void *)) */
	static_assert(sizeof(void *) == 1U << log_min_alignment, "log_min_alignment does not match sizeof(void *)");

protected:
	virtual void add_managed_region_to_heap(byte_span r_heap) = 0;
public:
	template <bool B>
		void write_hist(const byte_span pool_) const
		{
			static bool suppress = false;
			if ( ! suppress )
			{
				hop_hash_log<B>::write(LOG_LOCATION, "pool ", ::base(pool_));
				std::size_t lower_bound = 0;
				auto limit = std::min(std::size_t(hist_report_upper_bound), _hist_alloc.data().size());
				for ( unsigned i = log_min_alignment; i != limit; ++i )
				{
					const std::size_t upper_bound = 1ULL << i;
					hop_hash_log<B>::write(LOG_LOCATION
						, "[", lower_bound, "..", upper_bound, "): "
						, _hist_alloc.data()[i], " ", _hist_inject.data()[i], " ", _hist_free.data()[i]
						, " "
					);
					lower_bound = upper_bound;
				}
				suppress = true;
			}
		}

	using common::log_source::debug_level;
	explicit heap_mm_ephemeral(
		unsigned debug_level
		, nupm::region_descriptor managed_regions
	);
	virtual ~heap_mm_ephemeral() {}

	virtual std::size_t allocated() const = 0;
	virtual std::size_t capacity() const = 0;

	virtual void allocate(persistent_t<void *> &p, std::size_t sz, std::size_t alignment) = 0;
	virtual std::size_t free(persistent_t<void *> &p_, std::size_t sz_) = 0;
	virtual void free_tracked(const void *p, std::size_t sz) = 0;

	nupm::region_descriptor get_managed_regions() const { return _managed_regions; }
	nupm::region_descriptor set_managed_regions(nupm::region_descriptor n)
	{
		using std::swap;
		swap(n, _managed_regions);
		return n;
	}
	void add_managed_region(byte_span r_full, byte_span r_heap);
	void reconstitute_managed_region(
		const byte_span r_full
		, const byte_span r_heap
		, ccpm::ownership_callback_t f
	);
	virtual void reconstitute_managed_region_to_heap(byte_span r_heap, ccpm::ownership_callback_t f) = 0;
	virtual bool is_crash_consistent() const = 0;
	virtual bool can_reconstitute() const = 0;
};

#endif

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas