text:
#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include "regional.h"

void init_weights(t_stochasticity_grid *g) {
    /* The weights will be normalized s.t. the sum of squared weights is one */
    g->weights[0] = 1.0;
    double sum = 1.0;
    for (int i = 1; i<g->depth; i++) {
        double val = g->weights[i-1] / g->weight_factor;
        g->weights[i] = val;
        sum += val*val;
    }
    double N = sqrt(sum);
    for (int i = 0; i < g->depth; i++) {
        g->weights[i] /= N;
    }
}

/*
 'arr' should point to a double array of size at least height*width.
*/
t_stochasticity_grid *create_stochasticity_grid(double *arr, int height, int width, double weight_factor, double mean, double stddev, long int seed) {
    t_stochasticity_grid *g = malloc( sizeof(t_stochasticity_grid) );
    g->arr = arr;
    g->height = height;
    g->width = width;

    /* It is best to use array dimensions of power of 2 */

    int d = height > width ? height : width; // Take the max
    g->depth = ceil(log2(d)) + 1;
    g->stddev = stddev;
    g->mean = mean;
    g->weight_factor = weight_factor;

    g->rng = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(g->rng, seed);

    g->weights = malloc( sizeof(double) * g->depth );
    init_weights(g);

    return g;
}

void free_stochasticity_grid(t_stochasticity_grid *g) {
    gsl_rng_free(g->rng);
    free(g->weights);
    g->weights = NULL;
    g->arr = NULL;
    free(g);
}

static inline double gaussian(const gsl_rng *rng, double stddev) {
    return gsl_ran_gaussian_ziggurat(rng, stddev);
}

void fill_quadrant(const t_stochasticity_grid *g, int y0, int x0, int y1, int x1, int level) {
    double weight = g->weights[level];
    double rval = gaussian(g->rng, g->stddev);
    for(int i = y0; i < y1; i++) {
        for(int j = x0; j < x1; j++) {
            int offset = i * g->width + j;
            double val = rval*weight;
            g->arr[offset] += val;
        }
    }
}

void recursive_division(const t_stochasticity_grid *g, int level, int y0, int x0, int y1, int x1) {
    if (level >= g->depth) return;
    if (x1-x0 == 1 && y1-y0 == 1) return; // Skip 1x1 grids

    if (y0 >= y1) return;
    if (x0 >= x1) return;
    //printf("Depth: %i. (%i, %i, %i, %i)\n", level, y0, x0, y1, x1);

    fill_quadrant(g, y0, x0, y1, x1, level);

    if (level + 1 == g->depth) return; // Avoid unnecessary recursion; skip level with 1x1 grids

    int ydist = (y1-y0)/2.0 + 0.5;
    int xdist = (x1-x0)/2.0 + 0.5;
    int ymid = y0+ydist;
    int xmid = x0+xdist;

    /* Upper left */
    recursive_division(g, level+1, y0, x0, ymid, xmid);

    /* Upper right */
    recursive_division(g, level+1, ymid, x0, y1, xmid);

    /* Bottom left */
    recursive_division(g, level+1, y0, xmid, ymid, x1);

    /* Bottom right */
    recursive_division(g, level+1, ymid, xmid, y1, x1);
}

/**
 Generate spatially autocorrelated coefficients using a recursive quadtree subdivision.
 The pointer 'arr' is a 2-dimensional output array of size 'height' * 'width'.
 The variable 'variance' gives the variance of the log-norm distribution with mean and cut-off at 1.0.
*/
void generate_stochasticity(t_stochasticity_grid *g) {
    /* Initialize the array by generating the last level (i.e., 1x1 grids) */
    for(int i = 0; i < g->height; i++) {
        for(int j = 0; j < g->width; j++) {
            double weight = g->weights[g->depth-1];
            g->arr[i*g->width + j] = weight*gaussian(g->rng, g->stddev);
        }
    }

    recursive_division(g, 0, 0, 0, g->height, g->width);

    for(int i = 0; i < g->height; i++) {
        for(int j = 0; j < g->width; j++) {
            double y = g->arr[i*g->width + j];
            y = exp(y + g->mean);
            if (y > 1.0) y = 1.0;
            g->arr[i*g->width + j] = y;
        }
    }

}

ext: c

max_stars_repo_name: rybicki/corridor-spom

max_issues_repo_name: rybicki/corridor-spom

max_forks_repo_name: rybicki/corridor-spom


text:
#pragma once

#if SEAL_COMPILER == SEAL_COMPILER_GCC

// We require GCC >= 6
#if (__GNUC__ < 6) || not defined(__cplusplus)
#error "SEAL requires __GNUC__ >= 6"
#endif

// Read in config.h
#include "config.h"

// Are we using MSGSL?
#ifdef SEAL_USE_MSGSL
#include <gsl/gsl>
#endif

// Are intrinsics enabled?
#ifdef SEAL_USE_INTRIN
#include <x86intrin.h>

#ifdef SEAL_USE___BUILTIN_CLZLL
#define SEAL_MSB_INDEX_UINT64(result, value) {                                      \
    *result = 63 - __builtin_clzll(value);                                          \
}
#endif

#ifdef SEAL_USE___INT128
#define SEAL_MULTIPLY_UINT64_HW64(operand1, operand2, hw64) {                       \
    *hw64 = static_cast<uint64_t>((static_cast<unsigned __int128>(operand1)         \
            * static_cast<unsigned __int128>(operand2)) >> 64);                     \
}

#define SEAL_MULTIPLY_UINT64(operand1, operand2, result128) {                       \
    unsigned __int128 product = static_cast<unsigned __int128>(operand1) * operand2;\
    result128[0] = static_cast<uint64_t>(product);                                  \
    result128[1] = product >> 64;                                                   \
}
#endif

#ifdef SEAL_USE__ADDCARRY_U64
#define SEAL_ADD_CARRY_UINT64(operand1, operand2, carry, result) _addcarry_u64(     \
    carry,                                                                          \
    static_cast<unsigned long long>(operand1),                                      \
    static_cast<unsigned long long>(operand2),                                      \
    reinterpret_cast<unsigned long long*>(result))
#endif

#ifdef SEAL_USE__SUBBORROW_U64
#if ((__GNUC__ == 7) && (__GNUC_MINOR__ >= 2)) || (__GNUC__  >= 8)
// The inverted arguments problem was fixed in GCC-7.2
// (https://patchwork.ozlabs.org/patch/784309/)
#define SEAL_SUB_BORROW_UINT64(operand1, operand2, borrow, result) _subborrow_u64(  \
    borrow,                                                                         \
    static_cast<unsigned long long>(operand1),                                      \
    static_cast<unsigned long long>(operand2),                                      \
    reinterpret_cast<unsigned long long*>(result))
#else
// Warning: Note the inverted order of operand1 and operand2
#define SEAL_SUB_BORROW_UINT64(operand1, operand2, borrow, result) _subborrow_u64(  \
    borrow,                                                                         \
    static_cast<unsigned long long>(operand2),                                      \
    static_cast<unsigned long long>(operand1),                                      \
    reinterpret_cast<unsigned long long*>(result))
#endif //(__GNUC__ == 7) && (__GNUC_MINOR__ >= 2)
#endif

#endif //SEAL_USE_INTRIN

#endif

ext: h


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

#include <config.h>
#include <gsl/gsl_fft.h>

#include "complex_internal.h"
#include "bitreverse.h"

static int
FUNCTION(fft_complex,bitreverse_order) (BASE data[],
                                        const size_t stride,
                                        const size_t n,
                                        size_t logn)
{
  /* This is the Goldrader bit-reversal algorithm */

  size_t i;
  size_t j = 0;

  logn = 0 ; /* not needed for this algorithm */

  for (i = 0; i < n - 1; i++)
    {
      size_t k = n / 2 ;

      if (i < j)
        {
          const BASE tmp_real = REAL(data,stride,i);
          const BASE tmp_imag = IMAG(data,stride,i);
          REAL(data,stride,i) = REAL(data,stride,j);
          IMAG(data,stride,i) = IMAG(data,stride,j);
          REAL(data,stride,j) = tmp_real;
          IMAG(data,stride,j) = tmp_imag;
        }

      while (k <= j)
        {
          j = j - k ;
          k = k / 2 ;
        }

      j += k ;
    }

  return 0;
}

static int
FUNCTION(fft_real,bitreverse_order) (BASE data[],
                                const size_t stride,
                                const size_t n,
                                size_t logn)
{
  /* This is the Goldrader bit-reversal algorithm */

  size_t i;
  size_t j = 0;

  logn = 0 ; /* not needed for this algorithm */

  for (i = 0; i < n - 1; i++)
    {
      size_t k = n / 2 ;

      if (i < j)
        {
          const BASE tmp = VECTOR(data,stride,i);
          VECTOR(data,stride,i) = VECTOR(data,stride,j);
          VECTOR(data,stride,j) = tmp;
        }

      while (k <= j)
        {
          j = j - k ;
          k = k / 2 ;
        }

      j += k ;
    }

  return 0;
}

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 "schedule.h"
#include <boost/random/uniform_01.hpp>
#include <boost/random/uniform_int_distribution.hpp>
#include <gsl/gsl-lite.hpp>
#include <pcg_random.hpp>
#include <random>

namespace angonoka::stun {
using RandomEngine = pcg32;

/**
    Miscellaneous random number generators.
*/
class RandomUtils {
public:
    /**
        Default constructor.
    */
    RandomUtils();

    /**
        Constructor with a fixed PRNG seed.

        @param seed Random engine seed
    */
    RandomUtils(gsl::index seed);

    /**
        Uniformally distributed real value between 0 and 1.

        @return Random number
    */
    float uniform_01() noexcept;

    /**
        Uniformally distributed discrete value between 0 and max.

        @param max Maximum value

        @return Random number
    */
    int16 uniform_int(int16 max) noexcept;

private:
    RandomEngine generator{
        pcg_extras::seed_seq_from<std::random_device>{}};
    boost::random::uniform_01<float> uniform_01_;
    boost::random::uniform_int_distribution<std::int_fast16_t>
        uniform_int_;
};
} // namespace angonoka::stun

max_stars_repo_name: coffee-lord/angonoka

max_issues_repo_name: coffee-lord/angonoka

max_forks_repo_name: coffee-lord/angonoka


text:
#undef FP_T
#undef FP_ID
#undef VECTOR_T
#undef VECTOR_ID
#undef MATRIX_T
#undef MATRIX_ID

#ifdef GSL_FLOAT
#define MATLAB_NAMESPACE matlab_float
#define BCT_NAMESPACE bct_float
#define FP_T float
#define FP_ID(id) id##_##float
#define VECTOR_T gsl_vector_float
#define VECTOR_ID(id) gsl_vector_float##_##id
#define MATRIX_T gsl_matrix_float
#define MATRIX_ID(id) gsl_matrix_float##_##id
#include <gsl/gsl_matrix_float.h>
#include <gsl/gsl_vector_float.h>
#endif

#ifdef GSL_DOUBLE
#define MATLAB_NAMESPACE matlab
#define BCT_NAMESPACE bct
#define FP_T double
#define FP_ID(id) id##_##double
#define VECTOR_T gsl_vector
#define VECTOR_ID(id) gsl_vector##_##id
#define MATRIX_T gsl_matrix
#define MATRIX_ID(id) gsl_matrix##_##id
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#endif

#ifdef GSL_LONG_DOUBLE
#define MATLAB_NAMESPACE matlab_long_double
#define BCT_NAMESPACE bct_long_double
#define FP_T long double
#define FP_ID(id) id##_##long_double
#define VECTOR_T gsl_vector_long_double
#define VECTOR_ID(id) gsl_vector_long_double##_##id
#define MATRIX_T gsl_matrix_long_double
#define MATRIX_ID(id) gsl_matrix_long_double##_##id
#include <gsl/gsl_matrix_long_double.h>
#include <gsl/gsl_vector_long_double.h>
#endif

max_stars_repo_name: devuci/bct-cpp

max_issues_repo_name: devuci/bct-cpp

max_forks_repo_name: devuci/bct-cpp


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

#ifndef __GSL_VECTOR_LONG_H__
#define __GSL_VECTOR_LONG_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_long.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

typedef struct
{
  size_t size;
  size_t stride;
  long *data;
  gsl_block_long *block;
  int owner;
}
gsl_vector_long;

typedef struct
{
  gsl_vector_long vector;
} _gsl_vector_long_view;

typedef _gsl_vector_long_view gsl_vector_long_view;

typedef struct
{
  gsl_vector_long vector;
} _gsl_vector_long_const_view;

typedef const _gsl_vector_long_const_view gsl_vector_long_const_view;

/* Allocation */

GSL_EXPORT gsl_vector_long *gsl_vector_long_alloc (const size_t n);
GSL_EXPORT gsl_vector_long *gsl_vector_long_calloc (const size_t n);

GSL_EXPORT gsl_vector_long *gsl_vector_long_alloc_from_block (gsl_block_long * b,
                                                                const size_t offset,
                                                                const size_t n,
                                                                const size_t stride);

GSL_EXPORT gsl_vector_long *gsl_vector_long_alloc_from_vector (gsl_vector_long * v,
                                                                 const size_t offset,
                                                                 const size_t n,
                                                                 const size_t stride);

GSL_EXPORT void gsl_vector_long_free (gsl_vector_long * v);

/* Views */

GSL_EXPORT
_gsl_vector_long_view
gsl_vector_long_view_array (long *v, size_t n);

GSL_EXPORT
_gsl_vector_long_view
gsl_vector_long_view_array_with_stride (long *base,
                                         size_t stride,
                                         size_t n);

GSL_EXPORT
_gsl_vector_long_const_view
gsl_vector_long_const_view_array (const long *v, size_t n);

GSL_EXPORT
_gsl_vector_long_const_view
gsl_vector_long_const_view_array_with_stride (const long *base,
                                               size_t stride,
                                               size_t n);

GSL_EXPORT
_gsl_vector_long_view
gsl_vector_long_subvector (gsl_vector_long *v,
                            size_t i,
                            size_t n);

GSL_EXPORT
_gsl_vector_long_view
gsl_vector_long_subvector_with_stride (gsl_vector_long *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_EXPORT
_gsl_vector_long_const_view
gsl_vector_long_const_subvector (const gsl_vector_long *v,
                                  size_t i,
                                  size_t n);

GSL_EXPORT
_gsl_vector_long_const_view
gsl_vector_long_const_subvector_with_stride (const gsl_vector_long *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_EXPORT long gsl_vector_long_get (const gsl_vector_long * v, const size_t i);
GSL_EXPORT void gsl_vector_long_set (gsl_vector_long * v, const size_t i, long x);

GSL_EXPORT long *gsl_vector_long_ptr (gsl_vector_long * v, const size_t i);
GSL_EXPORT const long *gsl_vector_long_const_ptr (const gsl_vector_long * v, const size_t i);

GSL_EXPORT void gsl_vector_long_set_zero (gsl_vector_long * v);
GSL_EXPORT void gsl_vector_long_set_all (gsl_vector_long * v, long x);
GSL_EXPORT int gsl_vector_long_set_basis (gsl_vector_long * v, size_t i);

GSL_EXPORT int gsl_vector_long_fread (FILE * stream, gsl_vector_long * v);
GSL_EXPORT int gsl_vector_long_fwrite (FILE * stream, const gsl_vector_long * v);
GSL_EXPORT int gsl_vector_long_fscanf (FILE * stream, gsl_vector_long * v);
GSL_EXPORT int gsl_vector_long_fprintf (FILE * stream, const gsl_vector_long * v,
                                         const char *format);

GSL_EXPORT int gsl_vector_long_memcpy (gsl_vector_long * dest, const gsl_vector_long * src);

GSL_EXPORT int gsl_vector_long_reverse (gsl_vector_long * v);

GSL_EXPORT int gsl_vector_long_swap (gsl_vector_long * v, gsl_vector_long * w);
GSL_EXPORT int gsl_vector_long_swap_elements (gsl_vector_long * v, const size_t i, const size_t j);

GSL_EXPORT long gsl_vector_long_max (const gsl_vector_long * v);
GSL_EXPORT long gsl_vector_long_min (const gsl_vector_long * v);
GSL_EXPORT void gsl_vector_long_minmax (const gsl_vector_long * v, long * min_out, long * max_out);

GSL_EXPORT size_t gsl_vector_long_max_index (const gsl_vector_long * v);
GSL_EXPORT size_t gsl_vector_long_min_index (const gsl_vector_long * v);
GSL_EXPORT void gsl_vector_long_minmax_index (const gsl_vector_long * v, size_t * imin, size_t * imax);

GSL_EXPORT int gsl_vector_long_add (gsl_vector_long * a, const gsl_vector_long * b);
GSL_EXPORT int gsl_vector_long_sub (gsl_vector_long * a, const gsl_vector_long * b);
GSL_EXPORT int gsl_vector_long_mul (gsl_vector_long * a, const gsl_vector_long * b);
GSL_EXPORT int gsl_vector_long_div (gsl_vector_long * a, const gsl_vector_long * b);
GSL_EXPORT int gsl_vector_long_scale (gsl_vector_long * a, const double x);
GSL_EXPORT int gsl_vector_long_add_constant (gsl_vector_long * a, const double x);

GSL_EXPORT int gsl_vector_long_isnull (const gsl_vector_long * v);

#ifdef HAVE_INLINE

extern inline
long
gsl_vector_long_get (const gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_long_set (gsl_vector_long * v, const size_t i, long x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
long *
gsl_vector_long_ptr (gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (long *) (v->data + i * v->stride);
}

extern inline
const long *
gsl_vector_long_const_ptr (const gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const long *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_LONG_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef bfkl_H_INCLUDED
#define bfkl_H_INCLUDED

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "cubature.h"
#include "string.h"
#include "alphas.h"
#include "wf.h"

//number of tabulated points in z and x for BFKL Kernel
#define Nz 94
#define Nx 1201
#define Nhar 24
#define ZMIN (0.15)
#define ZMAX (4.8)
#define XMAX (100.00)

void ReadInBFKL();
double BFKLfunc(double z, double x, int n);
double d2N_BFKL(double pT, double qT, double phi, double yp, double yq, double rts);

struct bfkl_params {
	double pT;
	double qT;
	double phi;
	double x1, x2;
    double dy;
	} ;

#endif

max_stars_repo_name: kdusling/mpc

max_issues_repo_name: kdusling/mpc

max_forks_repo_name: kdusling/mpc


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

#ifndef __GSL_RANDIST_H__
#define __GSL_RANDIST_H__
#include <gsl/gsl_rng.h>
#include <gsl/gsl_types.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

GSL_EXPORT unsigned int gsl_ran_bernoulli (const gsl_rng * r, double p);
GSL_EXPORT double gsl_ran_bernoulli_pdf (const unsigned int k, double p);

GSL_EXPORT double gsl_ran_beta (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_beta_pdf (const double x, const double a, const double b);

GSL_EXPORT unsigned int gsl_ran_binomial (const gsl_rng * r, double p, unsigned int n);
GSL_EXPORT unsigned int gsl_ran_binomial_knuth (const gsl_rng * r, double p, unsigned int n);
GSL_EXPORT unsigned int gsl_ran_binomial_tpe (const gsl_rng * r, double p, unsigned int n);
GSL_EXPORT double gsl_ran_binomial_pdf (const unsigned int k, const double p, const unsigned int n);

GSL_EXPORT double gsl_ran_exponential (const gsl_rng * r, const double mu);
GSL_EXPORT double gsl_ran_exponential_pdf (const double x, const double mu);

GSL_EXPORT double gsl_ran_exppow (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_exppow_pdf (const double x, const double a, const double b);

GSL_EXPORT double gsl_ran_cauchy (const gsl_rng * r, const double a);
GSL_EXPORT double gsl_ran_cauchy_pdf (const double x, const double a);

GSL_EXPORT double gsl_ran_chisq (const gsl_rng * r, const double nu);
GSL_EXPORT double gsl_ran_chisq_pdf (const double x, const double nu);

GSL_EXPORT void gsl_ran_dirichlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
GSL_EXPORT double gsl_ran_dirichlet_pdf (const size_t K, const double alpha[], const double theta[]);
GSL_EXPORT double gsl_ran_dirichlet_lnpdf (const size_t K, const double alpha[], const double theta[]);

GSL_EXPORT double gsl_ran_erlang (const gsl_rng * r, const double a, const double n);
GSL_EXPORT double gsl_ran_erlang_pdf (const double x, const double a, const double n);

GSL_EXPORT double gsl_ran_fdist (const gsl_rng * r, const double nu1, const double nu2);
GSL_EXPORT double gsl_ran_fdist_pdf (const double x, const double nu1, const double nu2);

GSL_EXPORT double gsl_ran_flat (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_flat_pdf (double x, const double a, const double b);

GSL_EXPORT double gsl_ran_gamma (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a);
GSL_EXPORT double gsl_ran_gamma_pdf (const double x, const double a, const double b);
GSL_EXPORT double gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b);

GSL_EXPORT double gsl_ran_gaussian (const gsl_rng * r, const double sigma);
GSL_EXPORT double gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma);
GSL_EXPORT double gsl_ran_gaussian_ziggurat (const gsl_rng * r, const double sigma);
GSL_EXPORT double gsl_ran_gaussian_pdf (const double x, const double sigma);

GSL_EXPORT double gsl_ran_ugaussian (const gsl_rng * r);
GSL_EXPORT double gsl_ran_ugaussian_ratio_method (const gsl_rng * r);
GSL_EXPORT double gsl_ran_ugaussian_pdf (const double x);

GSL_EXPORT double gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma);
GSL_EXPORT double gsl_ran_gaussian_tail_pdf (const double x, const double a, const double sigma);

GSL_EXPORT double gsl_ran_ugaussian_tail (const gsl_rng * r, const double a);
GSL_EXPORT double gsl_ran_ugaussian_tail_pdf (const double x, const double a);

GSL_EXPORT void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y);
GSL_EXPORT double gsl_ran_bivariate_gaussian_pdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);

GSL_EXPORT double gsl_ran_landau (const gsl_rng * r);
GSL_EXPORT double gsl_ran_landau_pdf (const double x);

GSL_EXPORT unsigned int gsl_ran_geometric (const gsl_rng * r, const double p);
GSL_EXPORT double gsl_ran_geometric_pdf (const unsigned int k, const double p);

GSL_EXPORT unsigned int gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2, unsigned int t);
GSL_EXPORT double gsl_ran_hypergeometric_pdf (const unsigned int k, const unsigned int n1, const unsigned int n2, unsigned int t);

GSL_EXPORT double gsl_ran_gumbel1 (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_gumbel1_pdf (const double x, const double a, const double b);

GSL_EXPORT double gsl_ran_gumbel2 (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_gumbel2_pdf (const double x, const double a, const double b);

GSL_EXPORT double gsl_ran_logistic (const gsl_rng * r, const double a);
GSL_EXPORT double gsl_ran_logistic_pdf (const double x, const double a);

GSL_EXPORT double gsl_ran_lognormal (const gsl_rng * r, const double zeta, const double sigma);
GSL_EXPORT double gsl_ran_lognormal_pdf (const double x, const double zeta, const double sigma);

GSL_EXPORT unsigned int gsl_ran_logarithmic (const gsl_rng * r, const double p);
GSL_EXPORT double gsl_ran_logarithmic_pdf (const unsigned int k, const double p);

GSL_EXPORT void gsl_ran_multinomial (const gsl_rng * r, const size_t K,
                                     const unsigned int N, const double p[],
                                     unsigned int n[] );
GSL_EXPORT double gsl_ran_multinomial_pdf (const size_t K,
                                           const double p[], const unsigned int n[] );
GSL_EXPORT double gsl_ran_multinomial_lnpdf (const size_t K,
                                             const double p[], const unsigned int n[] );

GSL_EXPORT unsigned int gsl_ran_negative_binomial (const gsl_rng * r, double p, double n);
GSL_EXPORT double gsl_ran_negative_binomial_pdf (const unsigned int k, const double p, double n);

GSL_EXPORT unsigned int gsl_ran_pascal (const gsl_rng * r, double p, unsigned int n);
GSL_EXPORT double gsl_ran_pascal_pdf (const unsigned int k, const double p, unsigned int n);

GSL_EXPORT double gsl_ran_pareto (const gsl_rng * r, double a, const double b);
GSL_EXPORT double gsl_ran_pareto_pdf (const double x, const double a, const double b);

GSL_EXPORT unsigned int gsl_ran_poisson (const gsl_rng * r, double mu);
GSL_EXPORT void gsl_ran_poisson_array (const gsl_rng * r, size_t n, unsigned int array[],
                                       double mu);
GSL_EXPORT double gsl_ran_poisson_pdf (const unsigned int k, const double mu);

GSL_EXPORT double gsl_ran_rayleigh (const gsl_rng * r, const double sigma);
GSL_EXPORT double gsl_ran_rayleigh_pdf (const double x, const double sigma);

GSL_EXPORT double gsl_ran_rayleigh_tail (const gsl_rng * r, const double a, const double sigma);
GSL_EXPORT double gsl_ran_rayleigh_tail_pdf (const double x, const double a, const double sigma);

GSL_EXPORT double gsl_ran_tdist (const gsl_rng * r, const double nu);
GSL_EXPORT double gsl_ran_tdist_pdf (const double x, const double nu);

GSL_EXPORT double gsl_ran_laplace (const gsl_rng * r, const double a);
GSL_EXPORT double gsl_ran_laplace_pdf (const double x, const double a);

GSL_EXPORT double gsl_ran_levy (const gsl_rng * r, const double c, const double alpha);
GSL_EXPORT double gsl_ran_levy_skew (const gsl_rng * r, const double c, const double alpha, const double beta);

GSL_EXPORT double gsl_ran_weibull (const gsl_rng * r, const double a, const double b);
GSL_EXPORT double gsl_ran_weibull_pdf (const double x, const double a, const double b);

GSL_EXPORT void gsl_ran_dir_2d (const gsl_rng * r, double * x, double * y);
GSL_EXPORT void gsl_ran_dir_2d_trig_method (const gsl_rng * r, double * x, double * y);
GSL_EXPORT void gsl_ran_dir_3d (const gsl_rng * r, double * x, double * y, double * z);
GSL_EXPORT void gsl_ran_dir_nd (const gsl_rng * r, size_t n, double * x);

GSL_EXPORT void gsl_ran_shuffle (const gsl_rng * r, void * base, size_t nmembm, size_t size);
GSL_EXPORT int gsl_ran_choose (const gsl_rng * r, void * dest, size_t k, void * src, size_t n, size_t size) ;
GSL_EXPORT void gsl_ran_sample (const gsl_rng * r, void * dest, size_t k, void * src, size_t n, size_t size) ;

typedef struct {                /* struct for Walker algorithm */
    size_t K;
    size_t *A;
    double *F;
} gsl_ran_discrete_t;

GSL_EXPORT gsl_ran_discrete_t * gsl_ran_discrete_preproc (size_t K, const double *P);
GSL_EXPORT void gsl_ran_discrete_free(gsl_ran_discrete_t *g);
GSL_EXPORT size_t gsl_ran_discrete (const gsl_rng *r, const gsl_ran_discrete_t *g);
GSL_EXPORT double gsl_ran_discrete_pdf (size_t k, const gsl_ran_discrete_t *g);

__END_DECLS

#endif /* __GSL_RANDIST_H__ */


text:
/* sys/ldfrexp.c
 *
 * Copyright (C) 2002, Gert Van den Eynde
 *
 * 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.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>

double
gsl_ldexp (const double x, const int e)
{
  double p2 = pow (2.0, (double)e);
  return x * p2;
}

double
gsl_frexp (const double x, int *e)
{
  if (x == 0.0)
    {
      *e = 0;
      return 0.0;
    }
  else
    {
      double ex = ceil (log (fabs (x)) / M_LN2);
      int ei = (int) ex;
      double f = gsl_ldexp (x, -ei);

      while (fabs (f) >= 1.0)
        {
          ei++;
          f /= 2.0;
        }

      while (fabs (f) < 0.5)
        {
          ei--;
          f *= 2.0;
        }

      *e = ei;
      return f;
    }
}


text:
/* histogram/gsl-histogram.c
 *
 * 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.
 */

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <gsl/gsl_histogram.h>

void
print_help(void)
{
  fprintf (stderr, "Usage: gsl-histogram [-u] xmin xmax [n]\n");
  fprintf (stderr, "Computes a histogram of the data on stdin using n bins from xmin to xmax.\n"
                   "If n is unspecified then bins of integer width are used.\n"
                   "If -u is given, histogram is normalized so that sum of all bins is unity.\n");
}

int
main (int argc, char **argv)
{
  double a = 0.0, b = 1.0;
  size_t n = 10;
  int unit = 0;
  int c;

  while ((c = getopt(argc, argv, "u")) != (-1))
    {
      switch (c)
        {
          case 'u':
            unit = 1;
            break;

          default:
            print_help();
            exit(0);
            break;
        }
    }

  if (argc - optind < 2)
    {
      print_help();
      exit(0);
    }

  a = atof (argv[optind++]);
  b = atof (argv[optind++]);

  if (argc - optind > 0)
    n = atoi (argv[optind++]);
  else
    n = (int)(b - a) ;

  {
    double x;
    gsl_histogram *h = gsl_histogram_alloc (n);

    gsl_histogram_set_ranges_uniform (h, a, b);

    while (fscanf(stdin, "%lg", &x) == 1)
      {
        gsl_histogram_increment(h, x);
      }

#ifdef DISPLAY_STATS
    {
      double mean = gsl_histogram_mean (h);
      double sigma = gsl_histogram_sigma (h);
      fprintf (stdout, "# mean = %g\n", mean);
      fprintf (stdout, "# sigma = %g\n", sigma);
    }
#endif

    /* normalize histogram if needed */
    if (unit)
      {
        double sum = gsl_histogram_sum(h);

        if (sum > 0.0)
          gsl_histogram_scale(h, 1.0 / sum);
      }

    gsl_histogram_fprintf (stdout, h, "%g", "%g");

    gsl_histogram_free (h);
  }

  return 0;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <linux/types.h>
#include <linux/types.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <math.h>

#include "sweeny_dc.h"
#include "timeseries.h"
#include "extract_args.h"

#define BUF_LEN 512
char fn[BUF_LEN];
char verbose = 0;

#define TS_FILE_D  "../time_series/dc/simulation_l%u_q%.4f_b%.4f_c%.4f_s%u.hdf5"
char impl_title[] = "Dynamic Connectivity implementation";
__u32 DX;
__u32 seed;
__u32 cutoff;
__u32 steps;
double beta;
double coupling;
double q;
double rnd_num;

__u32 *num_bonds, *num_cluster, *size_giant;
__u64 *sec_cs_moment,*four_cs_moment;

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

    if(!extractArgs(argc, argv,impl_title))
      return EXIT_FAILURE;
    snprintf(fn,BUF_LEN,TS_FILE_D,DX,q,beta,coupling,seed);
    if(!init_observables())
        return EXIT_FAILURE;
    if(!init_sweeny_dc(q,DX,beta,coupling,cutoff,steps,seed,num_bonds,num_cluster,size_giant,
                sec_cs_moment,four_cs_moment))
        return EXIT_FAILURE;
    if(!simulate_sweeny_dc())
        return EXIT_FAILURE;
    destroy_sweeny_dc();
    save_timeseries();
    destroy_observables();
    return EXIT_SUCCESS;

}

max_issues_repo_name: ernmeel/sweeny


text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <halley/utils/utils.h>
#include <halley/support/exception.h>
#include <gsl/span>
#include <cstdint>

namespace Halley {
	class MT199937AR;

	class Random {
	public:
		static Random& getGlobal();

		Random();
		Random(uint32_t seed);
		Random(uint64_t seed);
		Random(gsl::span<const gsl::byte> data);
		~Random();

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

		int32_t getInt(int32_t min, int32_t max); // [min, max]
		uint32_t getInt(uint32_t min, uint32_t max); // [min, max]
		int64_t getInt(int64_t min, int64_t max); // [min, max]
		uint64_t getInt(uint64_t min, uint64_t max); // [min, max]
		size_t getSizeT(size_t min, size_t max); // [min, max]
		float getFloat(float min, float max); // [min, max)
		double getDouble(double min, double max); // [min, max)

		template <typename T>
		size_t getRandomIndex(const T& vec)
		{
			size_t size = std::end(vec) - std::begin(vec);
			if (size == 0) {
				throw Exception("Can't get random index of empty sequence.", HalleyExceptions::Utils);
			}
			return getSizeT(size_t(0), size_t(size - 1));
		}

		template <typename T>
		auto getRandomElement(T& vec) -> decltype(vec[0])&
		{
			return vec[getRandomIndex(vec)];
		}

		void getBytes(gsl::span<gsl::byte> dst);
		void getBytes(gsl::span<Byte> dst);
		void setSeed(uint32_t seed);
		void setSeed(uint64_t seed);
		void setSeed(gsl::span<const gsl::byte> data);
		void setSeed(gsl::span<Byte> data);

		uint32_t getRawInt();
		float getRawFloat();
		double getRawDouble();

	private:
		std::unique_ptr<MT199937AR> generator;
	};

}

max_stars_repo_name: moonantonio/halley

max_issues_repo_name: moonantonio/halley

max_forks_repo_name: moonantonio/halley


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_sf.h>
#include "adkGSL.h"
#include "AFS.h"
#include "cs.h"
#include "adkCSparse.h"
#include "AFS_ctmc.h"
#include <pthread.h>

void coalMarkovChainTopologyMatrix_pthread(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int *dim1, int *dim2, int start, int stop, int *nnz);
void coalMarkovChainTopologyMatrix_sparse_pthread(afsStateSpace *S,double *topol, int *moveType, int *dim1, int *dim2, int start, int stop, int *nnz);

typedef struct matrixThreadObject
{
	int nnz,start, stop;
	int *dim1,*dim2;
	struct island_lik_params *island_params;
}
matrixThreadObject;

typedef struct matrixThreadObject_sparse
{
	int nnz,start, stop, nnzA;
	int *dim1,*dim2, *moveType;
	double *topol;
	struct im_lik_params *im_params;
}
matrixThreadObject_sparse;

//gonna try to use a global declared in island.c
void *threadDoRows(void* arg){
	matrixThreadObject *obj = (matrixThreadObject*)(arg);

	coalMarkovChainTopologyMatrix_pthread(obj->island_params->stateSpace,obj->island_params->topol, obj->island_params->moveType,
	obj->dim1,obj->dim2,obj->start,obj->stop,&(obj->nnz));
	//obj->nnz -= 1;

	pthread_exit((void* )obj);
}

//sparse matrix version
void *threadDoRows_sparse(void* arg){
	matrixThreadObject_sparse *obj = (matrixThreadObject_sparse*)(arg);

	coalMarkovChainTopologyMatrix_sparse_pthread(obj->im_params->stateSpace,obj->topol, obj->moveType,
	obj->dim1,obj->dim2,obj->start,obj->stop,&(obj->nnz));
	//obj->nnz -= 1;

	pthread_exit((void* )obj);
}

//sparse matrix version on reduced state space
void *threadDoRows_sparse_reduced(void* arg){
	matrixThreadObject_sparse *obj = (matrixThreadObject_sparse*)(arg);

	coalMarkovChainTopologyMatrix_sparse_pthread(obj->im_params->reducedStateSpace,obj->topol, obj->moveType,
	obj->dim1,obj->dim2,obj->start,obj->stop,&(obj->nnzA));
	//obj->nnz -= 1;

	pthread_exit((void* )obj);
}

void coalMarkovChainTopologyMatrix_pthread(afsStateSpace *S,gsl_matrix *topol, gsl_matrix_int *moveType, int *dim1, int *dim2, int start, int stop, int *nnz){
	int i,j,k,steps,x,y;
	double top,bottom;
	afsObject *delta;
	int **activePos1, **activePos2; //first dimension here relates to 2 popns
	int nlin1,nlin2,compat;
	int nonZeroCount, tot, tCount;
	int size = 200;
	nlin1 = nlin2 = 0;
	nonZeroCount = 0;
	if(S->nstates != topol->size1 || S->nstates != moveType->size1){
		fprintf(stderr,"Error: StateSpace and matrices are not of equal size\n");
		exit(1);
	}
	*nnz = 0;
	//arrays for finding positions of entries
	activePos1 = malloc(size*sizeof(int *));
	activePos2 = malloc(size*sizeof(int *));
	for(i=0;i<size;i++){
		activePos1[i] = malloc(2*sizeof(int));
		activePos1[i][0] = activePos1[i][1] = 0;
		activePos2[i] = malloc(2*sizeof(int));
		activePos2[i][0] = activePos2[i][1] = 0;

	}

	tot = S->nstates * S->nstates;
	tCount = 0;
	delta = afsObjectNewFrom(S->states[0]);
	for(i=start;i<stop;i++){
		for(j=0;j<S->nstates;j++){
			gsl_matrix_int_set(moveType,i,j,666);

			//is ith state root state?
			if(S->states[i]->nalleles == 1){
				//set correct absorbing conditions
				if(i==j){
					gsl_matrix_set(topol,i,j,1.0);
					gsl_matrix_int_set(moveType,i,j,-1);
					dim1[nonZeroCount] = i;
					dim2[nonZeroCount] = j;
					nonZeroCount += 1;
				}
				else{
					gsl_matrix_set(topol,i,j,0.0);
					gsl_matrix_int_set(moveType,i,j,-1);
				}
			}
			else{ //ith not root; what is the move?
				steps = S->states[i]->nalleles - S->states[j]->nalleles ;
				if(steps < 2){
				//	delta = afsObjectDelta(S->states[i],S->states[j]);
					afsObjectDeltaPre(S->states[i],S->states[j],delta);
					switch(steps){
						case 0:
						// 0 "steps", so no changes in nalleles between two
						//check counts and positions
						if(abs(matrixSum(delta->popMats[0])) == 1){
							nonZeroEntries(delta->popMats[0], activePos1, &nlin1);
							nonZeroEntries(delta->popMats[1], activePos2, &nlin2);

							//migration?
							if(nlin1 == nlin2 && activePos1[0][0] == activePos2[0][0] &&
								(gsl_matrix_int_min(delta->popMats[0]) >= 0 || gsl_matrix_int_min(delta->popMats[1]) >= 0))
							{

								//migration popn1?
								if(matrixSum(delta->popMats[0]) == 1){
									gsl_matrix_set(topol,i,j,
										(double) gsl_matrix_int_get(S->states[i]->popMats[0],
										activePos1[0][0], activePos1[0][1]) / S->states[i]->aCounts[0]);
									gsl_matrix_int_set(moveType,i,j,2);
									dim1[nonZeroCount] = i;
									dim2[nonZeroCount] = j;
									nonZeroCount += 1;
								}
								else{ //migration popn2
									gsl_matrix_set(topol,i,j,
										(double) gsl_matrix_int_get(S->states[i]->popMats[1],
										activePos2[0][0],activePos2[0][1]) / S->states[i]->aCounts[1]);
									gsl_matrix_int_set(moveType,i,j,3);
									dim1[nonZeroCount] = i;
									dim2[nonZeroCount] = j;
									nonZeroCount += 1;
								}
							}
						}
						break;

						case 1://coalescent?
						//coal popn1?
						if(matrixSum(delta->popMats[0]) == 1 && countNegativeEntries(delta->popMats[1]) == 0 && S->states[j]->aCounts[0] >= 1 ){
							entriesGreaterThan(delta->popMats[0], activePos1, &nlin1,0);//old lineages
							entriesLessThan(delta->popMats[0], activePos2, &nlin2,0);//new lineages
							if(nlin2 != 0 && nlin2 < 2){

								compat=0;
								if(nlin1 == 1){
									if(2*activePos1[0][0] == activePos2[0][0] && 2*activePos1[0][1] == activePos2[0][1])
										compat = 1;
								}
								if(nlin1==2){
									if(activePos1[0][0] + activePos1[1][0] == activePos2[0][0] &&
										activePos1[0][1] + activePos1[1][1] == activePos2[0][1])
										compat = 1;
								}
								if(compat != 1){
									gsl_matrix_set(topol,i,j,0.0);
									gsl_matrix_int_set(moveType,i,j,666);
								}
								else{

									top = 1.0;
									bottom = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1) /2.0 ;
									for(k=0;k<nlin1;k++){
										x = activePos1[k][0];
										y = activePos1[k][1];
										top *= (float) xchoosey(gsl_matrix_int_get(S->states[i]->popMats[0],x,y),
											gsl_matrix_int_get(delta->popMats[0],x,y));
									}
									gsl_matrix_set(topol,i,j,top/bottom);
									gsl_matrix_int_set(moveType,i,j,0);
									dim1[nonZeroCount] = i;
									dim2[nonZeroCount] = j;
									nonZeroCount += 1;
								}
							}
						}
						else{
							if(matrixSum(delta->popMats[1]) == 1 && countNegativeEntries(delta->popMats[0]) == 0 && S->states[j]->aCounts[1] >= 1 ){
								entriesGreaterThan(delta->popMats[1], activePos1, &nlin1,0);//old lineages
								entriesLessThan(delta->popMats[1], activePos2, &nlin2,0);//new lineages
								if(nlin2 != 0 && nlin2 < 2){
									compat=0;
									if(nlin1 == 1){
										if(2*activePos1[0][0] == activePos2[0][0] && 2*activePos1[0][1] == activePos2[0][1])
											compat = 1;
									}
									if(nlin1==2){
										if(activePos1[0][0] + activePos1[1][0] == activePos2[0][0] &&
											activePos1[0][1] + activePos1[1][1] == activePos2[0][1])
											compat = 1;
									}
									if(compat != 1){

										gsl_matrix_set(topol,i,j,0.0);
										gsl_matrix_int_set(moveType,i,j,666);
									}
									else{
										top = 1.0;

										bottom = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1) /2.0 ;
										for(k=0;k<nlin1;k++){
											x = activePos1[k][0];
											y = activePos1[k][1];

											top *= xchoosey(gsl_matrix_int_get(S->states[i]->popMats[1],x,y),
												gsl_matrix_int_get(delta->popMats[1],x,y));
										}

										gsl_matrix_set(topol,i,j,top/bottom);
										gsl_matrix_int_set(moveType,i,j,1);
										dim1[nonZeroCount] = i;
										dim2[nonZeroCount] = j;
										nonZeroCount += 1;
									}
								}
							}
						}
						break;
					}
					//afsObjectFree(delta);
				}
			}
		}

	}
	//cleanup and free
	for(i=0;i<size;i++){
		free(activePos1[i]);
		free(activePos2[i]);
	}
	free(activePos1);
	free(activePos2);
	*nnz = nonZeroCount;
}

void coalMarkovChainTopologyMatrix_sparse_pthread(afsStateSpace *S,double *topol, int *moveType, int *dim1, int *dim2, int start, int stop, int *nnz){
	int i,j,k,steps,x,y;
	double top,bottom;
	afsObject *delta;
	int **activePos1, **activePos2; //first dimension here relates to 2 popns
	int nlin1,nlin2,compat;
	int nonZeroCount, tot, tCount;
	int size = 400;
	nlin1 = nlin2 = 0;

	nonZeroCount = 0;
	*nnz = 0;
	//arrays for finding positions of entries
	activePos1 = malloc(size*sizeof(int *));
	activePos2 = malloc(size*sizeof(int *));
	for(i=0;i<size;i++){
		activePos1[i] = malloc(2*sizeof(int));
		activePos1[i][0] = activePos1[i][1] = 0;
		activePos2[i] = malloc(2*sizeof(int));
		activePos2[i][0] = activePos2[i][1] = 0;
	}

	tot = S->nstates * S->nstates;
	tCount = 0;
	delta = afsObjectNewFrom(S->states[0]);
	for(i=start;i<stop;i++){
		for(j=0;j<S->nstates;j++){
		//	printf("nonzerocount=%d\n",nonZeroCount);

			if(S->states[i]->nalleles == 1){
				//set correct absorbing conditions
				if(i==j){
					topol[nonZeroCount] = 1.0;
					moveType[nonZeroCount] = -1;
					dim1[nonZeroCount] = i;
					dim2[nonZeroCount] = j;
					nonZeroCount += 1;
				}
			}

			else{ //ith not root; what is the move?
				steps = S->states[i]->nalleles - S->states[j]->nalleles ;
				if(steps < 2){
					delta = afsObjectDelta(S->states[i],S->states[j]);
					switch(steps){
						case 0:
						// 0 "steps", so no changes in nalleles between two
						//check counts and positions
						if(abs(matrixSum(delta->popMats[0])) == 1){
							nonZeroEntries(delta->popMats[0], activePos1, &nlin1);
							nonZeroEntries(delta->popMats[1], activePos2, &nlin2);

							//migration?
							if(nlin1 == nlin2 && activePos1[0][0] == activePos2[0][0] &&
								(gsl_matrix_int_min(delta->popMats[0]) >= 0 || gsl_matrix_int_min(delta->popMats[1]) >= 0)){

								//migration popn1?
								if(matrixSum(delta->popMats[0]) == 1){
									topol[nonZeroCount] = (double) gsl_matrix_int_get(S->states[i]->popMats[0],
										activePos1[0][0], activePos1[0][1]) / S->states[i]->aCounts[0];
									moveType[nonZeroCount] = 2;
									dim1[nonZeroCount] = i;
									dim2[nonZeroCount] = j;
									nonZeroCount += 1;
								}
								else{ //migration popn2
									topol[nonZeroCount] = (double) gsl_matrix_int_get(S->states[i]->popMats[1],
										activePos2[0][0],activePos2[0][1]) / S->states[i]->aCounts[1];
									moveType[nonZeroCount] = 3;

									dim1[nonZeroCount] = i;
									dim2[nonZeroCount] = j;
									nonZeroCount += 1;
								}
							}
						}
						break;

						case 1://coalescent?
								//coal popn1?
						if(matrixSum(delta->popMats[0]) == 1 && countNegativeEntries(delta->popMats[1]) == 0 && S->states[j]->aCounts[0] >= 1 ){
							entriesGreaterThan(delta->popMats[0], activePos1, &nlin1,0);//old lineages
							entriesLessThan(delta->popMats[0], activePos2, &nlin2,0);//new lineages
							if(nlin2 != 0 && nlin2 < 2){

								compat=0;
								if(nlin1 == 1){
									if(2*activePos1[0][0] == activePos2[0][0] && 2*activePos1[0][1] == activePos2[0][1])
										compat = 1;
								}
								if(nlin1==2){
									if(activePos1[0][0] + activePos1[1][0] == activePos2[0][0] &&
										activePos1[0][1] + activePos1[1][1] == activePos2[0][1])
										compat = 1;
								}
								if(compat == 1){
									top = 1.0;
									bottom = S->states[i]->aCounts[0] * (S->states[i]->aCounts[0]-1) /2.0 ;
									for(k=0;k<nlin1;k++){
										x = activePos1[k][0];
										y = activePos1[k][1];
										top *= (float) xchoosey(gsl_matrix_int_get(S->states[i]->popMats[0],x,y),
											gsl_matrix_int_get(delta->popMats[0],x,y));
									}

									topol[nonZeroCount] = top/bottom;
									moveType[nonZeroCount] = 0;
									dim1[nonZeroCount] = i;
									dim2[nonZeroCount] = j;
									nonZeroCount += 1;
								}
							}
						}
						else{
							if(matrixSum(delta->popMats[1]) == 1 && countNegativeEntries(delta->popMats[0]) == 0 && S->states[j]->aCounts[1] >= 1 ){
								entriesGreaterThan(delta->popMats[1], activePos1, &nlin1,0);//old lineages
								entriesLessThan(delta->popMats[1], activePos2, &nlin2,0);//new lineages
								if(nlin2 != 0 && nlin2 < 2){
									compat=0;
									if(nlin1 == 1){
										if(2*activePos1[0][0] == activePos2[0][0] && 2*activePos1[0][1] == activePos2[0][1])
											compat = 1;
									}
									if(nlin1==2){
										if(activePos1[0][0] + activePos1[1][0] == activePos2[0][0] &&
											activePos1[0][1] + activePos1[1][1] == activePos2[0][1])
											compat = 1;
									}
									if(compat == 1){
										top = 1.0;
										bottom = S->states[i]->aCounts[1] * (S->states[i]->aCounts[1]-1) /2.0 ;
										for(k=0;k<nlin1;k++){
											x = activePos1[k][0];
											y = activePos1[k][1];

											top *= xchoosey(gsl_matrix_int_get(S->states[i]->popMats[1],x,y),
												gsl_matrix_int_get(delta->popMats[1],x,y));
										}

										topol[nonZeroCount] = top/bottom;
										moveType[nonZeroCount] = 1;
										dim1[nonZeroCount] = i;
										dim2[nonZeroCount] = j;
										nonZeroCount += 1;
									}
								}
							}
						}
						break;
					}
					//afsObjectFree(delta);
				}
			}
		}
	}
	//cleanup and free
	for(i=0;i<size;i++){
		free(activePos1[i]);
		free(activePos2[i]);
	}
	free(activePos1);
	free(activePos2);
	*nnz = nonZeroCount;
}

max_stars_repo_name: dortegadelv/IMaDNA

max_issues_repo_name: dortegadelv/IMaDNA


text:
#pragma once

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include <time.h>

void random_matrix(double *matrix, int N);
void print_vector(double *vector, int N);
void print_matrix(double *matrix, int N);
void matrix_vector_mult(double *matrix, double *vector,
                        double *result, int N);
void matrix_matrix_mult(double *mat_1, double *mat_2,
                        double *result, int N);
void matrix_power(double *matrix, double *result,
                  int N, int power);


text:
/** ***************************************************************************
*
*  Copyright 2018 Jose Fernando Lopez Fernandez - All Rights Reserved
*  See the file AUTHORS included with the application distribution for
*  specific information.
*
*  Redistribution and use in source and binary forms, with or without
*  modification, are permitted provided that the following
*  conditions are met:
*
*   * Redistributions of source code must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*
*   * Redistributions in binary form must reproduce the above
*     copyright notice, this list of conditions and the following
*     disclaimer in the documentation and/or other materials provided
*     with the distribution.
*
*     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
*     CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
*     INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
*     MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
*     DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
*     CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
*     SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
*     NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
*     LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
*     HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
*     CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
*     OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
*     EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*  ***************************************************************************
*
*  @author Jose Fernando Lopez Fernandez
*  @date 15-June-2018
*  @brief This is the main keygen driver file
*
*  Tasks:
*      1. @todo Add file details
*      2. @todo Add file details
*
*  **************************************************************************/

#ifndef KEYGEN_INCLUDE_KEYGEN_H_
#define KEYGEN_INCLUDE_KEYGEN_H_

#define TRUE  1
#define FALSE 0

/** @def _CRT_SECURE_NO_WARNINGS
*  @brief This preprocessor definition disables Visual Studio's warnings
*  when not using the safe string functions provided by the MSVC runtime.
*
*  Disable the safe string library due to its lack of compatibility with
*  *nix systems.
*
*/

#define _CRT_SECURE_NO_WARNINGS

#ifndef _CRT_SECURE_NO_WARNINGS
#error "The safe string library is not compatible with *nix systems"
#endif

/** Include Microsoft's Guideline Support Library (GSL).
*
*  The Guideline Support Library (GSL) contains functions and types that are
*  suggested for use by the C++ Core Guidelines maintained by the Standard C++
*  Foundation. This repo contains Microsoft's implementation of GSL.
*
*  The library includes types like span<T>, string_span, owner<> and others.
*
*  The entire implementation is provided inline in the headers under the gsl
*  directory. The implementation generally assumes a platform that implements
*  C++14 support. There are specific workarounds to support MSVC 2015.
*
*  While some types have been broken out into their own headers
*  (e.g. gsl/span), it is simplest to just include gsl/gsl and gain access to
*  the entire library.
*
*/

#include <gsl/gsl>

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cstddef>
#include <cstdint>
#include <cerrno>
#include <csignal>
#include <cmath>
#include <ctime>

/** C++ Standard Library Headers
*
*/

#include <algorithm>
#include <chrono>
#include <iostream>
#include <limits>
#include <ratio>
#include <sstream>
#include <string>
#include <string_view>
#include <vector>
#include <utility>

/** @file config.h
 *  @brief This file contains customization settings and project metadata such
 *  as the current release version, compilation settings, etc.
 *
 */

#include "config.h"
#include "password.h"

#endif // KEYGEN_KEYGEN_H_

max_stars_repo_name: lopezfjose/KeyGen

max_issues_repo_name: lopezfjose/KeyGen

max_forks_repo_name: lopezfjose/KeyGen


text:
/**
 *
 * @file qwrapper_zlantr.c
 *
 *  PLASMA core_blas quark wrapper
 *  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
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

void
CORE_zlantr_quark(Quark *quark);
void
CORE_zlantr_f1_quark(Quark *quark);

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlantr(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                       const PLASMA_Complex64_t *A, int LDA, int szeA,
                       int szeW, double *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;
    QUARK_Insert_Task(quark, CORE_zlantr_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(PLASMA_enum),                &diag,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(double)*szeW,                 NULL,          SCRATCH,
        sizeof(double),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlantr_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                          const PLASMA_Complex64_t *A, int LDA, int szeA,
                          int szeW, double *result,
                          double *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_zlantr_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_zlantr_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT,
            sizeof(double)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlantr_quark = PCORE_zlantr_quark
#define CORE_zlantr_quark PCORE_zlantr_quark
#endif
void CORE_zlantr_quark(Quark *quark)
{
    double *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    double *work;

    quark_unpack_args_9(quark, norm, uplo, diag, M, N, A, LDA, work, normA);
    *normA = LAPACKE_zlantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlantr_f1_quark = PCORE_zlantr_f1_quark
#define CORE_zlantr_f1_quark PCORE_zlantr_f1_quark
#endif
void CORE_zlantr_f1_quark(Quark *quark)
{
    double *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    double *work;
    double *fake;

    quark_unpack_args_10(quark, norm, uplo, diag, M, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_zlantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}


text:
/*
Copyright (C) 2019-2020 JingWeiZhangHuai <jingweizhanghuai@163.com>
Licensed under the Apache License, Version 2.0; 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.
*/

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

#include <cblas.h>
#include "morn_tensor.h"

struct TensorBatchNormPara
{
    MLayer *prev;

    int res_valid;

    float rate;
    float decay;
    float momentum;
};

void *mTensorBatchNormPara(MList *ini,char *name)
{
    struct TensorBatchNormPara *para = (struct TensorBatchNormPara *)mMalloc(sizeof(struct TensorBatchNormPara));

    char *value = mINIRead(ini,name,"prev");
    para->prev = mNetworkLayer(ini,value);
    mException((para->prev == NULL),EXIT,"invalid prev");

    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    value = mINIRead(ini,name,"rate");
    if(value != NULL) para->rate = atof(value);
    else
    {
        value = mINIRead(ini,"para","rate");
        if(value != NULL) para->rate = atof(value);
        else              para->rate = 0.001;
    }

    value = mINIRead(ini,name,"decay");
    if(value != NULL) para->decay = atof(value);
    else
    {
        value = mINIRead(ini,"para","decay");
        if(value != NULL) para->decay = atof(value);
        else              para->decay = 0.01;
    }
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    value = mINIRead(ini,name,"momentum");
    if(value != NULL) para->momentum = atof(value);
    else
    {
        value = mINIRead(ini,"para","momentum");
        if(value != NULL) para->momentum = atof(value);
        else              para->momentum = 0.9;
    }

    return para;
}

struct HandleTensorBatchNorm
{
    float *k;
    float *b;

    float *k_update;
    float *b_update;

    double *mean;
    double *var;
    float *roll_mean;
    float *roll_var;
};
void endTensorBatchNorm(void *info)
{
    struct HandleTensorBatchNorm *handle = (struct HandleTensorBatchNorm *)info;
    if(handle->k        != NULL) mFree(handle->k);
    if(handle->b        != NULL) mFree(handle->b);
    if(handle->mean     != NULL) mFree(handle->mean);
    if(handle->var      != NULL) mFree(handle->var );
    if(handle->roll_mean!= NULL) mFree(handle->roll_mean);
    if(handle->roll_var != NULL) mFree(handle->roll_var );
}
#define HASH_TensorBatchNorm 0xde8952f6

void TensorBatchNormSet(MLayer *layer)
{
    if(layer->state != DFLT) return;
    struct TensorBatchNormPara *para = (struct TensorBatchNormPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->tns;

    MHandle *hdl=mHandle(out,TensorBatchNorm);
    struct HandleTensorBatchNorm *handle = (struct HandleTensorBatchNorm *)(hdl->handle);

    if(handle->k != NULL) {mFree(handle->k);} handle->k =(float *)mMalloc(in->channel*sizeof(float));
    if(handle->b != NULL) {mFree(handle->b);} handle->b =(float *)mMalloc(in->channel*sizeof(float));

    if(handle->mean!= NULL) {mFree(handle->mean);} handle->mean=(double *)mMalloc(in->channel*sizeof(double));
    if(handle->var != NULL) {mFree(handle->var );} handle->var =(double *)mMalloc(in->channel*sizeof(double));

    if(handle->roll_mean!= NULL) {mFree(handle->roll_mean);} handle->roll_mean=(float *)mMalloc(in->channel*sizeof(float));
    if(handle->roll_var != NULL) {mFree(handle->roll_var );} handle->roll_var =(float *)mMalloc(in->channel*sizeof(float));

    mTensorRedefine(out,in->batch,in->channel,in->height,in->width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,out->data);
        else                    mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);

        handle->k_update = handle->roll_mean;
        handle->b_update = handle->roll_var ;
        memset(handle->k_update,0,in->channel*sizeof(float));
        memset(handle->b_update,0,in->channel*sizeof(float));
    }

    if(morn_network_parafile==NULL)
    {
        for(int i=0;i<in->channel;i++){handle->k[i] = 1.0f;handle->b[i] = 1.0f;}
        // printf("handle->k[0] is %f\n",handle->k[0]);
        // printf("handle->b[0] is %f\n",handle->b[0]);
    }
    else
    {
        mNetworkParaRead(layer,"scale",handle->k,in->channel*sizeof(float));
        mNetworkParaRead(layer,"bias" ,handle->b,in->channel*sizeof(float));
        // printf("handle->k[0] is %f\n",handle->k[0]);
        // printf("handle->b[0] is %f\n",handle->b[0]);
    }
    hdl->valid = 1;
}

void mTensorBatchNormForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("BatchNorm",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorBatchNormPara *para = (struct TensorBatchNormPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;

    TensorBatchNormSet(layer);

    MHandle *hdl=mHandle(out,TensorBatchNorm);
    struct HandleTensorBatchNorm *handle = (struct HandleTensorBatchNorm *)(hdl->handle);

    int in_size = in->height*in->width;

    memset(handle->mean,0,in->channel*sizeof(double));
    memset(handle->var ,0,in->channel*sizeof(double));
    for(int b=0;b<in->batch;b++)for(int c=0;c<in->channel;c++)for(int i=0;i<in_size;i++)
    {
        double data = (double)(in->data[b][c*in_size+i]);
        handle->mean[c] += data;
        handle->var [c] += data*data;
    }
    for(int c=0;c<in->channel;c++)
    {
        handle->mean[c] = handle->mean[c]/(double)(in->batch*in_size);
        handle->var [c] = handle->var [c]/(double)(in->batch*in_size);
        handle->var [c]-= (handle->mean[c]*handle->mean[c]);
        handle->var [c] = (handle->var[c]<=0.0f)?0.00001f:sqrt(handle->var[c]);
        // printf("handle->mean[c] is %f,handle->var[c] is %f\n",handle->mean[c],handle->var[c]);
    }

    if(morn_network_flag == MORN_PREDICT)
    {
        for(int c=0;c<in->channel;c++)
        {
            if(layer->state == DFLT)
            {
                handle->roll_mean[c] = handle->mean[c];
                handle->roll_var [c] = handle->var [c];
            }

            handle->mean[c] = handle->roll_mean[c]*0.99f+handle->mean[c]*0.01;
            handle->var [c] = handle->roll_var [c]*0.99f+handle->var [c]*0.01;
            handle->roll_mean[c] = handle->mean[c];
            handle->roll_var [c] = handle->var [c];
        }
    }

    for(int b=0;b<in->batch;b++)for(int c=0;c<in->channel;c++)for(int i=0;i<in_size;i++)
    {
        float data = in->data[b][c*in_size+i];
        out->data[b][c*in_size+i] = ((data-handle->mean[c])/handle->var[c])*handle->k[c]+handle->b[c];
    }

    layer->state = MORN_FORWARD;
}

void mTensorBatchNormBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("BatchNorm",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorBatchNormPara *para = (struct TensorBatchNormPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out= layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorBatchNorm);
    struct HandleTensorBatchNorm *handle = (struct HandleTensorBatchNorm *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    mNetworkParaWrite(layer,"scale",handle->k,in->channel*sizeof(float));
    mNetworkParaWrite(layer,"bias" ,handle->b,in->channel*sizeof(float));

    int in_size = in->height*in->width;

    for(int c=0;c<in->channel;c++)
    {
        handle->k_update[c] *= para->momentum;
        handle->b_update[c] *= para->momentum;
    }
    for(int b=0;b<in->batch;b++)for(int c=0;c<in->channel;c++)for(int i=0;i<in_size;i++)
    {
        int idx = c*in_size+i;
        handle->k_update[c] += out->data[b][idx]*((in->data[b][idx]-handle->mean[c])/handle->var[c]);
        handle->b_update[c] += out->data[b][idx];
    }
    for(int c=0;c<in->channel;c++)
    {
        handle->k[c] = (handle->k[c]*(1.0f-(para->decay*para->rate)))-(handle->k_update[c]*para->rate/(float)(in->batch*in_size));
        handle->b[c] = (handle->b[c]*(1.0f-(para->decay*para->rate)))-(handle->b_update[c]*para->rate/(float)(in->batch*in_size));
    }

    if(para->res_valid==0) return;

    for(int b=0;b<in->batch;b++)for(int c=0;c<in->channel;c++)for(int i=0;i<in_size;i++)
    {
        int idx = c*in_size+i;
        float data;
        data = (in->data[b][idx] - handle->mean[c])/handle->var[c];
        data = (out->data[b][idx]*handle->k[c]/handle->var[c])*(1.0+(1.0f-data*data)/((float)(in->batch*in_size)));

        res->data[b][idx]=(para->prev->state == MORN_FORWARD)?data:(res->data[b][idx]+data);
    }

    para->prev->state = MORN_BACKWARD;
}

max_stars_repo_name: ishine/Morn

max_issues_repo_name: ishine/Morn

max_forks_repo_name: ishine/Morn


text:
#ifndef STATEMENT_INCLUDE
#define STATEMENT_INCLUDE

#include <string>

#include <gsl/gsl>

#include "executionContent.h"
#include "yaml.h"

namespace execHelper {
namespace test {
namespace baseUtils {
using StatementKey = std::string;
using StatementCollection = ExecutionContent::ConfigCommand;

class Statement {
  public:
    virtual ~Statement() = default;

    virtual StatementKey getKey() const noexcept = 0;
    virtual void write(gsl::not_null<YamlWriter*> yaml,
                       const std::string& command) const noexcept = 0;

    virtual inline unsigned int getNumberOfExecutions() const noexcept {
        return m_execution.getNumberOfExecutions();
    }

    virtual inline void resetExecutions() noexcept { m_execution.clear(); }

  protected:
    Statement(ReturnCode returnCode) noexcept : m_execution(returnCode) {}

    ExecutionContent m_execution;
};
} // namespace baseUtils
} // namespace test
} // namespace execHelper

#endif /* STATEMENT_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
#ifndef __GSL_MANAGED_H__
#define __GSL_MANAGED_H__

#include "cuda.h"
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>

gsl_matrix *gsl_matrix_alloc_managed(const size_t n1, const size_t n2);
gsl_matrix *gsl_matrix_calloc_managed(const size_t n1, const size_t n2);
void gsl_matrix_free_managed(gsl_matrix *m);
__host__ __device__ double *gsl_matrix_ptr_managed(gsl_matrix * m, const size_t i, const size_t j);

gsl_vector *gsl_vector_alloc_managed(const size_t n, bool readmostly);
gsl_vector *gsl_vector_calloc_managed(const size_t n, bool readmostly);
__host__ __device__ double gsl_vector_get_managed(gsl_vector *vec, const size_t i);
void gsl_vector_free_managed(gsl_vector *m);

#endif /* __GSL_MANAGED_H__ */

max_stars_repo_name: artis-mcrt/artis

max_issues_repo_name: artis-mcrt/artis

max_forks_repo_name: artis-mcrt/artis


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_hemm (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int side = 141;
   int uplo = 121;
   int M = 1;
   int N = 2;
   float alpha[2] = {0.0f, 0.1f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.126f, 0.079f };
   int lda = 1;
   float B[] = { -0.954f, -0.059f, 0.296f, -0.988f };
   int ldb = 2;
   float C[] = { -0.859f, -0.731f, 0.737f, 0.593f };
   int ldc = 2;
   float C_expected[] = { 0.0723566f, -0.0738796f, -0.0717488f, 0.0699704f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1550) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1550) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 141;
   int uplo = 121;
   int M = 1;
   int N = 2;
   float alpha[2] = {0.0f, 0.1f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { 0.652f, 0.584f };
   int lda = 1;
   float B[] = { -0.983f, -0.734f, -0.422f, -0.825f };
   int ldb = 1;
   float C[] = { 0.387f, 0.341f, -0.734f, 0.632f };
   int ldc = 1;
   float C_expected[] = { 0.0137568f, -0.0253916f, -0.00941f, -0.100914f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1551) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1551) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 142;
   int uplo = 121;
   int M = 1;
   int N = 2;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.78f, 0.885f, 0.507f, 0.765f, 0.911f, -0.461f, 0.707f, 0.508f };
   int lda = 2;
   float B[] = { -0.905f, 0.633f, 0.85f, -0.943f };
   int ldb = 2;
   float C[] = { 0.045f, -0.237f, 0.078f, -0.252f };
   int ldc = 2;
   float C_expected[] = { 0.589611f, -0.759345f, 0.960095f, -0.09013f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1552) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1552) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 142;
   int uplo = 121;
   int M = 1;
   int N = 2;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.947f, 0.939f, -0.267f, -0.819f, -0.827f, -0.937f, 0.991f, 0.838f };
   int lda = 2;
   float B[] = { 0.871f, -0.988f, -0.232f, -0.434f };
   int ldb = 1;
   float C[] = { -0.261f, 0.927f, -0.351f, -0.203f };
   int ldc = 1;
   float C_expected[] = { 1.0551f, 0.496359f, 0.780145f, -1.67298f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1553) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1553) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 141;
   int uplo = 122;
   int M = 1;
   int N = 2;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {0.0f, 0.0f};
   float A[] = { -0.593f, -0.9f };
   int lda = 1;
   float B[] = { -0.861f, 0.747f, -0.984f, 0.595f };
   int ldb = 2;
   float C[] = { -0.589f, -0.671f, -0.011f, -0.417f };
   int ldc = 2;
   float C_expected[] = { -0.510573f, 0.442971f, -0.583512f, 0.352835f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1554) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1554) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 141;
   int uplo = 122;
   int M = 1;
   int N = 2;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {0.0f, 0.0f};
   float A[] = { -0.79f, 0.132f };
   int lda = 1;
   float B[] = { -0.243f, -0.12f, 0.633f, -0.556f };
   int ldb = 1;
   float C[] = { -0.658f, -0.74f, -0.47f, 0.481f };
   int ldc = 1;
   float C_expected[] = { -0.19197f, -0.0948f, 0.50007f, -0.43924f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1555) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1555) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 142;
   int uplo = 122;
   int M = 1;
   int N = 2;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { -0.114f, -0.515f, -0.513f, -0.527f, -0.995f, 0.986f, 0.229f, -0.076f };
   int lda = 2;
   float B[] = { 0.084f, 0.522f, 0.61f, 0.694f };
   int ldb = 2;
   float C[] = { 0.802f, 0.136f, -0.161f, -0.364f };
   int ldc = 2;
   float C_expected[] = { 0.269101f, 0.716492f, 0.237088f, 0.0290666f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1556) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1556) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 142;
   int uplo = 122;
   int M = 1;
   int N = 2;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { 0.798f, -0.324f, -0.693f, -0.893f, -0.223f, 0.749f, 0.102f, -0.357f };
   int lda = 2;
   float B[] = { -0.572f, -0.569f, -0.391f, -0.938f };
   int ldb = 1;
   float C[] = { 0.152f, -0.834f, -0.633f, -0.473f };
   int ldc = 1;
   float C_expected[] = { 1.08642f, -0.113853f, 0.234826f, -0.48289f };
   cblas_chemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "chemm(case 1557) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "chemm(case 1557) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 141;
   int uplo = 121;
   int M = 1;
   int N = 2;
   double alpha[2] = {0, 0.1};
   double beta[2] = {0, 0.1};
   double A[] = { -0.359, 0.089 };
   int lda = 1;
   double B[] = { -0.451, -0.337, -0.901, -0.871 };
   int ldb = 2;
   double C[] = { 0.729, 0.631, 0.364, 0.246 };
   int ldc = 2;
   double C_expected[] = { -0.0751983, 0.0890909, -0.0558689, 0.0687459 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1558) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1558) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 141;
   int uplo = 121;
   int M = 1;
   int N = 2;
   double alpha[2] = {0, 0.1};
   double beta[2] = {0, 0.1};
   double A[] = { 0.044, -0.496 };
   int lda = 1;
   double B[] = { -0.674, 0.281, 0.366, 0.888 };
   int ldb = 1;
   double C[] = { -0.9, 0.919, 0.857, -0.049 };
   int ldc = 1;
   double C_expected[] = { -0.0931364, -0.0929656, 0.0009928, 0.0873104 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1559) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1559) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 142;
   int uplo = 121;
   int M = 1;
   int N = 2;
   double alpha[2] = {0, 0};
   double beta[2] = {0, 0.1};
   double A[] = { -0.314, 0.115, 0.114, 0.878, 0.961, -0.224, 0.973, 0.771 };
   int lda = 2;
   double B[] = { 0.5, -0.016, -0.5, 0.149 };
   int ldb = 2;
   double C[] = { -0.054, 0.064, 0.02, 0.245 };
   int ldc = 2;
   double C_expected[] = { -0.0064, -0.0054, -0.0245, 0.002 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1560) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1560) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 142;
   int uplo = 121;
   int M = 1;
   int N = 2;
   double alpha[2] = {0, 0};
   double beta[2] = {0, 0.1};
   double A[] = { 0.186, 0.578, 0.797, -0.957, -0.539, -0.969, -0.21, 0.354 };
   int lda = 2;
   double B[] = { 0.641, -0.968, 0.15, -0.569 };
   int ldb = 1;
   double C[] = { -0.556, -0.9, 0.197, 0.31 };
   int ldc = 1;
   double C_expected[] = { 0.09, -0.0556, -0.031, 0.0197 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1561) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1561) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 141;
   int uplo = 122;
   int M = 1;
   int N = 2;
   double alpha[2] = {1, 0};
   double beta[2] = {1, 0};
   double A[] = { 0.323, 0.641 };
   int lda = 1;
   double B[] = { -0.188, 0.091, -0.235, 0.523 };
   int ldb = 2;
   double C[] = { 0.919, 0.806, 0.823, -0.94 };
   int ldc = 2;
   double C_expected[] = { 0.858276, 0.835393, 0.747095, -0.771071 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1562) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1562) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 141;
   int uplo = 122;
   int M = 1;
   int N = 2;
   double alpha[2] = {1, 0};
   double beta[2] = {1, 0};
   double A[] = { -0.688, 0.915 };
   int lda = 1;
   double B[] = { 0.914, -0.204, 0.205, -0.476 };
   int ldb = 1;
   double C[] = { 0.27, -0.628, -0.079, 0.507 };
   int ldc = 1;
   double C_expected[] = { -0.358832, -0.487648, -0.22004, 0.834488 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1563) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1563) imag");
     };
   };
  };

  {
   int order = 101;
   int side = 142;
   int uplo = 122;
   int M = 1;
   int N = 2;
   double alpha[2] = {0, 1};
   double beta[2] = {0, 0.1};
   double A[] = { 0.681, 0.574, -0.425, -0.64, 0.792, 0.661, -0.009, 0.005 };
   int lda = 2;
   double B[] = { -0.221, 0.554, -0.465, -0.95 };
   int ldb = 2;
   double C[] = { 0.331, -0.958, -0.826, -0.972 };
   int ldc = 2;
   double C_expected[] = { 0.778291, 0.142269, -0.496199, 0.112747 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1564) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1564) imag");
     };
   };
  };

  {
   int order = 102;
   int side = 142;
   int uplo = 122;
   int M = 1;
   int N = 2;
   double alpha[2] = {0, 1};
   double beta[2] = {0, 0.1};
   double A[] = { 0.959, 0.34, -0.23, 0.064, 0.516, -0.275, 0.714, 0.899 };
   int lda = 2;
   double B[] = { -0.502, -0.987, -0.134, 0.215 };
   int ldb = 1;
   double C[] = { 0.929, 0.181, -0.16, -0.921 };
   int ldc = 1;
   double C_expected[] = { 0.986459, -0.371458, -0.320548, -0.059384 };
   cblas_zhemm(order, side, uplo, M, N, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zhemm(case 1565) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zhemm(case 1565) imag");
     };
   };
  };

}


text:
/* specfunc/bessel_K0.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * Copyright (C) 2016 Pavel Holoborodko, 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

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

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/*
 Minimax rational approximation for [0,1), peak relative error = 2.04*GSL_DBL_EPSILON.
 Source: http://www.advanpix.com/?p=3812
*/
static double k0_poly[8] = {
   1.1593151565841244842077226e-01,
   2.7898287891460317300886539e-01,
   2.5248929932161220559969776e-02,
   8.4603509072136578707676406e-04,
   1.4914719243067801775856150e-05,
   1.6271068931224552553548933e-07,
   1.2082660336282566759313543e-09,
   6.6117104672254184399933971e-12
};

static double i0_poly[7] = {
   1.0000000000000000044974165e+00,
   2.4999999999999822316775454e-01,
   2.7777777777892149148858521e-02,
   1.7361111083544590676709592e-03,
   6.9444476047072424198677755e-05,
   1.9288265756466775034067979e-06,
   3.9908220583262192851839992e-08
};

/*
 Chebyshev expansion for [1,8], peak relative error = 1.28*GSL_DBL_EPSILON.
 Source: Pavel Holoborodko.
*/
static double ak0_data[24] = {
  -3.28737867094650101e-02,
  -4.49369057710236880e-02,
  +2.98149992004308095e-03,
  -3.03693649396187920e-04,
  +3.91085569307646836e-05,
  -5.86872422399215952e-06,
  +9.82873709937322009e-07,
  -1.78978645055651171e-07,
  +3.48332306845240957e-08,
  -7.15909210462546599e-09,
  +1.54019930048919494e-09,
  -3.44555485579194210e-10,
  +7.97356101783753023e-11,
  -1.90090968913069735e-11,
  +4.65295609304114621e-12,
  -1.16614287433470780e-12,
  +2.98554375218596891e-13,
  -7.79276979512292169e-14,
  +2.07027467168948402e-14,
  -5.58987860393825313e-15,
  +1.53202965950646914e-15,
  -4.25737536712188186e-16,
  +1.19840238501357389e-16,
  -3.41407346762502397e-17
};

static cheb_series ak0_cs = {
  ak0_data,
  23,
  -1, 1,
  10
};

/*
 Chebyshev expansion for [8,inf), peak relative error = 1.25*GSL_DBL_EPSILON.
 Source: SLATEC/dbsk0e.f
*/
static double ak02_data[14] = {
  -.1201869826307592240E-1,
  -.9174852691025695311E-2,
  +.1444550931775005821E-3,
  -.4013614175435709729E-5,
  +.1567831810852310673E-6,
  -.7770110438521737710E-8,
  +.4611182576179717883E-9,
  -.3158592997860565771E-10,
  +.2435018039365041128E-11,
  -.2074331387398347898E-12,
  +.1925787280589917085E-13,
  -.1927554805838956104E-14,
  +.2062198029197818278E-15,
  -.2341685117579242403E-16
};

static cheb_series ak02_cs = {
  ak02_data,
  13,
  -1, 1,
  8
};

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

int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 1.0) {
    const double lx = log(x);
    const double ex = exp(x);
    const double x2 = x*x;
    result->val  = ex * (gsl_poly_eval(k0_poly,8,x2)-lx*(1.0+0.25*x2*gsl_poly_eval(i0_poly,7,0.25*x2)));
    result->err  = ex * (1.6+fabs(lx)*0.6) * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= 8.0) {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak0_cs, (16.0/x-9.0)/7.0, &c);
    result->val  = (1.203125 + c.val) / sx; /* 1.203125 = 77/64 */
    result->err  = c.err / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak02_cs, 16.0/x-1.0, &c);
    result->val  = (1.25 + c.val) / sx;
    result->err  = (c.err + GSL_DBL_EPSILON) / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 1.0) {
    const double lx = log(x);
    const double x2 = x*x;
    result->val  = gsl_poly_eval(k0_poly,8,x2)-lx*(1.0+0.25*x2*gsl_poly_eval(i0_poly,7,0.25*x2));
    result->err  = (1.6+fabs(lx)*0.6) * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result K0_scaled;
    int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled);
    int stat_e  = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x),
                                           K0_scaled.val, K0_scaled.err,
                                           result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K0);
  }
}

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

#include "eval.h"

double gsl_sf_bessel_K0_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K0_scaled_e(x, &result));
}

double gsl_sf_bessel_K0(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K0_e(x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game