text:
/*
** Copyright (C) 2004 Jonathan G. Underwood <j.underwood@open.ac.uk>
**
** 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_SF_WIGNER_H__
#define  __GSL_SF_WIGNER_H__ 1
#include <gsl/gsl_sf_result.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
  double gsl_sf_wigner_3j (const int two_j1, const int two_j2,
                           const int two_j3, const int two_m1,
                           const int two_m2, const int two_m3);
int gsl_sf_wigner_3j_e (const int two_j1, const int two_j2, const int two_j3,
                        const int two_m1, const int two_m2, const int two_m3,
                        gsl_sf_result * result);

double gsl_sf_wigner_6j (const int two_j1, const int two_j2, const int two_j3,
                         const int two_j4, const int two_j5,
                         const int two_j6);
int gsl_sf_wigner_6j_e (const int two_j1, const int two_j2, const int two_j3,
                        const int two_j4, const int two_j5, const int two_j6,
                        gsl_sf_result * result);

double gsl_sf_wigner_9j (const int two_j1, const int two_j2, const int two_j3,
                         const int two_j4, const int two_j5, const int two_j6,
                         const int two_j7, const int two_j8,
                         const int two_j9);
int gsl_sf_wigner_9j_e (const int two_j1, const int two_j2, const int two_j3,
                        const int two_j4, const int two_j5, const int two_j6,
                        const int two_j7, const int two_j8, const int two_j9,
                        gsl_sf_result * result);

double gsl_sf_wigner_drot (const int two_j, const int two_m1,
                           const int two_m2, const double theta);
int gsl_sf_wigner_drot_e (const int two_j, const int two_m1, const int two_m2,
                          const double theta, gsl_sf_result * result);

__END_DECLS
#endif /* __GSL_SF_WIGNER_H__ */

ext: h

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
#pragma once

#include <optional>

#include <gsl/span>

#include "Ds4InputData.h"

/**
 * \brief Serialized input report from a \c Ds4Device
 * \sa Ds4Device
 */
class Ds4Input
{
public:
	Ds4Input() = default;

	/**
	 * \brief Indicates if button press state has changed since the last poll.
	 */
	bool buttonsChanged = false; // TODO: private set

	/**
	 * \brief Indicates if touch state has changed since last poll.
	 */
	bool touchChanged = false; // TODO: private set

	/**
	 * \brief Buttons currently held.
	 * \sa Ds4Buttons, Ds4Buttons_t
	 */
	Ds4Buttons_t heldButtons = 0; // TODO: private set

	/**
	 * \brief Buttons pressed since last poll.
	 * \sa Ds4Buttons, Ds4Buttons_t
	 */
	Ds4Buttons_t pressedButtons = 0; // TODO: private set

	/**
	 * \brief Buttons released since last poll.
	 * \sa Ds4Buttons, Ds4Buttons_t
	 */
	Ds4Buttons_t releasedButtons = 0; // TODO: private set

	/**
	 * \brief Each axis that has changed since the last poll.
	 * \sa Ds4Axes, Ds4Axes_t
	 */
	Ds4Axes_t axes = 0; // TODO: private set

	Ds4InputData data {};

	/**
	 * \brief Updates serialized data using the given buffer.
	 * \param buffer Buffer containing raw input report data.
	 */
	void update(const gsl::span<uint8_t>& buffer);

	/**
	 * \brief Updates button change states since last poll.
	 */
	void updateChangedState();

	/**
	 * \brief Get the magnitude of an axis.
	 * \param axis The axis to retrieve.
	 * \param polarity The desired polarity of the axis, or \c std::nullopt for both positive and negative.
	 * \return The magnitude of the axis. If it does not align with the desired \a polarity, \c 0.0f is returned.
	 */
	[[nodiscard]] float getAxis(Ds4Axes_t axis, const std::optional<AxisPolarity>& polarity) const;

private:
	Ds4Buttons_t lastHeldButtons = 0;
	uint8_t lastTouchFrame {};

	void addButton(bool pressed, Ds4Buttons_t buttons);
	void updateButtons();
	void updateAxes(const Ds4InputData& last);
};


text:
/*! \file Exceptions.h
 *  \brief header file that defines exceptions
 */
// Function.h

#ifndef NBODYEXCEPTION_H
#define NBODYEXCEPTION_H

#include <exception>
#include <string>
#include <utility>

#include <gsl/gsl_errno.h>

namespace Math
{

class gsl_error : public std::exception {

private:
	int gsl_errno;
	std::string errmsg;

public:
	gsl_error(int gsl_errno) :
		gsl_errno(gsl_errno),
		errmsg(gsl_strerror(gsl_errno))
	{
	}

	int get_gsl_errno() {
		return gsl_errno;
	}

	const char *what() const noexcept
	{
		return errmsg.c_str();
	}

};

template<typename Func, typename ... Args>
void gsl_invoke(Func &&f, Args ... args)
{
	int status = f(std::forward<Args>(args)...);
	if (status != GSL_SUCCESS) {
		throw gsl_error(status);
	}
}

}

#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void add_MCMC_DS(MCMC_info *MCMC, const char *name, int n_D, int *D, double *DS, double *dDS, void *params, int n_arrays, ...) {
    int           i_array;
    int           i_D;
    int           n_M;
    MCMC_DS_info *new_DS;
    va_list       vargs;
    va_start(vargs, n_arrays);

    SID_log("Adding data set {%s} to MCMC structure...", SID_LOG_OPEN, name);

#if USE_CFITSIO == 0
    if(n_D > 1)
        SID_exit_error("You can only use 2D datasets if you compile with USE_CFITSIO=1 at the moment.", SID_ERROR_LOGIC);
#endif

    // Initialize new dataset
    init_MCMC_DS(&new_DS, name, n_D, D, DS, dDS, params, n_arrays, vargs);

    // Add new dataset to linked list
    if(MCMC->DS == NULL)
        MCMC->DS = new_DS;
    else
        MCMC->last->next = new_DS;
    MCMC->last = new_DS;

    MCMC->n_DS++;

    SID_log("Done.", SID_LOG_CLOSE);

    va_end(vargs);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
//===--- Sudoku/Board.h                                                 ---===//
//
// Data container for sudoku puzzles
//===----------------------------------------------------------------------===//
//
// accepts only squared dimensions
// gives full, row, column and block access
// no processing functionality
//
//===----------------------------------------------------------------------===//
#pragma once

#include "Board_Iterators.h"
#include "Board_Section.h"
#include "Location.h"
#include "Location_Utilities.h"
#include "Size.h"
#include "exceptions.h"

#include <gsl/gsl>

#include <array>
#include <initializer_list>

#include <algorithm>

#include "Board.fwd.h" // Forward declarations

#include <cassert>

namespace Sudoku
{
template<typename T, int N = 3>
class Board
{
	static_assert(N > 1, "Board.h base_size value too small");
	static constexpr auto Size = size_t{full_size<N>};

	using index       = gsl::index;
	using Location    = ::Sudoku::Location<N>;
	using Row         = Board_Section::Row<T, N>;
	using const_Row   = Board_Section::const_Row<T, N>;
	using Col         = Board_Section::Col<T, N>;
	using const_Col   = Board_Section::const_Col<T, N>;
	using Block       = Board_Section::Block<T, N>;
	using const_Block = Board_Section::const_Block<T, N>;

public:
	using value_type             = T;
	using size_type              = size_t;
	using difference_type        = std::ptrdiff_t;
	using reference              = value_type&;
	using const_reference        = value_type const&;
	using pointer                = value_type*;
	using const_pointer          = value_type const*;
	using iterator               = Board_iterator<T, N>;
	using const_iterator         = const_Board_iterator<T, N>;
	using reverse_iterator       = reverse_Board_iterator<T, N>;
	using const_reverse_iterator = const_reverse_Board_iterator<T, N>;

	constexpr Board() noexcept;
	explicit constexpr Board(const T& default_value);
	// NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
	constexpr Board(std::array<T, Size> const& list) // NOLINT(runtime/explicit)
		: board_(list)
	{
	}
	Board(std::initializer_list<T>); // construct from initializer_list

	void clear();

	// Query properties
	[[nodiscard]] static constexpr size_t size() noexcept { return Size; }
	[[nodiscard]] static constexpr size_t max_size() noexcept { return size(); }
	[[nodiscard]] static constexpr bool empty() noexcept { return false; }
	[[nodiscard]] bool operator==(const Board&) const;

	// Element access
	[[nodiscard]] constexpr T& front() noexcept { return board_.front(); }
	[[nodiscard]] constexpr T const& front() const noexcept
	{
		return board_.front();
	}
	[[nodiscard]] constexpr T& back() noexcept { return board_.back(); }
	[[nodiscard]] constexpr T const& back() const noexcept
	{
		return board_.back();
	}
	// Checked
	constexpr T& at(Location);
	[[nodiscard]] constexpr T const& at(Location) const;
	constexpr T& at(index row, index col);
	[[nodiscard]] constexpr T const& at(index row, index col) const;
	[[deprecated]] T& at(index elem);
	[[deprecated, nodiscard]] const T& at(index elem) const;
	// Unchecked
	constexpr T& operator[](Location) noexcept;
	constexpr T const& operator[](Location) const noexcept;

	// Element Selection Operator (using a proxy object)
	//   usable as [row][col] where col is processed by the (const_)Row
	Row operator[](index row_id) noexcept { return row(row_id); }
	constexpr const_Row operator[](index row_id) const noexcept
	{
		return row(row_id);
	}

	// Iterators
	constexpr iterator begin() noexcept;
	constexpr iterator end() noexcept;
	[[nodiscard]] constexpr const_iterator cbegin() const noexcept;
	[[nodiscard]] constexpr const_iterator cend() const noexcept;
	[[nodiscard]] constexpr const_iterator begin() const noexcept
	{
		return cbegin();
	}
	[[nodiscard]] constexpr const_iterator end() const noexcept
	{
		return cend();
	}
	constexpr reverse_iterator rbegin() noexcept;
	constexpr reverse_iterator rend() noexcept;
	[[nodiscard]] constexpr const_reverse_iterator crbegin() const noexcept;
	[[nodiscard]] constexpr const_reverse_iterator crend() const noexcept;
	[[nodiscard]] constexpr const_reverse_iterator rbegin() const noexcept
	{
		return crbegin();
	}
	[[nodiscard]] constexpr const_reverse_iterator rend() const noexcept
	{
		return crend();
	}

	// Sections
	// clang-format off
	Row row(index id) noexcept         { return Row(*this, id); }
	[[nodiscard]] constexpr const_Row   row(index id) const noexcept
									   { return const_Row(*this, id); }
	Row row(Location loc) noexcept     { return Row(*this, loc); }
	[[nodiscard]] constexpr const_Row   row(Location loc) const noexcept
									   { return const_Row(*this, loc); }
	Col col(index id) noexcept         { return Col(*this, id); }
	[[nodiscard]] constexpr const_Col   col(index id) const noexcept
									   { return const_Col(*this, id); }
	Col col(Location loc) noexcept     { return Col(*this, loc); }
	[[nodiscard]] constexpr const_Col   col(Location loc) const noexcept
									   { return const_Col(*this, loc); }
	Block block(index id) noexcept     { return Block(*this, id); }
	[[nodiscard]] constexpr const_Block block(index id) const noexcept
									   { return const_Block(*this, id); }
	Block block(Location loc) noexcept { return Block(*this, loc); }
	[[nodiscard]] constexpr const_Block block(Location loc) const noexcept
									   { return const_Block(*this, loc); }
	// clang-format on

private:
	std::array<T, Size> board_{};
}; // class Board

//===--- free-functions ---------------------------------------------------===//

template<typename T, int N>
constexpr bool operator!=(Board<T, N> const& left, Board<T, N> const& right)
{
	return !(left == right);
}

//===----------------------------------------------------------------------===//
// Board - member-functions
//===----------------------------------------------------------------------===//
// Board - Constructors
template<typename T, int N>
constexpr Board<T, N>::Board() noexcept
{
	valid_dimensions<N>();
}

template<typename T, int N>
constexpr Board<T, N>::Board(const T& default_value)
{
	valid_dimensions<N>();
	board_.fill(default_value);
}

template<typename T, int N>
Board<T, N>::Board(std::initializer_list<T> list)
{
	valid_dimensions<N>();
	if (list.size() != size_t{full_size<N>})
		throw std::length_error{"Invalid length initializer_list"};
	std::copy(std::begin(list), std::end(list), std::begin(board_));
}

//===----------------------------------------------------------------------===//
// properties
template<typename T, int N>
bool Board<T, N>::operator==(const Board& other) const
{
	return board_ == other.board_;
}

//===----------------------------------------------------------------------===//
// Board - Operations
template<typename T, int N>
inline void Board<T, N>::clear()
{ // all elements to the empty value
	board_.fill(T{});
}

//===----------------------------------------------------------------------===//
// Board - element access
template<typename T, int N>
constexpr T& Board<T, N>::at(const Location loc)
{
	if (!is_valid<N>(loc))
	{
		throw error::invalid_Location{"Board::at(Location)"};
	}
	// NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
	return board_[gsl::narrow_cast<size_t>(loc.element())];
}

template<typename T, int N>
constexpr T const& Board<T, N>::at(const Location loc) const
{
	if (!is_valid<N>(loc))
	{
		throw error::invalid_Location{"Board::at(Location) const"};
	}
	// NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
	return board_[gsl::narrow_cast<size_t>(loc.element())];
}

template<typename T, int N>
constexpr T& Board<T, N>::at(const index row, const index col)
{
	if (!is_valid_size<N>(row, col))
	{
		throw error::invalid_Location{"Board::at(int row, col)"}; // <stdexcept>
	}
	// NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
	return board_[gsl::narrow_cast<size_t>(Location(row, col).element())];
}

template<typename T, int N>
constexpr T const& Board<T, N>::at(const index row, const index col) const
{
	if (!is_valid_size<N>(row, col))
	{
		throw error::invalid_Location{"Board::at(row, col) const"};
	}
	// NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
	return board_[gsl::narrow_cast<size_t>(Location(row, col).element())];
}

// deprecated
template<typename T, int N>
T& Board<T, N>::at(const index elem)
{
	if (!is_valid_size<N>(elem))
	{
		throw error::invalid_Location{"Board::at(int)"};
	}
	return board_.at(gsl::narrow_cast<size_t>(elem));
}

// deprecated
template<typename T, int N>
const T& Board<T, N>::at(const index elem) const
{
	if (!is_valid_size<N>(elem))
	{
		throw error::invalid_Location{"Board::at(int) const"};
	}
	return board_.at(gsl::narrow_cast<size_t>(elem));
}

template<typename T, int N>
constexpr T& Board<T, N>::operator[](const Location loc) noexcept
{
	// NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
	return board_[gsl::narrow_cast<size_t>(loc.element())];
}

template<typename T, int N>
constexpr T const& Board<T, N>::operator[](const Location loc) const noexcept
{
	// NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-constant-array-index)
	return board_[gsl::narrow_cast<size_t>(loc.element())];
}

//===----------------------------------------------------------------------===//
// Board - Iterators
template<typename T, int N>
constexpr typename Board<T, N>::iterator Board<T, N>::begin() noexcept
{
	return Board_iterator<T, N>(gsl::not_null<Board*>{this});
}
template<typename T, int N>
constexpr typename Board<T, N>::iterator Board<T, N>::end() noexcept
{
	return Board_iterator<T, N>(
		gsl::not_null<Board*>{this}, Location{full_size<N>});
}
template<typename T, int N>
constexpr typename Board<T, N>::const_iterator
	Board<T, N>::cbegin() const noexcept
{
	return const_Board_iterator<T, N>(gsl::not_null<Board const*>{this});
}
template<typename T, int N>
constexpr typename Board<T, N>::const_iterator
	Board<T, N>::cend() const noexcept
{
	return const_Board_iterator<T, N>(
		gsl::not_null<Board const*>{this}, Location{full_size<N>});
}
template<typename T, int N>
constexpr typename Board<T, N>::reverse_iterator Board<T, N>::rbegin() noexcept
{
	return reverse_Board_iterator<T, N>(gsl::not_null<Board*>{this});
}
template<typename T, int N>
constexpr typename Board<T, N>::reverse_iterator Board<T, N>::rend() noexcept
{
	return reverse_Board_iterator<T, N>(
		gsl::not_null<Board*>{this}, Location{-1});
}
template<typename T, int N>
constexpr typename Board<T, N>::const_reverse_iterator
	Board<T, N>::crbegin() const noexcept
{
	return const_reverse_Board_iterator<T, N>(
		gsl::not_null<Board const*>{this});
}
template<typename T, int N>
constexpr typename Board<T, N>::const_reverse_iterator
	Board<T, N>::crend() const noexcept
{
	return const_reverse_Board_iterator<T, N>(
		gsl::not_null<Board const*>{this}, Location{-1});
}

} // namespace Sudoku

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
/* randist/beta.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 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 <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>

/* The beta distribution has the form

   p(x) dx = (Gamma(a + b)/(Gamma(a) Gamma(b))) x^(a-1) (1-x)^(b-1) dx

   The method used here is the one described in Knuth */

double
gsl_ran_beta (const gsl_rng * r, const double a, const double b)
{
  if ( (a <= 1.0) && (b <= 1.0) )
    {
      double U, V, X, Y;
      while (1)
        {
          U = gsl_rng_uniform_pos(r);
          V = gsl_rng_uniform_pos(r);
          X = pow(U, 1.0/a);
          Y = pow(V, 1.0/b);
          if ((X + Y ) <= 1.0)
            {
              if (X + Y > 0)
                {
                  return X/ (X + Y);
                }
              else
                {
                  double logX = log(U)/a;
                  double logY = log(V)/b;
                  double logM = logX > logY ? logX: logY;
                  logX -= logM;
                  logY -= logM;
                  return exp(logX - log(exp(logX) + exp(logY)));
                }
            }
        }
    }
  else
    {
      double x1 = gsl_ran_gamma (r, a, 1.0);
      double x2 = gsl_ran_gamma (r, b, 1.0);
      return x1 / (x1 + x2);
    }
}

double
gsl_ran_beta_pdf (const double x, const double a, const double b)
{
  if (x < 0 || x > 1)
    {
      return 0 ;
    }
  else
    {
      double p;

      double gab = gsl_sf_lngamma (a + b);
      double ga = gsl_sf_lngamma (a);
      double gb = gsl_sf_lngamma (b);

      if (x == 0.0 || x == 1.0)
        {
	  if (a > 1.0 && b > 1.0)
	    {
	      p = 0.0;
	    }
	  else
	    {
	      p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
	    }
        }
      else
        {
          p = exp (gab - ga - gb + log(x) * (a - 1)  + log1p(-x) * (b - 1));
        }

      return p;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
# ifndef __HEADER__GA__
# define __HEADER__GA__

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include "utils.h"

# define UL_SIZE sizeof(unsigned long)

# define GENES_C1 3
# define GENES_C2 11

// we are going to use a two chromosome genome,
typedef struct Genome {
	unsigned long c1[GENES_C1];
	unsigned long c2[GENES_C2];
	double fitness;
} Genome;

void generate_genome(Genome * genome);

Genome * generate_population(int individuals);
int next_generation(
	Genome * parents, Genome * children,
	int n_elitism, int n_select, int n_cross, int n_new, double p_mutation, int mutation_bit
);

/*
 * Copies the information in the input genome to the output genome.
 * Note that the fitness function is also copied over.
 */
void copy_genome(Genome * in, Genome * out);

// a.k.a crossover
void tinder(Genome * population, int pop_size, Genome * out);
// actual crossover
void crossover_genomes(
	Genome * gen1in,
	Genome * gen2in,
	Genome * out
);

void mutate_genome(Genome * genome, double p_mut);

/*
 * Implementation of elitism operator.
 * This selects the best of the best individuals to be passed over to the next generation.
 *
 * @param population the initial population with fitness values calculated
 * @param pop_size the total population size
 * @param number_elitism the number of best individuals to be selected
 * @param out the population to be filled with the best individuals
 *
 * @return the position of the best individual
 */

int extinction(int ek, Genome * population, Genome * survivors, int pop_size, int number_survivors);

void elitism(Genome * population, int pop_size, int number_elitism, Genome * out);

/*
 * Implementation of roulette wheel method.
 * This method chooses the individuals according to their fidelity, individuals with
 * lower fidelity have a lower change of being selected.
 *
 * @param population the initial population with fitness values calculated
 * @param pop_size the total population size
 * @param best_genomes the number of best individuals to choose
 * @param out the population to be filled with the best individuals
 */
int casting(Genome * population, int pop_size, int best_genomes, Genome * out);

/*
 * Add `n_new` individuals to the start of the passed genome array.
 */
void migration(
	Genome * genome,
	int n_new
);

# endif


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

/* Author:  G. Jungman */

#ifndef __GSL_SF_BESSEL_H__
#define __GSL_SF_BESSEL_H__

#include <stdlib.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_precision.h>
#include <gsl/gsl_sf_result.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

/* Regular Bessel Function J_0(x)
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_bessel_J0_e(const double x,  gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_J0(const double x);

/* Regular Bessel Function J_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_J1_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_J1(const double x);

/* Regular Bessel Function J_n(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Jn_e(int n, double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Jn(const int n, const double x);

/* Regular Bessel Function J_n(x),  nmin <= n <= nmax
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array);

/* Irregular Bessel function Y_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Y0_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Y0(const double x);

/* Irregular Bessel function Y_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Y1_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Y1(const double x);

/* Irregular Bessel function Y_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Yn_e(int n,const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Yn(const int n,const double x);

/* Irregular Bessel function Y_n(x), nmin <= n <= nmax
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Yn_array(const int nmin, const int nmax, const double x, double * result_array);

/* Regular modified Bessel function I_0(x)
 *
 * exceptions: GSL_EOVRFLW
 */
GSL_EXPORT int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_I0(const double x);

/* Regular modified Bessel function I_1(x)
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_I1(const double x);

/* Regular modified Bessel function I_n(x)
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_In_e(const int n, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_In(const int n, const double x);

/* Regular modified Bessel function  I_n(x) for n=nmin,...,nmax
 *
 * nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_In_array(const int nmin, const int nmax, const double x, double * result_array);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_0(x)
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_I0_scaled(const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_I1_scaled(const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_n(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_In_scaled_e(int n, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_In_scaled(const int n, const double x);

/* Scaled regular modified Bessel function
 *  exp(-|x|) I_n(x)  for n=nmin,...,nmax
 *
 * nmin >=0, nmax >= nmin
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_In_scaled_array(const int nmin, const int nmax, const double x, double * result_array);

/* Irregular modified Bessel function K_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_K0(const double x);

/* Irregular modified Bessel function K_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_K1(const double x);

/* Irregular modified Bessel function K_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Kn_e(const int n, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Kn(const int n, const double x);

/* Irregular modified Bessel function  K_n(x)  for n=nmin,...,nmax
 *
 * x > 0.0, nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EOVRFLW, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Kn_array(const int nmin, const int nmax, const double x, double * result_array);

/* Scaled irregular modified Bessel function
 *  exp(x) K_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM
 */
GSL_EXPORT int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_K0_scaled(const double x);

/* Scaled irregular modified Bessel function
 *  exp(x) K_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_K1_scaled(const double x);

/* Scaled irregular modified Bessel function
 *  exp(x) K_n(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Kn_scaled_e(int n, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Kn_scaled(const int n, const double x);

/* Scaled irregular modified Bessel function  exp(x) K_n(x)  for n=nmin,...,nmax
 *
 * x > 0.0, nmin >=0, nmax >= nmin
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Kn_scaled_array(const int nmin, const int nmax, const double x, double * result_array);

/* Regular spherical Bessel function j_0(x) = sin(x)/x
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_bessel_j0_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_j0(const double x);

/* Regular spherical Bessel function j_1(x) = (sin(x)/x - cos(x))/x
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_j1_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_j1(const double x);

/* Regular spherical Bessel function j_2(x) = ((3/x^2 - 1)sin(x) - 3cos(x)/x)/x
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_j2_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_j2(const double x);

/* Regular spherical Bessel function j_l(x)
 *
 * l >= 0, x >= 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_jl_e(const int l, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_jl(const int l, const double x);

/* Regular spherical Bessel function j_l(x) for l=0,1,...,lmax
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_jl_array(const int lmax, const double x, double * result_array);

/* Regular spherical Bessel function j_l(x) for l=0,1,...,lmax
 * Uses Steed's method.
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_jl_steed_array(const int lmax, const double x, double * jl_x_array);

/* Irregular spherical Bessel function y_0(x)
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_bessel_y0_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_y0(const double x);

/* Irregular spherical Bessel function y_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_y1_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_y1(const double x);

/* Irregular spherical Bessel function y_2(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_y2_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_y2(const double x);

/* Irregular spherical Bessel function y_l(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_yl_e(int l, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_yl(const int l, const double x);

/* Irregular spherical Bessel function y_l(x) for l=0,1,...,lmax
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_yl_array(const int lmax, const double x, double * result_array);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_0(x)
 *
 * exceptions: none
 */
GSL_EXPORT int gsl_sf_bessel_i0_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_i0_scaled(const double x);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_1(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_i1_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_i1_scaled(const double x);

/* Regular scaled modified spherical Bessel function
 *
 * Exp[-|x|] i_2(x)
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_i2_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_i2_scaled(const double x);

/* Regular scaled modified spherical Bessel functions
 *
 * Exp[-|x|] i_l(x)
 *
 * i_l(x) = Sqrt[Pi/(2x)] BesselI[l+1/2,x]
 *
 * l >= 0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_il_scaled_e(const int l, double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_il_scaled(const int l, const double x);

/* Regular scaled modified spherical Bessel functions
 *
 * Exp[-|x|] i_l(x)
 * for l=0,1,...,lmax
 *
 * exceptions: GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array);

/* Irregular scaled modified spherical Bessel function
 * Exp[x] k_0(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_k0_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_k0_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_1(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_EXPORT int gsl_sf_bessel_k1_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_k1_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_2(x)
 *
 * x > 0.0
 * exceptions: GSL_EDOM, GSL_EUNDRFLW, GSL_EOVRFLW
 */
GSL_EXPORT int gsl_sf_bessel_k2_scaled_e(const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_k2_scaled(const double x);

/* Irregular modified spherical Bessel function
 * Exp[x] k_l[x]
 *
 * k_l(x) = Sqrt[Pi/(2x)] BesselK[l+1/2,x]
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_kl_scaled(const int l, const double x);

/* Irregular scaled modified spherical Bessel function
 * Exp[x] k_l(x)
 *
 * for l=0,1,...,lmax
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_kl_scaled_array(const int lmax, const double x, double * result_array);

/* Regular cylindrical Bessel function J_nu(x)
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Jnu_e(const double nu, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Jnu(const double nu, const double x);

/* Irregular cylindrical Bessel function Y_nu(x)
 *
 * exceptions:
 */
GSL_EXPORT int gsl_sf_bessel_Ynu_e(double nu, double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Ynu(const double nu, const double x);

/* Regular cylindrical Bessel function J_nu(x)
 * evaluated at a series of x values. The array
 * contains the x values. They are assumed to be
 * strictly ordered and positive. The array is
 * over-written with the values of J_nu(x_i).
 *
 * exceptions: GSL_EDOM, GSL_EINVAL
 */
GSL_EXPORT int gsl_sf_bessel_sequence_Jnu_e(double nu, gsl_mode_t mode, size_t size, double * v);

/* Scaled modified cylindrical Bessel functions
 *
 * Exp[-|x|] BesselI[nu, x]
 * x >= 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
GSL_EXPORT int gsl_sf_bessel_Inu_scaled_e(double nu, double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Inu_scaled(double nu, double x);

/* Modified cylindrical Bessel functions
 *
 * BesselI[nu, x]
 * x >= 0, nu >= 0
 *
 * exceptions: GSL_EDOM, GSL_EOVRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Inu_e(double nu, double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Inu(double nu, double x);

/* Scaled modified cylindrical Bessel functions
 *
 * Exp[+|x|] BesselK[nu, x]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
GSL_EXPORT int gsl_sf_bessel_Knu_scaled_e(const double nu, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Knu_scaled(const double nu, const double x);

/* Modified cylindrical Bessel functions
 *
 * BesselK[nu, x]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM, GSL_EUNDRFLW
 */
GSL_EXPORT int gsl_sf_bessel_Knu_e(const double nu, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_Knu(const double nu, const double x);

/* Logarithm of modified cylindrical Bessel functions.
 *
 * Log[BesselK[nu, x]]
 * x > 0, nu >= 0
 *
 * exceptions: GSL_EDOM
 */
GSL_EXPORT int gsl_sf_bessel_lnKnu_e(const double nu, const double x, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_lnKnu(const double nu, const double x);

/* s'th positive zero of the Bessel function J_0(x).
 *
 * exceptions:
 */
GSL_EXPORT int gsl_sf_bessel_zero_J0_e(unsigned int s, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_zero_J0(unsigned int s);

/* s'th positive zero of the Bessel function J_1(x).
 *
 * exceptions:
 */
GSL_EXPORT int gsl_sf_bessel_zero_J1_e(unsigned int s, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_zero_J1(unsigned int s);

/* s'th positive zero of the Bessel function J_nu(x).
 *
 * exceptions:
 */
GSL_EXPORT int gsl_sf_bessel_zero_Jnu_e(double nu, unsigned int s, gsl_sf_result * result);
GSL_EXPORT double gsl_sf_bessel_zero_Jnu(double nu, unsigned int s);

__END_DECLS

#endif /* __GSL_SF_BESSEL_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/**
 *
 * @file core_blas.h
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Jakub Kurzak
 * @author Hatem Ltaief
 * @date 2010-11-15
 *
 **/
#ifndef _PLASMA_CORE_BLAS_H_
#define _PLASMA_CORE_BLAS_H_

#ifdef USE_MKL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "plasmatypes.h"
#include "descriptor.h"

#include "core_dblas.h"

#ifdef __cplusplus
extern "C" {
#endif

  /*
   * Coreblas Error
   */
#define coreblas_error(k, str) fprintf(stderr, "%s: Parameter %d / %s\n", __func__, k, str);

 /** ****************************************************************************
  *  LAPACK Constants
  **/
extern char *plasma_lapack_constants[];
#define lapack_const(plasma_const) plasma_lapack_constants[plasma_const][0]

  /*
   * CBlas enum
   */
#define CBLAS_TRANSPOSE int
#define CBLAS_UPLO      int
#define CBLAS_DIAG      int
#define CBLAS_SIDE      int
//#define CBLAS_TRANSPOSE enum CBLAS_TRANSPOSE
//#define CBLAS_UPLO      enum CBLAS_UPLO
//#define CBLAS_DIAG      enum CBLAS_DIAG
//#define CBLAS_SIDE      enum CBLAS_SIDE

#ifdef __cplusplus
}
#endif

#endif /* _PLASMA_CORE_BLAS_H_ */


text:
/* Author: Romain "Artefact2" Dal Maso <artefact2@gmail.com> */

/* This program is free software. It comes without any warranty, to the
 * extent permitted by applicable law. You can redistribute it and/or
 * modify it under the terms of the Do What The Fuck You Want To Public
 * License, Version 2, as published by Sam Hocevar. See
 * http://sam.zoy.org/wtfpl/COPYING for more details. */

#define warn(...) do {	  \
		char msg[256]; \
		snprintf(msg, 256, __VA_ARGS__); \
		fprintf(stderr, "%s(): %s\n", __func__, msg); \
	} while(0)

#define fatal(...) do {	  \
		warn(__VA_ARGS__); \
		exit(1); \
	} while(0)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>

#define MCTS (1 << 15)

void read_rolls(unsigned int, unsigned int* rcount, unsigned int* n);
unsigned int ecdf_distance(unsigned int, unsigned int*);
void ecdf_gen_mc_table(unsigned int sides, unsigned int n, unsigned int count, unsigned int*);
bool find_in_sorted_array(unsigned int val, unsigned int count, unsigned int* array, unsigned int* low, unsigned int* high);
double chisq_test(unsigned int sides, unsigned int n, unsigned int*);
void confidence_test(unsigned int sides, unsigned int n, unsigned int*);

int main(int argc, char** argv) {
	if(argc != 2) {
		fprintf(stderr, "Usage: %s <num-sides> < rolls.txt\n", argv[0]);
		exit(2);
	}

	unsigned int sides = strtol(argv[1], NULL, 10);
	if(sides < 2) {
		fatal("must have at least 2 sides, got %d", sides);
	}

	unsigned int rcount[sides];
	unsigned int n;
	read_rolls(sides, rcount, &n);
	if(n % sides) {
		fatal("sample size (%d) not a multiple of sides (%d)", n, sides);
	}

	printf("SmpSize:   n=%d\n", n);

	unsigned int ecdf_table[MCTS];
	ecdf_gen_mc_table(sides, n, MCTS, ecdf_table);

	unsigned int dist = ecdf_distance(sides, rcount);
	unsigned int lo, hi;
	bool found;

	found = find_in_sorted_array(dist, MCTS, ecdf_table, &lo, &hi);
	printf("ECDF:      p%s%5.4f\n", found ? "=" : "<", found ? (1.0 - (double)(lo + hi) / (double)(2 * MCTS)) : ((double)(MCTS - lo) / (double)MCTS));

	printf("ChiSq:     p=%5.4f\n", chisq_test(sides, n, rcount));

	printf("ConfInt99:");
	confidence_test(sides, n, rcount);
	printf("\n");
}

void read_rolls(unsigned int sides, unsigned int* rcount, unsigned int* n) {
	char buf[16];
	unsigned int roll;

	memset(rcount, 0, sides * sizeof(unsigned int));
	*n = 0;

	while(!feof(stdin)) {
		if(fgets(buf, 16, stdin) == NULL) continue;
		if(buf[0] == '\n' || buf[0] == '\0') continue;

		roll = strtol(buf, NULL, 10);
		if(roll < 1 || roll > sides) {
			warn("ignoring roll %d", roll);
			continue;
		}

		++(rcount[roll - 1]);
		++(*n);
	}
}

unsigned int ecdf_distance(unsigned int sides, unsigned int* rcount) {
	unsigned int cumulative = 0, ideal_cumulative, ideal_step, i, distance, tdist;

	for(i = 0; i < sides; ++i) {
		cumulative += rcount[i];
	}

	ideal_cumulative = cumulative;
	ideal_step = ideal_cumulative / sides;
	distance = 0;

	for(i = 0; i < sides; ++i) {
		cumulative -= rcount[i];
		ideal_cumulative -= ideal_step;
		tdist = (cumulative >= ideal_cumulative) ? (cumulative - ideal_cumulative) : (ideal_cumulative - cumulative);
		if(tdist > distance) distance = tdist;
	}

	return distance;
}

int cmpuint(const void* a, const void* b) {
	return *(unsigned int*)a - *(unsigned int*)b;
}

void ecdf_gen_mc_table(unsigned int sides, unsigned int n, unsigned int count, unsigned int* table) {
	unsigned int rcount[sides];
	unsigned int i, j;

	for(i = 0; i < count; ++i) {
		memset(rcount, 0, sides * sizeof(unsigned int));

		for(j = 0; j < n; ++j) {
			++(rcount[rand() % sides]);
		}

		table[i] = ecdf_distance(sides, rcount);
	}

	qsort(table, count, sizeof(unsigned int), cmpuint);
}

bool find_in_sorted_array(unsigned int val, unsigned int count, unsigned int* array, unsigned int* low, unsigned int* high) {
	if(array[0] > val) {
		*low = *high = 0;
		return false;
	} else if(array[count - 1] < val) {
		*low = *high = count - 1;
		return false;
	}

	unsigned int mid;

	*low = 0;
	*high = count - 1;

	while(*high > *low) {
		mid = *low + (*high - *low) / 2;

		if(array[mid] == val) {
			*low = *high = mid;
			while(*low < count && array[*low] == val) --(*low);
			++(*low);
			while(*high < count && array[*high] == val) ++high;
			--(*high);
			return true;
		} else if(array[mid] > val) {
			*high = mid;
		} else {
			*low = mid;
		}
	}

	return false;
}

double chisq_test(unsigned int sides, unsigned int n, unsigned int* rcount) {
	unsigned int ideal = n / sides;
	double chisq = 0.0;

	for(unsigned int i = 0; i < sides; ++i) {
		chisq += (rcount[i] - ideal) * (rcount[i] - ideal);
	}
	chisq /= (double)ideal;

	return 1.0 - gsl_sf_gamma_inc_P((double)(sides - 1) / (double)2, chisq / (double)2);
}

void confidence_test(unsigned int sides, unsigned int n, unsigned int* rcount) {
	double f, amp, ideal = 1.0 / (double)sides;
	bool anomalies = false;

	for(unsigned int i = 0; i < sides; ++i) {
		f = (double)rcount[i] / (double)n;
		amp = 2.575 * sqrt(f * (1.0 - f) / (double)n);
		if(ideal > f + amp) {
			printf(" %d-", i+1);
			anomalies = true;
		} else if(ideal < f - amp) {
			printf(" %d+", i+1);
			anomalies = true;
		}
	}

	if(!anomalies) {
		printf(" OK");
	}
}


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

#ifndef __GSL_CHEBYSHEV_H__
#define __GSL_CHEBYSHEV_H__

#include <gsl/gsl_mode.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

/* data for a Chebyshev series over a given interval */

struct gsl_cheb_series_struct {

  double * c;   /* coefficients                */
  size_t order; /* order of expansion          */
  double a;     /* lower interval point        */
  double b;     /* upper interval point        */

  /* The following exists (mostly) for the benefit
   * of the implementation. It is an effective single
   * precision order, for use in single precision
   * evaluation. Users can use it if they like, but
   * only they know how to calculate it, since it is
   * specific to the approximated function. By default,
   * order_sp = order.
   * It is used explicitly only by the gsl_cheb_eval_mode
   * functions, which are not meant for casual use.
   */
  size_t order_sp;

  /* Additional elements not used by specfunc */

  double * f;   /* function evaluated at chebyschev points  */
};
typedef struct gsl_cheb_series_struct gsl_cheb_series;

/* Calculate a Chebyshev series of specified order over
 * a specified interval, for a given function.
 * Return 0 on failure.
 */
gsl_cheb_series * gsl_cheb_alloc(const size_t order);

/* Free a Chebyshev series previously calculated with gsl_cheb_alloc().
 */
void gsl_cheb_free(gsl_cheb_series * cs);

/* Calculate a Chebyshev series using the storage provided.
 * Uses the interval (a,b) and the order with which it
 * was initially created.
 *
 */
int gsl_cheb_init(gsl_cheb_series * cs, const gsl_function * func,
                  const double a, const double b);

/* Evaluate a Chebyshev series at a given point.
 * No errors can occur for a struct obtained from gsl_cheb_new().
 */
double gsl_cheb_eval(const gsl_cheb_series * cs, const double x);
int gsl_cheb_eval_err(const gsl_cheb_series * cs, const double x,
                      double * result, double * abserr);

/* Evaluate a Chebyshev series at a given point, to (at most) the given order.
 * No errors can occur for a struct obtained from gsl_cheb_new().
 */
double gsl_cheb_eval_n(const gsl_cheb_series * cs, const size_t order,
                       const double x);
int gsl_cheb_eval_n_err(const gsl_cheb_series * cs, const size_t order,
                        const double x, double * result, double * abserr);

/* Evaluate a Chebyshev series at a given point, using the default
 * order for double precision mode(s) and the single precision
 * order for other modes.
 * No errors can occur for a struct obtained from gsl_cheb_new().
 */
double gsl_cheb_eval_mode(const gsl_cheb_series * cs, double x, gsl_mode_t mode);
int gsl_cheb_eval_mode_e(const gsl_cheb_series * cs, const double x, gsl_mode_t mode, double * result, double * abserr);

/* Compute the derivative of a Chebyshev series.
 */
int gsl_cheb_calc_deriv(gsl_cheb_series * deriv, const gsl_cheb_series * cs);

/* Compute the integral of a Chebyshev series. The
 * integral is fixed by the condition that it equals zero at
 * the left end-point, ie it is precisely
 *       Integrate[cs(t; a,b), {t, a, x}]
 */
int gsl_cheb_calc_integ(gsl_cheb_series * integ, const gsl_cheb_series * cs);

__END_DECLS

#endif /* __GSL_CHEBYSHEV_H__ */


text:
/*
 * File:   particle.h
 * Author: adithya
 *
 * Created on December 15, 2014, 7:23 PM
 */

#ifndef PARTICLE_H
#define	PARTICLE_H

//#include "myVector.h"
//#include "quaternion.h"
//#include "species.h"
//#include "simulationParameters.h"
//
//#include <fstream>
//#include <sstream>
//#include <vector>
//
//#include <gsl/gsl_rng.h>
//#include <gsl/gsl_randist.h>

#include "myVector.h"
#include "quaternion.h"
#include "main.h"

class particle
{
public:
    particle();
    particle(const particle& orig);
    virtual ~particle();
public:
    int particleID;/**This ID helps identifying the particle in question*/
    int speciesType;/**Is the species type of the particle, which in turn holds information about the particles' physical properties.*/
    int particleType;/**Describes the type of the particle, 1 if a particle with a domain(eGFRD), 2 if LD*/
    //std::vector < std::pair <int, int> > complexList;

    std::vector < std::pair <int, int> > complexList;/**Gives the list of all other particles (if any)
                                                      connectected to this particle in a complex
                                                      and the patches of the particle that are occupied*/
    /**
     * @memberVariableOfSpeciesClass is a vector that contains the
     * the state of all patches
     * for a substrate
     * 0 is unphosphorylated
     * 1 is phosphorylated
     * for an enzyme its always 1
     */
    std::vector<int> patchState;

    myVector particlePosition;  /**3D position of the centre of the particle*/
    myVector force;             /**force on the particle*/

    quaternion q;               /**The orientation of the particle in 3 dimensions*/
    quaternion torque;          /**The torque of the particle*/

    std::vector<myVector> patchPos;          /** 3D position of the centre of the patch of the particle which is used to
                                 visualize the rotation*/

    /*GREENS FUNCTIONS RELATED VARIABLES*/
    myVector nextEscapePosition;/**The center of the particle when the particle escapes*/
    myVector burstPosition;/**center of the particle when the domain bursts*/
    double R;/**This is the radius of displacement of the particle from the original position, after bursting*/

     /*DOMAIN RELATED MEMBER VARIABLES*/
    int nextReactionType;/**if 0 then nothing (just Brownian step)
                            if 1 then monomolecular dissociation
                            if 2 hopping
                            if 3 phosphorylation/dephosphorylation*/
    int nextEventType;/**If 0 its escape(from the domain)
                         if 1 its reaction
                         */
    double nextEscapeTime;/**The time at which the next escape occurs*/
    double nextReactionTime;/**The time at which the next reaction occurs
                             minimum of dissociation time, hopping time and phosphorylation time*/
    double nextEventTime;/**The time at which the next event occurs, the min of the above two values*/
    double timeOfConstructionOfDomain;/**Global time at which the domain is constructed*/
    double domainRadius;/**radius of the domain around the particle*/

    bool operator < (const particle temp) const
    {
        return (this->nextEventTime < temp.nextEventTime);
    }

     /**Function prints out the particle class
     */
    friend std::ostream& operator<< (std::ostream &out, particle &part);

    /**
     This function draws a random position at the sphere of domain radius if the particle escapes
     */
    void drawNextEscapePosition();

    /**
     This function draws a new position for the particle if the domain is prematurely burst
     */
    void drawBurstPosition (double time);

    /**
     This function gets a next event time and type using greens functions for a single domain
     */
    void drawNextEventTimeAndType ();

    /**
     This function builds a domain on the particle which calls the function

     @param nearestNeighborDistance radius of the domain to be built
     */
    void buildADomain(double nearestNeighborDistance);

    /**
     This function bursts the given domain it can be used for bursting a domain, moving a particle or converting a particle to a LD particle

     @param instruction escape, burst or convert
     */
    void burstTheDomain(std::string instruction);

    /**
     This dissociates a particle in a domain
     */
    void dissociateParticleFromDomain();

};/**Particle holds all information related to the particle in the simulation*/

#endif	/* PARTICLE_H */


text:
#ifndef _GP_H_
#define _GP_H_

#include "armadillo"
using namespace arma;
#include "myTypes.h"

#include "KernelFunction.h"
#include "MeanFunction.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_multimin.h>

double f_eval_mean(const gsl_vector *x, void *param);
void df_eval_mean(const gsl_vector *x, void *param, gsl_vector *g);
void fdf_eval_mean(const gsl_vector *x, void *param, double *f, gsl_vector *g);

double f_eval_kernel(const gsl_vector *x, void *param);
void df_eval_kernel(const gsl_vector *x, void *param, gsl_vector *g);
void fdf_eval_kernel(const gsl_vector *x, void *param, double *f, gsl_vector *g);

double f_eval_noise(const gsl_vector *x, void *param);
void df_eval_noise(const gsl_vector *x, void *param, gsl_vector *g);
void fdf_eval_noise(const gsl_vector *x, void *param, double *f, gsl_vector *g);

class GP {
public:
  GP(REAL s2_n, KernelFunction *kernel, MeanFunction *mean);

  void SetTraining(const Mat<REAL>& X, const Row<REAL> &y);
  void AddTraining(const Mat<REAL>& X, const Row<REAL> &y);
  void Predict(const Mat<REAL> &Xs, Row<REAL> &mu);
  void Predict(const Mat<REAL> &Xs, Row<REAL> &mu, Row<REAL> &var);
  void Predict(const Col<REAL> &Xs, REAL &mu);
  void Predict(const Col<REAL> &Xs, REAL &mu, REAL &var);

  void PredictGradient(const Col<REAL> &Xs, Col<REAL> &grad);
  void PredictGradient(const Col<REAL> &Xs, Col<REAL> &grad, Col<REAL> &vargrad);

  void ComputeAlpha();
  void ComputeChol();
  void ComputeW();
  REAL ComputeLikelihood();

  void SetKernelFuncParams(const Col<REAL>& param);
  void SetMeanFuncParams(const Col<REAL>& param);
  void SetNoise(const REAL &s2_n);

  void GradLikelihoodMeanParams(Col<REAL> &grad);
  void GradLikelihoodKernelParams(Col<REAL> &grad);
  void GradLikelihoodNoise(Col<REAL> &grad);

  void HessianLikelihoodMeanParams(Mat<REAL> &hessian);

  void OptimizeNoiseParam(REAL &noise_param, int max_iterations=10);
  void OptimizeMeanParam(Col<REAL> &mean_param, int max_iterations=10);
  void OptimizeKernelParam(Col<REAL> &kernel_param, int max_iterations=10);

  inline Mat<REAL> &GetTrainingData() { return X; }
  inline KernelFunction *GetKernelFunction() { return this->kernel; }
  inline MeanFunction *GetMeanFunction() { return this->mean; }

protected:
  KernelFunction *kernel;
  MeanFunction *mean;
  Mat<REAL> K;
  Mat<REAL> W;
  Mat<REAL> X;
  Mat<REAL> L;
  Row<REAL> y;
  Row<REAL> meanvals;
  Col<REAL> alpha;
  REAL s2_n;
  REAL loglikelihood;

  bool need_to_compute_alpha;
  bool need_to_compute_chol;
  bool need_to_compute_w;

  void MatrixMap(Mat<REAL> &matrix, const Mat<REAL> &a, const Mat<REAL> &b);
};

#endif

max_stars_repo_name: jonfink/GP

max_issues_repo_name: RobustFieldAutonomyLab/GP

max_forks_repo_name: RobustFieldAutonomyLab/GP


text:
#ifndef PyGSL_COMPLEX_HELPERS_H
#define PyGSL_COMPLEX_HELPERS_H 1
#include <pygsl/utils.h>
#include <pygsl/intern.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_errno.h>

/*
 *  These macros convert a PyObect to  a complex if numeric
 *  Input recieved. If complex, accessing the data, otherwise
 *  reverting to a function.
 */
/* -------------------------------------------------------------------------
   Helper Functions
   ------------------------------------------------------------------------- */
PyGSL_API_EXTERN int
PyGSL_PyComplex_to_gsl_complex(PyObject * src, gsl_complex * mycomplex);

PyGSL_API_EXTERN int
PyGSL_PyComplex_to_gsl_complex_float(PyObject * src,
				     gsl_complex_float * mycomplex);
PyGSL_API_EXTERN int
PyGSL_PyComplex_to_gsl_complex_long_double(PyObject * src,
					   gsl_complex_long_double * mycomplex);

#ifndef _PyGSL_API_MODULE
#define PyGSL_PyComplex_to_gsl_complex \
(*(int (*) (PyObject *, gsl_complex *))              PyGSL_API[PyGSL_PyComplex_to_gsl_complex_NUM])
#define PyGSL_PyComplex_to_gsl_complex_float \
(*(int (*) (PyObject *, gsl_complex_float *))        PyGSL_API[PyGSL_PyComplex_to_gsl_complex_float_NUM])
#define PyGSL_PyComplex_to_gsl_complex_long_double \
(*(int (*) (PyObject *, gsl_complex_long_double *))  PyGSL_API[PyGSL_PyComplex_to_gsl_complex_long_double_NUM])
#endif /* _PyGSL_API_MODULE */

#define PyGSL_PyCOMPLEX_TO_gsl_complex(object, tmp)            \
(PyComplex_Check((object))) ? 				       \
     ((tmp)->dat[0] = ((PyComplexObject *)(object))->cval.real,\
      (tmp)->dat[1] = ((PyComplexObject *)(object))->cval.imag,\
      GSL_SUCCESS) 					       \
     : PyGSL_PyComplex_to_gsl_complex(object, tmp)

#define PyGSL_PyCOMPLEX_TO_gsl_complex_float(object, tmp )     \
(PyComplex_Check((object))) ?                                  \
     ((tmp)->dat[0] = ((PyComplexObject *)(object))->cval.real,\
      (tmp)->dat[1] = ((PyComplexObject *)(object))->cval.imag,\
      GSL_SUCCESS)                                             \
     : PyGSL_PyComplex_to_gsl_complex_float(object, tmp)

#define PyGSL_PyCOMPLEX_TO_gsl_complex_long_double(object, tmp ) \
(PyComplex_Check((object))) ?                                     \
     ((tmp)->dat[0] = ((PyComplexObject *)(object))->cval.real,  \
      (tmp)->dat[1] = ((PyComplexObject *)(object))->cval.imag,  \
      GSL_SUCCESS)                                               \
     : PyGSL_PyComplex_to_gsl_complex_long_double(object, tmp)

#endif  /* PyGSL_COMPLEX_HELPERS_H */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* $Header$ */

/* Purpose: Description (definition) of spherical geometry functions */

/* This file includes BSD-licensed code whose copyright is held by another author
   The copyright owner and license terms for the NCO modifications to that code are
   Copyright (C) 2018--present Charlie Zender
   This file is part of NCO, the netCDF Operators. NCO is free software.
   You may redistribute and/or modify NCO under the terms of the
   3-Clause BSD License with exceptions described in the LICENSE file */

/* This copyright statement and license terms for modification and redistribution
   of the original code were agreed to by the original author, Joseph O'Rourke, on 20190517:

   This code is described in "Computational Geometry in C" (Second Edition),
   Chapter 7.  It is not written to be comprehensible without the
   explanation in that book.

   Written by Joseph O'Rourke.
   Last modified: December 1997
   Questions to jorourke@smith.edu.
   -------------------------------------------------------------------------
   Copyright 1997 by Joseph O'Rourke <jorourke@smith.edu>

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions are
   met:

   1. Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.

   2. 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 HOLDER
   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.
   ------------------------------------------------------------------------- */

/* Usage:
   #include "nco_sph.h" *//* Spherical geometry intersections */

#ifndef NCO_SPH_H /* Contents have not yet been inserted in current source file */
#define NCO_SPH_H

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

#ifdef HAVE_CONFIG_H
# include <config.h> /* Autotools tokens */
#endif /* !HAVE_CONFIG_H */

#ifdef ENABLE_GSL
  #include <gsl/gsl_vector.h>
  #include <gsl/gsl_matrix.h>
  #include <gsl/gsl_permutation.h>
  #include <gsl/gsl_linalg.h>
#endif

/* Personal headers */
#include "nco.h" /* netCDF Operator (NCO) definitions */
#include "nco_mmr.h"     /* Memory management */
#include "nco_omp.h"     /* OpenMP utilities */
#include "nco_rgr.h"     /* Regridding */
#include "nco_sld.h"     /* Swath-Like Data */
#include "nco_sng_utl.h" /* String utilities */

#include "nco_crt.h" /* Cartesian geometry intersections */
#include "nco_ply.h" /* Polygon structure & utilities */

#define NBR_SPH (5)
#define NBR_RLL (5)

#define VP_MAX    1000            /* Max # of pts in polygon */

/* this is 1.0e-20 * PI / 180.0 */
#define ARC_MIN_LENGTH_RAD (1.0e-15)

/* smallest RADIAN */
#define SIGMA_RAD (1.0e-12)

#define SIGMA_TOLERANCE (1.0e-16)
#define DOT_TOLERANCE (1.0e-14)

/* this value plays nice with edges on grids/ne120np4_pentagons.100310.nc */
#define DIST_TOLERANCE (1.0e-14)

/* convert Degrees to Radians */
#define D2R(x)  ((x) * M_PI / 180.0)
/* convert Radians to degrees */
#define R2D(x)  ((x) * 180.0 / M_PI)

/* if true then longitude 0-360 */
/* we need this to convert 3D back to 2D */
#define IS_LON_360 (1)

// #define DEBUG_SPH (1)

#ifdef __cplusplus
/* Use C-bindings so C++-compiled and C-compiled libraries are compatible */
extern "C" {
#endif /* !__cplusplus */

/*---------------------------------------------------------------------
Structures
---------------------------------------------------------------------*/

/* vertex info  */
typedef struct {
  poly_vrl_flg_enm  in_flag;  /* if P_in - vertex from a, Q_in Vertex from b */
  int p_vrt;                  /* if -1 then genuine intersection */
  int q_vrt;
  double p0[NBR_SPH];           /* actual point - used for debuggging */
}vrt_info_sct;

/*---------------------------------------------------------------------
Function prototypes.
---------------------------------------------------------------------*/

int
nco_sph_intersect(poly_sct *P, poly_sct *Q, poly_sct *R, int *r, int flg_snp_to, const char *pq_pre);

char
nco_sph_seg_int_old(double *a, double *b, double *c, double *d, double *p, double *q);

nco_bool
nco_sph_seg_int(double *p0, double *p1, double *q0, double *q1, double *r0, double *r1, int *pq_cross, int flg_snp_to, char *codes);

int
nco_sph_mk_pqcross( double *p0, double *p1, double *pCross, double *q0, double  *q1, double *qCross, int pqCross[], nco_edg_typ_enm rgr_edg_typ,   nco_edg_typ_enm *p_edg_typ, nco_edg_typ_enm *q_edg_typ);

nco_bool
nco_sph_seg_parallel(double *p0, double *p1, double *q0, double *q1, double *r0, double *r1, poly_vrl_flg_enm *inflag, char *codes );

nco_bool
nco_sph_seg_smc(double *p0, double *p1, double *q0, double *q1, double *r0, double *r1, int *pq_cross, int flg_snp_to, char *codes);

nco_bool
nco_sph_seg_vrt_int(double *a, double *b, double *c);

int
nco_sph_lhs(double *Pi, double *Qi);

nco_bool
nco_sph_face(int iLHS, int iRHS, int jRHS);

double
nco_sph_dot(double *a, double *b);

double
nco_sph_dot_nm(double *a, double *b);

double
nco_sph_cross(double *a, double *b, double *c);

double
nco_sph_cross_sub(double *a, double *b, double *c);

double
nco_sph_cross2(double *a, double *b, double *c);

double
nco_sph_trp(double *a, double *b, double *c);

void
nco_sph_add(double *a, double *b, double *c);

void
nco_sph_sub(double *a, double *b, double *c);

void
nco_sph_mlt(double *a, double m);

double
nco_sph_dist(double *a, double *b);

double
nco_sph_rad(double *a);

double
nco_sph_rad2(double *a);

double
nco_sph_sxcross(double *a, double *b, double *c);

void
nco_sph_adi(double *a, double *b);

void
nco_sph_add_pnt(double **R, int *r, double *P);

bool
nco_sph_vrt_info_cmp( vrt_info_sct *info_a, vrt_info_sct *info_b);

void
nco_sph_add_pnt_chk( vrt_info_sct *vrt_info, poly_vrl_flg_enm inflag, int p_vrt, int q_vrt,  double **R, int *r, double *P);

nco_bool
nco_sph_between(double a, double b, double x);

void
nco_sph_prn_pnt(const char *sMsg, double *p, int style, nco_bool bRet);

nco_bool
nco_sph_is_convex(double **sP, int np);

void
nco_sph_prn(double **sR, int r, int istyle);

int
nco_sph_pnt_in_poly(double **sP, int n, double *pControl, double *pVertex);

nco_bool
nco_sph_poly_in_poly(poly_sct *sP,poly_sct *sQ);

void
nco_sph_set_domain(double lon_min_rad, double lon_max_rad, double lat_min_rad, double lat_max_rad);

void
nco_sph_add_lonlat(double *ds);

int
nco_sph_mk_control(poly_sct *sP, nco_bool bInside,  double* pControl  ); /* make a control point that is outside polygon */

nco_bool
nco_sph_intersect_pre(poly_sct *sP,poly_sct *sQ, char *sq_sng  );

int
nco_sph_process_pre(poly_sct *sQ, char *sq_sng, nco_bool *bGenuine);

void
nco_sph_centroid_mk(poly_sct *sP, double *pControl);

nco_bool
nco_sph_inside_mk(poly_sct *sP, double *pControl);

nco_bool
nco_sph_metric( double *p, double *q);

int
nco_sph_metric_int(double *c, double *d, double *Icross);

double
nco_sph_area_karney( double **sP, int np);

double
nco_sph_area_quadrature(double  **sP, int np);

/***************** nco_geo functions these manimpulate lat & lon  ***************************/
void
nco_geo_sph_2_lonlat(double *a, double *lon, double *lat, nco_bool bDeg);

void
nco_geo_lonlat_2_sph(double lon, double lat, double *b, nco_bool bSimple, nco_bool bDeg);

double
nco_geo_lat_correct(double lat1, double lon1, double lon2);

void
nco_geo_get_lat_correct(double lon1, double lat1, double lon2, double lat2, double *dp_min, double *dp_max, nco_bool bDeg);

/**************** functions for RLL grids ***************************************************/
int
nco_rll_intersect(poly_sct *P, poly_sct *Q, poly_sct *R, int *r);

char
nco_rll_seg_int(double *a, double *b, double *c, double *d, double *p, double *q);

nco_bool
nco_rll_seg_parallel(double *p0, double *p1, double *q0, double *q1, double *r0, double *r1, poly_vrl_flg_enm *inflag, char *codes );

void
nco_rll_area(poly_sct *pl);

nco_bool
nco_rll_is_lat_circle(double *p0, double *p1);

int
nco_rll_lhs(double *p0, double *QCross);

int
nco_rll_lhs_lat(double *p0, double *q0, double *q1);

void
nco_rll_add_pnt(double **R, int *r, double *P);

/*********************** functions for matrix*****************************************************/
void
nco_mat_mlt
(double mat[], double vec[], double vec_out[]);

nco_bool
nco_mat_inv(double *mat, double *mat_inv);

nco_bool
nco_mat_int_pl(const double *p0, const double *p1, const double *q0, const double *q1, double *r0);

#ifdef __cplusplus
} /* end extern "C" */
#endif /* !__cplusplus */

#endif /* NCO_SPH_H */

max_stars_repo_name: hbdch/nco

max_issues_repo_name: hbdch/nco

max_forks_repo_name: hbdch/nco


text:
/*  These functions compute linear preprocessing for
    the UE using LAPACKE and CBLAS modules of
    LAPACK libraries.
    MMSE and MMSE whitening filters are available.
    Functions are using RowMajor storage of the
    matrices, like in conventional C. Traditional
    Fortran functions of LAPACK employ ColumnMajor
    data storage. */

#include<stdio.h>
#include<math.h>
#include<complex.h>
#include <stdlib.h>
#include <cblas.h>
#include <string.h>
#include <linux/version.h>
#if RHEL_RELEASE_CODE >= 1796
    #include <lapacke/lapacke_utils.h>
    #include <lapacke/lapacke.h>
#else
    #include <lapacke_utils.h>
    #include <lapacke.h>
#endif
//#define DEBUG_PREPROC

void transpose(int N, float complex *A, float complex *Result)
{
    // COnputes C := alpha*op(A)*op(B) + beta*C,
    enum CBLAS_TRANSPOSE transa = CblasTrans;
    enum CBLAS_TRANSPOSE transb = CblasNoTrans;
    int rows_opA = N; // number of rows in op(A) and in C
    int col_opB = N; //number of columns of op(B) and in C
    int col_opA = N; //number of columns in op(A) and rows in op(B)
    int col_B; //number of columns in B
    float complex alpha = 1.0 + I * 0;
    int lda  = rows_opA;
    float complex beta = 0.0 + I * 0;
    int ldc = rows_opA;
    int i;
    float complex *B;

    int ldb = col_opB;

    if(transb == CblasNoTrans)
    {
        B = (float complex *)calloc(ldb * col_opB, sizeof(float complex));
        col_B = col_opB;
    }
    else
    {
        B = (float complex *)calloc(ldb * col_opA, sizeof(float complex));
        col_B = col_opA;
    }
    float complex *C = (float complex *)malloc(ldc * col_opB * sizeof(float complex));

    for(i = 0; i < lda * col_B; i += N + 1)
    {
        B[i] = 1.0 + I * 0;
    }

    cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, B, ldb, &beta, C, ldc);

    memcpy(Result, C, N * N * sizeof(float complex));

    free(B);
    free(C);
}

void conjugate_transpose(int N, float complex *A, float complex *Result)
{
    // Computes C := alpha*op(A)*op(B) + beta*C,
    enum CBLAS_TRANSPOSE transa = CblasConjTrans;
    enum CBLAS_TRANSPOSE transb = CblasNoTrans;
    int rows_opA = N; // number of rows in op(A) and in C
    int col_opB = N; //number of columns of op(B) and in C
    int col_opA = N; //number of columns in op(A) and rows in op(B)
    int col_B; //number of columns in B
    float complex alpha = 1.0 + I * 0;
    int lda  = rows_opA;
    float complex beta = 0.0 + I * 0;
    int ldc = rows_opA;
    int i;
    float complex *B;
    int ldb = col_opB;

    if(transb == CblasNoTrans)
    {
        B = (float complex *)calloc(ldb * col_opB, sizeof(float complex));
        col_B = col_opB;
    }
    else
    {
        B = (float complex *)calloc(ldb * col_opA, sizeof(float complex));
        col_B = col_opA;
    }
    float complex *C = (float complex *)malloc(ldc * col_opB * sizeof(float complex));

    for(i = 0; i < lda * col_B; i += N + 1)
    {
        B[i] = 1.0 + I * 0;
    }

    cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, B, ldb, &beta, C, ldc);

    memcpy(Result, C, N * N * sizeof(float complex));

    free(B);
    free(C);
}

void H_hermH_plus_sigma2I(int N, int M, float complex *A, float sigma2, float complex *Result)
{
    //C := alpha*op(A)*op(B) + beta*C,
    enum CBLAS_TRANSPOSE transa = CblasConjTrans;
    enum CBLAS_TRANSPOSE transb = CblasNoTrans;
    int rows_opA = N; // number of rows in op(A) and in C
    int col_opB = N; //number of columns of op(B) and in C
    int col_opA = N; //number of columns in op(A) and rows in op(B)
    int col_C = N; //number of columns in B
    float complex alpha = 1.0 + I * 0;
    int lda  = col_opA;
    float complex beta = 1.0 + I * 0;
    int ldc = col_opA;
    int i;

    float complex *C = (float complex *)calloc(ldc * col_opB, sizeof(float complex));

    for(i = 0; i < lda * col_C; i += N + 1)
    {
        C[i] = sigma2 * (1.0 + I * 0);
    }

    cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, A, lda, A, lda, &beta, C, ldc);

    memcpy(Result, C, N * M * sizeof(float complex));
    free(C);

}

void HH_herm_plus_sigma2I(int M, int N, float complex *A, float sigma2, float complex *Result)
{
    //C := alpha*op(A)*op(B) + beta*C,
    enum CBLAS_TRANSPOSE transa = CblasNoTrans;
    enum CBLAS_TRANSPOSE transb = CblasConjTrans;
    int k = N; //number of columns in op(A) and rows in op(B),k
    float complex alpha = 1.0 + I * 0;
    int lda  = N;
    int ldb  = N;
    int ldc = M;
    int i;

    float complex *C = (float complex *)calloc(M * M, sizeof(float complex));

    for(i = 0; i < M * M; i += M + 1)
    {
        C[i] = 1.0 + I * 0;
    }

    cblas_cgemm(CblasRowMajor, transa, transb, M, M, k, &alpha, A, lda, A, ldb, &sigma2, C, ldc);

    memcpy(Result, C, M * M * sizeof(float complex));
    free(C);

}

void eigen_vectors_values(int N, float complex *A, float complex *Vectors, float *Values_Matrix)
{
    // This function computes ORTHONORMAL eigenvectors and eigenvalues of matrix A,
    // where Values_Matrix is a diagonal matrix of eigenvalues.
    // A=Vectors*Values_Matrix*Vectors'
    char jobz = 'V';
    char uplo = 'U';
    int order_A = N;
    int lda = N;
    int i;
    float *Values = (float *)malloc(sizeof(float) * 1 * N);

    LAPACKE_cheev(LAPACK_ROW_MAJOR, jobz, uplo, order_A, A, lda, Values);

    memcpy(Vectors, A, N * N * sizeof(float complex));

    for(i = 0; i < lda; i += 1)
    {
        Values_Matrix[i * (lda + 1)] = Values[i];
    }

    free(Values);
}

void lin_eq_solver(int N, float complex *A, float complex *B, float complex *Result)
{
    int n = N;
    int lda = N;
    int ldb = N;
    int nrhs = N;

    char transa = 'N';
    int *IPIV = malloc(N * N * sizeof(int));

    // Compute LU-factorization
    LAPACKE_cgetrf(LAPACK_ROW_MAJOR, n, nrhs, A, lda, IPIV);

    // Solve AX=B
    LAPACKE_cgetrs(LAPACK_ROW_MAJOR, transa, n, nrhs, A, lda, IPIV, B, ldb);

    // cgetrs( "N", N, 4, A, lda, IPIV, B, ldb, INFO )

    memcpy(Result, B, N * N * sizeof(float complex));

    free(IPIV);

}

void mutl_matrix_matrix_row_based(float complex *M0, float complex *M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex *Result)
{
    enum CBLAS_TRANSPOSE transa = CblasNoTrans;
    enum CBLAS_TRANSPOSE transb = CblasNoTrans;
    int rows_opA = rows_M0; // number of rows in op(A) and in C
    int col_opB = col_M1; //number of columns of op(B) and in C
    int col_opA = col_M0; //number of columns in op(A) and rows in op(B)
    float complex alpha = 1.0;
    int lda  = col_M0;
    float complex beta = 0.0;
    int ldc = col_M1;
    int ldb = col_M1;

#ifdef DEBUG_PREPROC
    int i = 0;
    printf("rows_M0 %d, col_M0 %d, rows_M1 %d, col_M1 %d\n", rows_M0, col_M0, rows_M1, col_M1);

    for(i = 0; i < rows_M0 * col_M0; ++i)
    {
        printf(" rows_opA = %d, col_opB = %d, W_MMSE[%d] = (%f + i%f)\n", rows_opA, col_opB,  i, creal(M0[i]), cimag(M0[i]));
    }

    for(i = 0; i < rows_M1 * col_M1; ++i)
    {
        printf(" M1[%d] = (%f + i%f)\n", i, creal(M1[i]), cimag(M1[i]));
    }
#endif

    cblas_cgemm(CblasRowMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, M0, lda, M1, ldb, &beta, Result, ldc);

#ifdef DEBUG_PREPROC
    for(i = 0; i < rows_opA * col_opB; ++i)
    {
        printf(" result[%d] = (%f + i%f)\n", i, creal(Result[i]), cimag(Result[i]));
    }
#endif

}
void mutl_matrix_matrix_col_based(float complex *M0, float complex *M1, int rows_M0, int col_M0, int rows_M1, int col_M1, float complex *Result)
{
    enum CBLAS_TRANSPOSE transa = CblasNoTrans;
    enum CBLAS_TRANSPOSE transb = CblasNoTrans;
    int rows_opA = rows_M0; // number of rows in op(A) and in C
    int col_opB = col_M1; //number of columns of op(B) and in C
    int col_opA = col_M0; //number of columns in op(A) and rows in op(B)
    float complex alpha = 1.0;
    int lda  = col_M0;
    float complex beta = 0.0;
    int ldc = rows_M1;
    int ldb = rows_M1;

#ifdef DEBUG_PREPROC
    int i = 0;
    printf("rows_M0 %d, col_M0 %d, rows_M1 %d, col_M1 %d\n", rows_M0, col_M0, rows_M1, col_M1);

    for(i = 0; i < rows_M0 * col_M0; ++i)
    {
        printf(" rows_opA = %d, col_opB = %d, W_MMSE[%d] = (%f + i%f)\n", rows_opA, col_opB,  i, creal(M0[i]), cimag(M0[i]));
    }

    for(i = 0; i < rows_M1 * col_M1; ++i)
    {
        printf(" M1[%d] = (%f + i%f)\n", i, creal(M1[i]), cimag(M1[i]));
    }
#endif

    cblas_cgemm(CblasColMajor, transa, transb, rows_opA, col_opB, col_opA, &alpha, M0, lda, M1, ldb, &beta, Result, ldc);

#ifdef DEBUG_PREPROC
    for(i = 0; i < rows_opA * col_opB; ++i)
    {
        printf(" result[%d] = (%f + i%f)\n", i, creal(Result[i]), cimag(Result[i]));
    }
#endif
}

/*FILTERS */
void compute_MMSE(float complex *H, int order_H, float sigma2, float complex *W_MMSE)
{
    int N = order_H;
    float complex *H_hermH_sigmaI = malloc(N * N * sizeof(float complex));
    float complex *H_herm =  malloc(N * N * sizeof(float complex));

    H_hermH_plus_sigma2I(N, N, H, sigma2, H_hermH_sigmaI);

#ifdef DEBUG_PREPROC
    int i = 0;
    for(i = 0; i < N * N; i++)
    {
        printf(" H_hermH_sigmaI[%d] = (%f + i%f)\n", i, creal(H_hermH_sigmaI[i]), cimag(H_hermH_sigmaI[i]));
    }
#endif

    conjugate_transpose(N, H, H_herm);  //equals H_herm

#ifdef DEBUG_PREPROC
    for(i = 0; i < N * N; i++)
    {
        printf(" H_herm[%d] = (%f + i%f)\n", i, creal(H_herm[i]), cimag(H_herm[i]));
    }
#endif

    lin_eq_solver(N, H_hermH_sigmaI, H_herm, W_MMSE);

#ifdef DEBUG_PREPROC
    for(i = 0; i < N * N; i++)
    {
        printf(" W_MMSE[%d] = (%f + i%f)\n", i, creal(W_MMSE[i]), cimag(W_MMSE[i]));
    }
#endif

    free(H_hermH_sigmaI);
    free(H_herm);
}

float sqrt_float(float x, float sqrt_x)
{
    sqrt_x = (float)(sqrt((double)(x)));
    return sqrt_x;
}