text:
#pragma once
/*
 * (C) Copyright 2020-2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \addtogroup ioda_cxx_attribute
 *
 * @{
 * \file Has_Attributes.h
 * \brief Interfaces for ioda::Has_Attributes and related classes.
 */

#include <gsl/gsl-lite.hpp>
#include <iostream>
#include <memory>
#include <string>
#include <utility>
#include <vector>

#include "ioda/Attributes/Attribute.h"
#include "ioda/Exception.h"
#include "ioda/Misc/Dimensions.h"
#include "ioda/Misc/Eigen_Compat.h"
#include "ioda/Types/Type.h"
#include "ioda/defs.h"

namespace ioda {
class Has_Attributes;
namespace detail {
class Has_Attributes_Backend;
class Has_Attributes_Base;
class Variable_Backend;
class Group_Backend;

/// \brief Describes the functions that can add attributes
/// \ingroup ioda_cxx_attribute
/// \details Using the (Curiously Recurring Template
/// Pattern)[https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern] to implement this
/// since we want to use the same functions for placeholder attribute construction.
/// \see Has_Attributes
/// \see Attribute_Creator
template <class DerivedHasAtts>
class CanAddAttributes {
public:
  /// @name Convenience functions for adding attributes
  /// @{

  /// \brief Create and write an Attribute, for arbitrary dimensions.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be created.
  /// \param data is a gsl::span (a pointer-length pair) that contains the data to be written.
  /// \param dimensions is an initializer list representing the size of the metadata. Each element
  /// is a dimension with a certain size.
  /// \returns Another instance of this Has_Attribute. Used for operation chaining.
  /// \throws jedi::xError if data.size() does not match the number of total elements
  ///   described by dimensions.
  /// \see gsl::span for details of how to make a span.
  template <class DataType>
  DerivedHasAtts add(const std::string& attrname, ::gsl::span<const DataType> data,
                     const ::std::vector<Dimensions_t>& dimensions) {
    auto derivedThis = static_cast<DerivedHasAtts*>(this);
    auto att         = derivedThis->template create<DataType>(attrname, dimensions);
    att.template write<DataType>(data);
    return *derivedThis;
  }

  /// \brief Create and write an Attribute, for arbitrary dimensions.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be created.
  /// \param data is an initializer list that contains the data to be written.
  /// \param dimensions is an initializer list representing the size of the metadata. Each
  ///   element is a dimension with a certain size.
  /// \returns Another instance of this Has_Attribute. Used for operation chaining.
  /// \throws jedi::xError if data.size() does not match the number of total
  ///   elements described by dimensions.
  template <class DataType>
  DerivedHasAtts add(const std::string& attrname, ::std::initializer_list<DataType> data,
                     const ::std::vector<Dimensions_t>& dimensions) {
    auto derivedThis = static_cast<DerivedHasAtts*>(this);
    auto att         = derivedThis->template create<DataType>(attrname, dimensions);
    att.template write<DataType>(data);
    return *derivedThis;
  }

  /// \brief Create and write an Attribute, for a single-dimensional span of 1-D data.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be created.
  /// \param data is a gsl::span (a pointer-length pair) that contains the data to be written.
  ///   The new Attribute will be one-dimensional and the length of the overall span.
  /// \returns Another instance of this Has_Attribute. Used for operation chaining.
  /// \see gsl::span for details of how to make a span.
  /// \see gsl::make_span
  template <class DataType>
  DerivedHasAtts add(const std::string& attrname, ::gsl::span<const DataType> data) {
    return add(attrname, data, {gsl::narrow<Dimensions_t>(data.size())});
  }

  /// \brief Create and write an Attribute, for a 1-D initializer list.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be created.
  /// \param data is an initializer list that contains the data to be written.
  ///   The new Attribute will be one-dimensional and the length of the overall span.
  /// \returns Another instance of this Has_Attribute. Used for operation chaining.
  template <class DataType>
  DerivedHasAtts add(const std::string& attrname, ::std::initializer_list<DataType> data) {
    return add(attrname, data, {gsl::narrow<Dimensions_t>(data.size())});
  }

  /// \brief Create and write a single datum of an Attribute.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be created.
  /// \param data is the data to be written. The new Attribute will be zero-dimensional and will
  /// contain only this datum.
  /// \note Even 0-dimensional data have a type, which may be a compound
  /// array (i.e. a single string of variable length).
  /// \returns Another instance of this Has_Attribute. Used for operation chaining.
  template <class DataType>
  DerivedHasAtts add(const std::string& attrname, const DataType& data) {
    return add(attrname, gsl::make_span(&data, 1), {1});
  }

  template <class EigenClass>
  DerivedHasAtts addWithEigenRegular(const std::string& attrname, const EigenClass& data,
                                     bool is2D = true) {
#if 1  //__has_include("Eigen/Dense")
    typedef typename EigenClass::Scalar ScalarType;
    // If d is already in Row Major form, then this is optimized out.
    Eigen::Array<ScalarType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> dout;
    dout.resize(data.rows(), data.cols());
    dout               = data;
    const auto& dconst = dout;  // To make some compilers happy.
    auto sp            = gsl::make_span(dconst.data(), static_cast<int>(data.rows() * data.cols()));

    if (is2D)
      return add(attrname, sp,
                 {gsl::narrow<Dimensions_t>(data.rows()), gsl::narrow<Dimensions_t>(data.cols())});
    else
      return add(attrname, sp);
#else
    static_assert(false, "The Eigen headers cannot be found, so this function cannot be used.");
#endif
  }

  template <class EigenClass>
  DerivedHasAtts addWithEigenTensor(const std::string& attrname, const EigenClass& data) {
#if 1  //__has_include("unsupported/Eigen/CXX11/Tensor")
    typedef typename EigenClass::Scalar ScalarType;
    ioda::Dimensions dims = detail::EigenCompat::getTensorDimensions(data);

    auto derived_this = static_cast<DerivedHasAtts*>(this);
    Attribute att     = derived_this->template create<ScalarType>(attrname, dims.dimsCur);

    att.writeWithEigenTensor(data);
    return *derived_this;
#else
    static_assert(
      false, "The Eigen unsupported/ headers cannot be found, so this function cannot be used.");
#endif
  }

  /// @}
};

/// \brief Describes the functions that can read attributes
/// \ingroup ioda_cxx_attribute
/// \details Uses the (CRTP)[https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern]
///   to implement this as we use the same functions for placeholder attribute construction.
/// \see Has_Attributes
template <class DerivedHasAtts>
class CanReadAttributes {
protected:
  CanReadAttributes() {}

public:
  virtual ~CanReadAttributes() {}

  /// @name Convenience functions for reading attributes
  /// @{

  /// \brief Open and read an Attribute, for expected dimensions.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be read.
  /// \param data is a pointer-size pair to the data buffer that is filled with the metadata's
  ///   contents. It should be pre-sized to accomodate all of the matadata. See
  ///   getDimensions().numElements. Data will be filled in row-major order.
  /// \throws jedi::xError on a size mismatch between Attribute dimensions and data.size().
  template <class DataType>
  const DerivedHasAtts read(const std::string& attrname, gsl::span<DataType> data) const {
    // Attribute att = this->open(attrname);
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    att.read(data);
    return *(static_cast<const DerivedHasAtts*>(this));
  }

  /// \brief Open and read an Attribute, with unknown dimensions.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be read.
  /// \param data is a vector acting as a data buffer that is filled with the metadata's
  ///   contents. It gets resized as needed. data will be filled in row-major order.
  template <class DataType>
  const DerivedHasAtts read(const std::string& attrname, std::vector<DataType>& data) const {
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    att.read(data);
    return *(static_cast<const DerivedHasAtts*>(this));
  }

  /// \brief Open and read an Attribute, with unknown dimensions.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be read.
  /// \param data is a valarray acting as a data buffer that is filled with the metadata's
  ///   contents. It gets resized as needed. data will be filled in row-major order.
  template <class DataType>
  const DerivedHasAtts read(const std::string& attrname, std::valarray<DataType>& data) const {
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    att.read(data);
    return *(static_cast<const DerivedHasAtts*>(this));
  }

  /// \brief Read a datum of an Attribute.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be read.
  /// \param data is a datum of type DataType.
  /// \throws jedi::xError if the underlying data have size > 1.
  template <class DataType>
  const DerivedHasAtts read(const std::string& attrname, DataType& data) const {
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    att.read<DataType>(data);
    return *(static_cast<const DerivedHasAtts*>(this));
  }

  /// \brief Read a datum of an Attribute.
  /// \tparam DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute to be read.
  /// \returns A datum of type DataType.
  /// \throws jedi::xError if the underlying data have size > 1.
  template <class DataType>
  DataType read(const std::string& attrname) const {
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    return att.read<DataType>();
  }

  template <class EigenClass, bool Resize = detail::EigenCompat::CanResize<EigenClass>::value>
  DerivedHasAtts readWithEigenRegular(const std::string& attrname, EigenClass& data) {
#if 1  // __has_include("Eigen/Dense")
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    att.readWithEigenRegular<EigenClass, Resize>(data);
    return *(static_cast<const DerivedHasAtts*>(this));
#else
    static_assert(false, "The Eigen headers cannot be found, so this function cannot be used.")
#endif
  }

  template <class EigenClass>
  DerivedHasAtts readWithEigenTensor(const std::string& attrname, EigenClass& data) {
#if 1  //__has_include("unsupported/Eigen/CXX11/Tensor")
    Attribute att = static_cast<const DerivedHasAtts*>(this)->open(attrname);
    att.readWithEigenTensor<EigenClass>(data);
    return *(static_cast<const DerivedHasAtts*>(this));
#else
    static_assert(
      false, "The Eigen unsupported/ headers cannot be found, so this function cannot be used.");
#endif
  }

  /// @}
};

/// \ingroup ioda_cxx_attribute
class IODA_DL Has_Attributes_Base {
private:
  /// Using an opaque object to implement the backend.
  std::shared_ptr<Has_Attributes_Backend> backend_;

  // Backend classes do occasionally need access to the Attribute_Backend object.
  friend class Group_Backend;
  friend class Variable_Backend;

protected:
  Has_Attributes_Base(std::shared_ptr<Has_Attributes_Backend>);

public:
  virtual ~Has_Attributes_Base();

  /// Query the backend and get the type provider.
  virtual detail::Type_Provider* getTypeProvider() const;

  /// @name General Functions
  /// @{
  ///

  /// List all attributes
  /// \returns an unordered vector of attribute names for an object.
  virtual std::vector<std::string> list() const;
  /// List all attributes
  /// \returns an unordered vector of attribute names for an object.
  inline std::vector<std::string> operator()() const { return list(); }

  /// \brief Does an Attribute with the specified name exist?
  /// \param attname is the name of the Attribute that we are looking for.
  /// \returns true if it exists.
  /// \returns false otherwise.
  virtual bool exists(const std::string& attname) const;
  /// \brief Delete an Attribute with the specified name.
  /// \param attname is the name of the Attribute that we are deleting.
  /// \throws jedi::xError if no such attribute exists.
  virtual void remove(const std::string& attname);
  /// \brief Open an Attribute by name
  /// \param name is the name of the Attribute to be opened.
  /// \returns An instance of an Attribute that can be queried (with getDimensions()) and read.
  virtual Attribute open(const std::string& name) const;
  /// \brief Open Attribute by name
  /// \param name is the name of the Attribute to be opened.
  /// \returns An instance of Attribute that can be queried (with getDimensions()) and read.
  inline Attribute operator[](const std::string& name) const { return open(name); }

  /// \brief Open all attributes in an object.
  /// \details This is a collective call, optimized for performance.
  /// \return A sequence of (name, Attribute) pairs.
  virtual std::vector<std::pair<std::string, Attribute>> openAll() const;

  /// \brief Create an Attribute without setting its data.
  /// \param attrname is the name of the Attribute.
  /// \param dimensions is a vector representing the size of the metadata.
  ///   Each element of the vector is a dimension with a certain size.
  /// \param in_memory_datatype is the runtime description of the Attribute's data type.
  /// \returns An instance of Attribute that can be written to.
  virtual Attribute create(const std::string& attrname, const Type& in_memory_dataType,
                           const std::vector<Dimensions_t>& dimensions = {1});

  /// Python compatability function
  inline Attribute _create_py(const std::string& attrname, BasicTypes dataType,
                              const std::vector<Dimensions_t>& dimensions = {1}) {
    return create(attrname, ::ioda::Type(dataType, getTypeProvider()), dimensions);
  }

  /// \brief Create an Attribute without setting its data.
  /// \tparam DataType is the type of the data. I.e. float, int, uint16_t, std::string, etc.
  /// \param attrname is the name of the Attribute.
  /// \param dimensions is a vector representing the size of the metadata. Each element of the
  ///   vector is a dimension with a certain size.
  /// \returns An instance of Attribute that can be written to.
  template <class DataType, class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Attribute create(const std::string& attrname,
                   const std::vector<Dimensions_t>& dimensions = {1}) {
    try {
      Type in_memory_dataType = TypeWrapper::GetType(getTypeProvider());
      auto att                = create(attrname, in_memory_dataType, dimensions);
      return att;
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// \brief Rename an Attribute
  /// \param oldName is the name of the Attribute to be changed.
  /// \param newName is the new name of the Attribute.
  /// \throws jedi::xError if oldName is not found.
  /// \throws jedi::xError if newName already exists.
  virtual void rename(const std::string& oldName, const std::string& newName);

  /// @}
};

class IODA_DL Has_Attributes_Backend : public Has_Attributes_Base {
protected:
  Has_Attributes_Backend();

public:
  virtual ~Has_Attributes_Backend();

  /// \brief Default implementation of Has_Attributes_Base::openAll.
  /// \see Has_Attributes_Base::openAll.
  std::vector<std::pair<std::string, Attribute>> openAll() const override;
};
}  // namespace detail

/** \brief This class exists inside of ioda::Group or ioda::Variable and provides
 *         the interface to manipulating Attributes.
 * \ingroup ioda_cxx_attribute
 *
 * \note It should only be constructed inside of a Group or Variable.
 *       It has no meaning elsewhere.
 * \see Attribute for the class that represents individual attributes.
 * \throws jedi::xError on all exceptions.
 **/
class IODA_DL Has_Attributes : public detail::CanAddAttributes<Has_Attributes>,
                               public detail::CanReadAttributes<Has_Attributes>,
                               public detail::Has_Attributes_Base {
public:
  Has_Attributes();
  Has_Attributes(std::shared_ptr<detail::Has_Attributes_Backend>);
  virtual ~Has_Attributes();
};

}  // namespace ioda

/// @} // End Doxygen block

max_stars_repo_name: Gibies/ioda

max_issues_repo_name: Gibies/ioda

max_forks_repo_name: Gibies/ioda


text:
/* rng/zuf.c
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* It is crucial that m == n-273 mod 607 at all times;
   For speed of execution, however, this is never enforced.
   Instead is is set in the initializer: note 607-273=334
   Note also that the state.u[607] is not initialized */

static inline unsigned long int zuf_get (void *vstate);
static double zuf_get_double (void *vstate);
static void zuf_set (void *state, unsigned long int s);

static const unsigned long int zuf_randmax = 16777216;	/* 2^24 */

typedef struct
  {
    int n;
    unsigned long int u[607];
  }
zuf_state_t;

/* The zufall package was implemented with float's, which is to say 24
   bits of precision.  Since I'm using long's instead, my RANDMAX
   reflects this. */

static inline unsigned long int
zuf_get (void *vstate)
{
  zuf_state_t *state = (zuf_state_t *) vstate;
  const int n = state->n;
  const int m = (n - 273 + 607) % 607;
  unsigned long int t = state->u[n] + state->u[m];

  while (t > zuf_randmax)
    t -= zuf_randmax;

  state->u[n] = t;

  if (n == 606)
    {
      state->n = 0;
    }
  else
    {
      state->n = n + 1;
    }

  return t;
}

static double
zuf_get_double (void *vstate)
{
  return zuf_get (vstate) / 16777216.0 ;
}

static void
zuf_set (void *vstate, unsigned long int s)
{
  /* A very elaborate seeding procedure is provided with the
     zufall package; this is virtually a copy of that procedure */

  /* Initialized data */

  long int kl = 9373;
  long int ij = 1802;

  /* Local variables */
  long int i, j, k, l, m;
  double x, y;
  long int ii, jj;

  zuf_state_t *state = (zuf_state_t *) vstate;

  state->n = 0;

/*  generates initial seed buffer by linear congruential */
/*  method. Taken from Marsaglia, FSU report FSU-SCRI-87-50 */
/*  variable seed should be 0 < seed <31328 */

  if (s == 0)
    s = 1802;	/* default seed is 1802 */

  ij = s;

  i = ij / 177 % 177 + 2;
  j = ij % 177 + 2;
  k = kl / 169 % 178 + 1;
  l = kl % 169;
  for (ii = 0; ii < 607; ++ii)
    {
      x = 0.0;
      y = 0.5;
      /* 24 bits?? */
      for (jj = 1; jj <= 24; ++jj)
	{
	  m = i * j % 179 * k % 179;
	  i = j;
	  j = k;
	  k = m;
	  l = (l * 53 + 1) % 169;
	  if (l * m % 64 >= 32)
	    {
	      x += y;
	    }
	  y *= 0.5;
	}
      state->u[ii] = (unsigned long int) (x * zuf_randmax);
    }
}

static const gsl_rng_type zuf_type =
{"zuf",				/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (zuf_state_t),
 &zuf_set,
 &zuf_get,
 &zuf_get_double};

const gsl_rng_type *gsl_rng_zuf = &zuf_type;

ext: c


text:
#pragma once

#include <gsl/gsl>

#include "exception.h"
#include <fmt/format.h>
#include "mathutil.h"
#include "stringutil.h"
#include "version.h"
#include "stopwatch.h"
#include "logging.h"

ext: h

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

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

#include "norm.h"
#include "constants.h"

int readFileVector( gsl_vector *vector, char *filePath, char *fileName )                //read file vals into vector
{
	char tmpFilepath[120] = "\0";
	char tmpReading[20];
	char *tmp2;
	double reading;                                 //WARNING this is similar to global var $readings, consider renaming

	strcat( tmpFilepath, filePath );                //TODO replace this and the first instance of strcat into strcpy, that should work
	strcat( tmpFilepath, fileName );        //combines the filepath and name into one searchable /path/file
	FILE *file = fopen( tmpFilepath, "rt");
	printf("%s\n", tmpFilepath);                    //TODO DELETE THIS
	if ( file == NULL ) {                   //check if file exists
		fprintf( stderr, "Error, unable to locate the Mu data file\n");
		exit(100);
	}

	for ( int i = 0; i < skip; i++)         //move past header of read files
	{
		fgets( tmpReading, 19, file );
	}

	for ( int i = 0; i < readings; i++ )
	{
		fgets( tmpReading, 19, file );

		tmp2 = strrchr( tmpReading, '\t' );
		tmp2 = tmp2 + 1;
		reading = atof(tmp2);
		gsl_vector_set( vector, i, reading);
	}
	fclose(file);
	return 0;
}

int readFileMatrix( gsl_matrix *matrix, char *filePath, int argc, char *argv[] )
{
	char *tmp2;
	double reading;
	for ( int j = 2; j < argc; j++ )
	{
		gsl_vector_view column;
		char tmpFilepath[120] = "\0";                           //TODO replace this and the first instance of strcat into strcpy, that should work
		char tmpReading[20];
		strcat( tmpFilepath, filePath );
		strcat( tmpFilepath, argv[j] );         //combines the filepath and name into one searchable /path/file
		FILE *file = fopen( tmpFilepath, "rt");
		printf("%s\n", tmpFilepath);                            //TODO DELETE THIS
		if ( file == NULL ) {                   //check if file exists
			fprintf( stderr, "Error, unable to locate a reference data file\n");
			exit(100);
		}
		for ( int i = 0; i < skip; i++)         //move past header of read files
		{
			fgets( tmpReading, 20, file);
		}
		for ( int i = 0; i < readings; i++ ) //fill the ith element of the jth column
		{
			fgets( tmpReading, 19, file );

			tmp2 = strrchr( tmpReading, '\t' );
			tmp2 = tmp2 + 1;
			reading = atof(tmp2);
			gsl_matrix_set( matrix, i, j-2, reading );
		}
		fclose(file);
		column = gsl_matrix_column( matrix, j-2 );                              //NORMALISE EACH REFERENCE VECOTR.
		normalizeVector( &column.vector );
	}

	return 0;
}

max_stars_repo_name: jakeinater/spectralILU

max_issues_repo_name: jakeinater/spectralILU

max_forks_repo_name: jakeinater/spectralILU


text:
/* specfunc/hyperg_2F1.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
 * Copyright (C) 2009 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_hyperg.h>

#include "error.h"

#define locEPS (1000.0*GSL_DBL_EPSILON)

/* Assumes c != negative integer.
 */
static int
hyperg_2F1_series(const double a, const double b, const double c,
                  const double x,
                  gsl_sf_result * result
                  )
{
  double sum_pos = 1.0;
  double sum_neg = 0.0;
  double del_pos = 1.0;
  double del_neg = 0.0;
  double del = 1.0;
  double k = 0.0;
  int i = 0;

  if(fabs(c) < GSL_DBL_EPSILON) {
    result->val = 0.0; /* FIXME: ?? */
    result->err = 1.0;
    GSL_ERROR ("error", GSL_EDOM);
  }

  do {
    if(++i > 30000) {
      result->val  = sum_pos - sum_neg;
      result->err  = del_pos + del_neg;
      result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
      result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
      GSL_ERROR ("error", GSL_EMAXITER);
    }
    del *= (a+k)*(b+k) * x / ((c+k) * (k+1.0));  /* Gauss series */

    if(del > 0.0) {
      del_pos  =  del;
      sum_pos +=  del;
    }
    else if(del == 0.0) {
      /* Exact termination (a or b was a negative integer).
       */
      del_pos = 0.0;
      del_neg = 0.0;
      break;
    }
    else {
      del_neg  = -del;
      sum_neg -=  del;
    }

    k += 1.0;
  } while(fabs((del_pos + del_neg)/(sum_pos-sum_neg)) > GSL_DBL_EPSILON);

  result->val  = sum_pos - sum_neg;
  result->err  = del_pos + del_neg;
  result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
  result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);

  return GSL_SUCCESS;
}

/* a = aR + i aI, b = aR - i aI */
static
int
hyperg_2F1_conj_series(const double aR, const double aI, const double c,
                       double x,
                       gsl_sf_result * result)
{
  if(c == 0.0) {
    result->val = 0.0; /* FIXME: should be Inf */
    result->err = 0.0;
    GSL_ERROR ("error", GSL_EDOM);
  }
  else {
    double sum_pos = 1.0;
    double sum_neg = 0.0;
    double del_pos = 1.0;
    double del_neg = 0.0;
    double del = 1.0;
    double k = 0.0;
    do {
      del *= ((aR+k)*(aR+k) + aI*aI)/((k+1.0)*(c+k)) * x;

      if(del >= 0.0) {
        del_pos  =  del;
        sum_pos +=  del;
      }
      else {
        del_neg  = -del;
        sum_neg -=  del;
      }

      if(k > 30000) {
        result->val  = sum_pos - sum_neg;
        result->err  = del_pos + del_neg;
        result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
        result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
        GSL_ERROR ("error", GSL_EMAXITER);
      }

      k += 1.0;
    } while(fabs((del_pos + del_neg)/(sum_pos - sum_neg)) > GSL_DBL_EPSILON);

    result->val  = sum_pos - sum_neg;
    result->err  = del_pos + del_neg;
    result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
    result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

/* Luke's rational approximation. The most accesible
 * discussion is in [Kolbig, CPC 23, 51 (1981)].
 * The convergence is supposedly guaranteed for x < 0.
 * You have to read Luke's books to see this and other
 * results. Unfortunately, the stability is not so
 * clear to me, although it seems very efficient when
 * it works.
 */
static
int
hyperg_2F1_luke(const double a, const double b, const double c,
                const double xin,
                gsl_sf_result * result)
{
  int stat_iter;
  const double RECUR_BIG = 1.0e+50;
  const int nmax = 20000;
  int n = 3;
  const double x  = -xin;
  const double x3 = x*x*x;
  const double t0 = a*b/c;
  const double t1 = (a+1.0)*(b+1.0)/(2.0*c);
  const double t2 = (a+2.0)*(b+2.0)/(2.0*(c+1.0));
  double F = 1.0;
  double prec;

  double Bnm3 = 1.0;                                  /* B0 */
  double Bnm2 = 1.0 + t1 * x;                         /* B1 */
  double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */

  double Anm3 = 1.0;                                                      /* A0 */
  double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
  double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */

  while(1) {
    double npam1 = n + a - 1;
    double npbm1 = n + b - 1;
    double npcm1 = n + c - 1;
    double npam2 = n + a - 2;
    double npbm2 = n + b - 2;
    double npcm2 = n + c - 2;
    double tnm1  = 2*n - 1;
    double tnm3  = 2*n - 3;
    double tnm5  = 2*n - 5;
    double n2 = n*n;
    double F1 =  (3.0*n2 + (a+b-6)*n + 2 - a*b - 2*(a+b)) / (2*tnm3*npcm1);
    double F2 = -(3.0*n2 - (a+b+6)*n + 2 - a*b)*npam1*npbm1/(4*tnm1*tnm3*npcm2*npcm1);
    double F3 = (npam2*npam1*npbm2*npbm1*(n-a-2)*(n-b-2)) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
    double E  = -npam1*npbm1*(n-c-1) / (2*tnm3*npcm2*npcm1);

    double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
    double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
    double r = An/Bn;

    prec = fabs((F - r)/F);
    F = r;

    if(prec < GSL_DBL_EPSILON || n > nmax) break;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An   /= RECUR_BIG;
      Bn   /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
      Anm3 /= RECUR_BIG;
      Bnm3 /= RECUR_BIG;
    }
    else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
      An   *= RECUR_BIG;
      Bn   *= RECUR_BIG;
      Anm1 *= RECUR_BIG;
      Bnm1 *= RECUR_BIG;
      Anm2 *= RECUR_BIG;
      Bnm2 *= RECUR_BIG;
      Anm3 *= RECUR_BIG;
      Bnm3 *= RECUR_BIG;
    }

    n++;
    Bnm3 = Bnm2;
    Bnm2 = Bnm1;
    Bnm1 = Bn;
    Anm3 = Anm2;
    Anm2 = Anm1;
    Anm1 = An;
  }

  result->val  = F;
  result->err  = 2.0 * fabs(prec * F);
  result->err += 2.0 * GSL_DBL_EPSILON * (n+1.0) * fabs(F);

  /* FIXME: just a hack: there's a lot of shit going on here */
  result->err *= 8.0 * (fabs(a) + fabs(b) + 1.0);

  stat_iter = (n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );

  return stat_iter;
}

/* Luke's rational approximation for the
 * case a = aR + i aI, b = aR - i aI.
 */
static
int
hyperg_2F1_conj_luke(const double aR, const double aI, const double c,
                     const double xin,
                     gsl_sf_result * result)
{
  int stat_iter;
  const double RECUR_BIG = 1.0e+50;
  const int nmax = 10000;
  int n = 3;
  const double x = -xin;
  const double x3 = x*x*x;
  const double atimesb = aR*aR + aI*aI;
  const double apb     = 2.0*aR;
  const double t0 = atimesb/c;
  const double t1 = (atimesb +     apb + 1.0)/(2.0*c);
  const double t2 = (atimesb + 2.0*apb + 4.0)/(2.0*(c+1.0));
  double F = 1.0;
  double prec;

  double Bnm3 = 1.0;                                  /* B0 */
  double Bnm2 = 1.0 + t1 * x;                         /* B1 */
  double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */

  double Anm3 = 1.0;                                                      /* A0 */
  double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
  double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */

  while(1) {
    double nm1 = n - 1;
    double nm2 = n - 2;
    double npam1_npbm1 = atimesb + nm1*apb + nm1*nm1;
    double npam2_npbm2 = atimesb + nm2*apb + nm2*nm2;
    double npcm1 = nm1 + c;
    double npcm2 = nm2 + c;
    double tnm1  = 2*n - 1;
    double tnm3  = 2*n - 3;
    double tnm5  = 2*n - 5;
    double n2 = n*n;
    double F1 =  (3.0*n2 + (apb-6)*n + 2 - atimesb - 2*apb) / (2*tnm3*npcm1);
    double F2 = -(3.0*n2 - (apb+6)*n + 2 - atimesb)*npam1_npbm1/(4*tnm1*tnm3*npcm2*npcm1);
    double F3 = (npam2_npbm2*npam1_npbm1*(nm2*nm2 - nm2*apb + atimesb)) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
    double E  = -npam1_npbm1*(n-c-1) / (2*tnm3*npcm2*npcm1);

    double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
    double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
    double r = An/Bn;

    prec = fabs(F - r)/fabs(F);
    F = r;

    if(prec < GSL_DBL_EPSILON || n > nmax) break;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An   /= RECUR_BIG;
      Bn   /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
      Anm3 /= RECUR_BIG;
      Bnm3 /= RECUR_BIG;
    }
    else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
      An   *= RECUR_BIG;
      Bn   *= RECUR_BIG;
      Anm1 *= RECUR_BIG;
      Bnm1 *= RECUR_BIG;
      Anm2 *= RECUR_BIG;
      Bnm2 *= RECUR_BIG;
      Anm3 *= RECUR_BIG;
      Bnm3 *= RECUR_BIG;
    }

    n++;
    Bnm3 = Bnm2;
    Bnm2 = Bnm1;
    Bnm1 = Bn;
    Anm3 = Anm2;
    Anm2 = Anm1;
    Anm1 = An;
  }

  result->val  = F;
  result->err  = 2.0 * fabs(prec * F);
  result->err += 2.0 * GSL_DBL_EPSILON * (n+1.0) * fabs(F);

  /* FIXME: see above */
  result->err *= 8.0 * (fabs(aR) + fabs(aI) + 1.0);

  stat_iter = (n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );

  return stat_iter;
}

/* Do the reflection described in [Moshier, p. 334].
 * Assumes a,b,c != neg integer.
 */
static
int
hyperg_2F1_reflect(const double a, const double b, const double c,
                   const double x, gsl_sf_result * result)
{
  const double d = c - a - b;
  const int intd  = floor(d+0.5);
  const int d_integer = ( fabs(d - intd) < locEPS );

  if(d_integer) {
    const double ln_omx = log(1.0 - x);
    const double ad = fabs(d);
    int stat_F2 = GSL_SUCCESS;
    double sgn_2;
    gsl_sf_result F1;
    gsl_sf_result F2;
    double d1, d2;
    gsl_sf_result lng_c;
    gsl_sf_result lng_ad2;
    gsl_sf_result lng_bd2;
    int stat_c;
    int stat_ad2;
    int stat_bd2;

    if(d >= 0.0) {
      d1 = d;
      d2 = 0.0;
    }
    else {
      d1 = 0.0;
      d2 = d;
    }

    stat_ad2 = gsl_sf_lngamma_e(a+d2, &lng_ad2);
    stat_bd2 = gsl_sf_lngamma_e(b+d2, &lng_bd2);
    stat_c   = gsl_sf_lngamma_e(c,    &lng_c);

    /* Evaluate F1.
     */
    if(ad < GSL_DBL_EPSILON) {
      /* d = 0 */
      F1.val = 0.0;
      F1.err = 0.0;
    }
    else {
      gsl_sf_result lng_ad;
      gsl_sf_result lng_ad1;
      gsl_sf_result lng_bd1;
      int stat_ad  = gsl_sf_lngamma_e(ad,   &lng_ad);
      int stat_ad1 = gsl_sf_lngamma_e(a+d1, &lng_ad1);
      int stat_bd1 = gsl_sf_lngamma_e(b+d1, &lng_bd1);

      if(stat_ad1 == GSL_SUCCESS && stat_bd1 == GSL_SUCCESS && stat_ad == GSL_SUCCESS) {
        /* Gamma functions in the denominator are ok.
         * Proceed with evaluation.
         */
        int i;
        double sum1 = 1.0;
        double term = 1.0;
        double ln_pre1_val = lng_ad.val + lng_c.val + d2*ln_omx - lng_ad1.val - lng_bd1.val;
        double ln_pre1_err = lng_ad.err + lng_c.err + lng_ad1.err + lng_bd1.err + GSL_DBL_EPSILON * fabs(ln_pre1_val);
        int stat_e;

        /* Do F1 sum.
         */
        for(i=1; i<ad; i++) {
          int j = i-1;
          term *= (a + d2 + j) * (b + d2 + j) / (1.0 + d2 + j) / i * (1.0-x);
          sum1 += term;
        }

        stat_e = gsl_sf_exp_mult_err_e(ln_pre1_val, ln_pre1_err,
                                       sum1, GSL_DBL_EPSILON*fabs(sum1),
                                       &F1);
        if(stat_e == GSL_EOVRFLW) {
          OVERFLOW_ERROR(result);
        }
      }
      else {
        /* Gamma functions in the denominator were not ok.
         * So the F1 term is zero.
         */
        F1.val = 0.0;
        F1.err = 0.0;
      }
    } /* end F1 evaluation */

    /* Evaluate F2.
     */
    if(stat_ad2 == GSL_SUCCESS && stat_bd2 == GSL_SUCCESS) {
      /* Gamma functions in the denominator are ok.
       * Proceed with evaluation.
       */
      const int maxiter = 2000;
      double psi_1 = -M_EULER;
      gsl_sf_result psi_1pd;
      gsl_sf_result psi_apd1;
      gsl_sf_result psi_bpd1;
      int stat_1pd  = gsl_sf_psi_e(1.0 + ad, &psi_1pd);
      int stat_apd1 = gsl_sf_psi_e(a + d1,   &psi_apd1);
      int stat_bpd1 = gsl_sf_psi_e(b + d1,   &psi_bpd1);
      int stat_dall = GSL_ERROR_SELECT_3(stat_1pd, stat_apd1, stat_bpd1);

      double psi_val = psi_1 + psi_1pd.val - psi_apd1.val - psi_bpd1.val - ln_omx;
      double psi_err = psi_1pd.err + psi_apd1.err + psi_bpd1.err + GSL_DBL_EPSILON*fabs(psi_val);
      double fact = 1.0;
      double sum2_val = psi_val;
      double sum2_err = psi_err;
      double ln_pre2_val = lng_c.val + d1*ln_omx - lng_ad2.val - lng_bd2.val;
      double ln_pre2_err = lng_c.err + lng_ad2.err + lng_bd2.err + GSL_DBL_EPSILON * fabs(ln_pre2_val);
      int stat_e;

      int j;

      /* Do F2 sum.
       */
      for(j=1; j<maxiter; j++) {
        /* values for psi functions use recurrence; Abramowitz+Stegun 6.3.5 */
        double term1 = 1.0/(double)j  + 1.0/(ad+j);
        double term2 = 1.0/(a+d1+j-1.0) + 1.0/(b+d1+j-1.0);
        double delta = 0.0;
        psi_val += term1 - term2;
        psi_err += GSL_DBL_EPSILON * (fabs(term1) + fabs(term2));
        fact *= (a+d1+j-1.0)*(b+d1+j-1.0)/((ad+j)*j) * (1.0-x);
        delta = fact * psi_val;
        sum2_val += delta;
        sum2_err += fabs(fact * psi_err) + GSL_DBL_EPSILON*fabs(delta);
        if(fabs(delta) < GSL_DBL_EPSILON * fabs(sum2_val)) break;
      }

      if(j == maxiter) stat_F2 = GSL_EMAXITER;

      if(sum2_val == 0.0) {
        F2.val = 0.0;
        F2.err = 0.0;
      }
      else {
        stat_e = gsl_sf_exp_mult_err_e(ln_pre2_val, ln_pre2_err,
                                       sum2_val, sum2_err,
                                       &F2);
        if(stat_e == GSL_EOVRFLW) {
          result->val = 0.0;
          result->err = 0.0;
          GSL_ERROR ("error", GSL_EOVRFLW);
        }
      }
      stat_F2 = GSL_ERROR_SELECT_2(stat_F2, stat_dall);
    }
    else {
      /* Gamma functions in the denominator not ok.
       * So the F2 term is zero.
       */
      F2.val = 0.0;
      F2.err = 0.0;
    } /* end F2 evaluation */

    sgn_2 = ( GSL_IS_ODD(intd) ? -1.0 : 1.0 );
    result->val  = F1.val + sgn_2 * F2.val;
    result->err  = F1.err + F2. err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(F1.val) + fabs(F2.val));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_F2;
  }
  else {
    /* d not an integer */

    gsl_sf_result pre1, pre2;
    double sgn1, sgn2;
    gsl_sf_result F1, F2;
    int status_F1, status_F2;

    /* These gamma functions appear in the denominator, so we
     * catch their harmless domain errors and set the terms to zero.
     */
    gsl_sf_result ln_g1ca,  ln_g1cb,  ln_g2a,  ln_g2b;
    double sgn_g1ca, sgn_g1cb, sgn_g2a, sgn_g2b;
    int stat_1ca = gsl_sf_lngamma_sgn_e(c-a, &ln_g1ca, &sgn_g1ca);
    int stat_1cb = gsl_sf_lngamma_sgn_e(c-b, &ln_g1cb, &sgn_g1cb);
    int stat_2a  = gsl_sf_lngamma_sgn_e(a, &ln_g2a, &sgn_g2a);
    int stat_2b  = gsl_sf_lngamma_sgn_e(b, &ln_g2b, &sgn_g2b);
    int ok1 = (stat_1ca == GSL_SUCCESS && stat_1cb == GSL_SUCCESS);
    int ok2 = (stat_2a  == GSL_SUCCESS && stat_2b  == GSL_SUCCESS);

    gsl_sf_result ln_gc,  ln_gd,  ln_gmd;
    double sgn_gc, sgn_gd, sgn_gmd;
    gsl_sf_lngamma_sgn_e( c, &ln_gc,  &sgn_gc);
    gsl_sf_lngamma_sgn_e( d, &ln_gd,  &sgn_gd);
    gsl_sf_lngamma_sgn_e(-d, &ln_gmd, &sgn_gmd);

    sgn1 = sgn_gc * sgn_gd  * sgn_g1ca * sgn_g1cb;
    sgn2 = sgn_gc * sgn_gmd * sgn_g2a  * sgn_g2b;

    if(ok1 && ok2) {
      double ln_pre1_val = ln_gc.val + ln_gd.val  - ln_g1ca.val - ln_g1cb.val;
      double ln_pre2_val = ln_gc.val + ln_gmd.val - ln_g2a.val  - ln_g2b.val + d*log(1.0-x);
      double ln_pre1_err = ln_gc.err + ln_gd.err + ln_g1ca.err + ln_g1cb.err;
      double ln_pre2_err = ln_gc.err + ln_gmd.err + ln_g2a.err  + ln_g2b.err;
      if(ln_pre1_val < GSL_LOG_DBL_MAX && ln_pre2_val < GSL_LOG_DBL_MAX) {
        gsl_sf_exp_err_e(ln_pre1_val, ln_pre1_err, &pre1);
        gsl_sf_exp_err_e(ln_pre2_val, ln_pre2_err, &pre2);
        pre1.val *= sgn1;
        pre2.val *= sgn2;
      }
      else {
        OVERFLOW_ERROR(result);
      }
    }
    else if(ok1 && !ok2) {
      double ln_pre1_val = ln_gc.val + ln_gd.val - ln_g1ca.val - ln_g1cb.val;
      double ln_pre1_err = ln_gc.err + ln_gd.err + ln_g1ca.err + ln_g1cb.err;
      if(ln_pre1_val < GSL_LOG_DBL_MAX) {
        gsl_sf_exp_err_e(ln_pre1_val, ln_pre1_err, &pre1);
        pre1.val *= sgn1;
        pre2.val = 0.0;
        pre2.err = 0.0;
      }
      else {
        OVERFLOW_ERROR(result);
      }
    }
    else if(!ok1 && ok2) {
      double ln_pre2_val = ln_gc.val + ln_gmd.val - ln_g2a.val - ln_g2b.val + d*log(1.0-x);
      double ln_pre2_err = ln_gc.err + ln_gmd.err + ln_g2a.err + ln_g2b.err;
      if(ln_pre2_val < GSL_LOG_DBL_MAX) {
        pre1.val = 0.0;
        pre1.err = 0.0;
        gsl_sf_exp_err_e(ln_pre2_val, ln_pre2_err, &pre2);
        pre2.val *= sgn2;
      }
      else {
        OVERFLOW_ERROR(result);
      }
    }
    else {
      pre1.val = 0.0;
      pre2.val = 0.0;
      UNDERFLOW_ERROR(result);
    }

    status_F1 = hyperg_2F1_series(  a,   b, 1.0-d, 1.0-x, &F1);
    status_F2 = hyperg_2F1_series(c-a, c-b, 1.0+d, 1.0-x, &F2);

    result->val  = pre1.val*F1.val + pre2.val*F2.val;
    result->err  = fabs(pre1.val*F1.err) + fabs(pre2.val*F2.err);
    result->err += fabs(pre1.err*F1.val) + fabs(pre2.err*F2.val);
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(pre1.val*F1.val) + fabs(pre2.val*F2.val));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;
  }
}

static int pow_omx(const double x, const double p, gsl_sf_result * result)
{
  double ln_omx;
  double ln_result;
  if(fabs(x) < GSL_ROOT5_DBL_EPSILON) {
    ln_omx = -x*(1.0 + x*(1.0/2.0 + x*(1.0/3.0 + x/4.0 + x*x/5.0)));
  }
  else {
    ln_omx = log(1.0-x);
  }
  ln_result = p * ln_omx;
  return gsl_sf_exp_err_e(ln_result, GSL_DBL_EPSILON * fabs(ln_result), result);
}

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

int
gsl_sf_hyperg_2F1_e(double a, double b, const double c,
                       const double x,
                       gsl_sf_result * result)
{
  const double d = c - a - b;
  const double rinta = floor(a + 0.5);
  const double rintb = floor(b + 0.5);
  const double rintc = floor(c + 0.5);
  const int a_neg_integer = ( a < 0.0  &&  fabs(a - rinta) < locEPS );
  const int b_neg_integer = ( b < 0.0  &&  fabs(b - rintb) < locEPS );
  const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );

  result->val = 0.0;
  result->err = 0.0;

   /* Handle x == 1.0 RJM */

  if (fabs (x - 1.0) < locEPS && (c - a - b) > 0 && c != 0 && !c_neg_integer) {
    gsl_sf_result lngamc, lngamcab, lngamca, lngamcb;
    double lngamc_sgn, lngamca_sgn, lngamcb_sgn;
    int status;
    int stat1 = gsl_sf_lngamma_sgn_e (c, &lngamc, &lngamc_sgn);
    int stat2 = gsl_sf_lngamma_e (c - a - b, &lngamcab);
    int stat3 = gsl_sf_lngamma_sgn_e (c - a, &lngamca, &lngamca_sgn);
    int stat4 = gsl_sf_lngamma_sgn_e (c - b, &lngamcb, &lngamcb_sgn);

    if (stat1 != GSL_SUCCESS || stat2 != GSL_SUCCESS
        || stat3 != GSL_SUCCESS || stat4 != GSL_SUCCESS)
      {
        DOMAIN_ERROR (result);
      }

    status =
      gsl_sf_exp_err_e (lngamc.val + lngamcab.val - lngamca.val - lngamcb.val,
                        lngamc.err + lngamcab.err + lngamca.err + lngamcb.err,
                        result);

    result->val *= lngamc_sgn / (lngamca_sgn * lngamcb_sgn);
      return status;
  }

  if(x < -1.0 || 1.0 <= x) {
    DOMAIN_ERROR(result);
  }

  if(c_neg_integer) {
    /* If c is a negative integer, then either a or b must be a
       negative integer of smaller magnitude than c to ensure
       cancellation of the series. */
    if(! (a_neg_integer && a > c + 0.1) && ! (b_neg_integer && b > c + 0.1)) {
      DOMAIN_ERROR(result);
    }
  }

  if(fabs(c-b) < locEPS || fabs(c-a) < locEPS) {
    return pow_omx(x, d, result);  /* (1-x)^(c-a-b) */
  }

  if(a >= 0.0 && b >= 0.0 && c >=0.0 && x >= 0.0 && x < 0.995) {
    /* Series has all positive definite
     * terms and x is not close to 1.
     */
    return hyperg_2F1_series(a, b, c, x, result);
  }

  if(fabs(a) < 10.0 && fabs(b) < 10.0) {
    /* a and b are not too large, so we attempt
     * variations on the series summation.
     */
    if(a_neg_integer) {
      return hyperg_2F1_series(rinta, b, c, x, result);
    }
    if(b_neg_integer) {
      return hyperg_2F1_series(a, rintb, c, x, result);
    }

    if(x < -0.25) {
      return hyperg_2F1_luke(a, b, c, x, result);
    }
    else if(x < 0.5) {
      return hyperg_2F1_series(a, b, c, x, result);
    }
    else {
      if(fabs(c) > 10.0) {
        return hyperg_2F1_series(a, b, c, x, result);
      }
      else {
        return hyperg_2F1_reflect(a, b, c, x, result);
      }
    }
  }
  else {
    /* Either a or b or both large.
     * Introduce some new variables ap,bp so that bp is
     * the larger in magnitude.
     */
    double ap, bp;
    if(fabs(a) > fabs(b)) {
      bp = a;
      ap = b;
    }
    else {
      bp = b;
      ap = a;
    }

    if(x < 0.0) {
      /* What the hell, maybe Luke will converge.
       */
      return hyperg_2F1_luke(a, b, c, x, result);
    }

    if(GSL_MAX_DBL(fabs(ap),1.0)*fabs(bp)*fabs(x) < 2.0*fabs(c)) {
      /* If c is large enough or x is small enough,
       * we can attempt the series anyway.
       */
      return hyperg_2F1_series(a, b, c, x, result);
    }

    if(fabs(bp*bp*x*x) < 0.001*fabs(bp) && fabs(ap) < 10.0) {
      /* The famous but nearly worthless "large b" asymptotic.
       */
      int stat = gsl_sf_hyperg_1F1_e(ap, c, bp*x, result);
      result->err = 0.001 * fabs(result->val);
      return stat;
    }

    /* We give up. */
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR ("error", GSL_EUNIMPL);
  }
}

int
gsl_sf_hyperg_2F1_conj_e(const double aR, const double aI, const double c,
                            const double x,
                            gsl_sf_result * result)
{
  const double ax = fabs(x);
  const double rintc = floor(c + 0.5);
  const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );

  result->val = 0.0;
  result->err = 0.0;

  if(ax >= 1.0 || c_neg_integer || c == 0.0) {
    DOMAIN_ERROR(result);
  }

  if(   (ax < 0.25 && fabs(aR) < 20.0 && fabs(aI) < 20.0)
     || (c > 0.0 && x > 0.0)
    ) {
    return hyperg_2F1_conj_series(aR, aI, c, x, result);
  }
  else if(fabs(aR) < 10.0 && fabs(aI) < 10.0) {
    if(x < -0.25) {
      return hyperg_2F1_conj_luke(aR, aI, c, x, result);
    }
    else {
      return hyperg_2F1_conj_series(aR, aI, c, x, result);
    }
  }
  else {
    if(x < 0.0) {
      /* What the hell, maybe Luke will converge.
       */
      return hyperg_2F1_conj_luke(aR, aI, c, x, result);
    }

    /* Give up. */
    result->val = 0.0;
    result->err = 0.0;
    GSL_ERROR ("error", GSL_EUNIMPL);
  }
}

int
gsl_sf_hyperg_2F1_renorm_e(const double a, const double b, const double c,
                              const double x,
                              gsl_sf_result * result
                              )
{
  const double rinta = floor(a + 0.5);
  const double rintb = floor(b + 0.5);
  const double rintc = floor(c + 0.5);
  const int a_neg_integer = ( a < 0.0  &&  fabs(a - rinta) < locEPS );
  const int b_neg_integer = ( b < 0.0  &&  fabs(b - rintb) < locEPS );
  const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );

  if(c_neg_integer) {
    if((a_neg_integer && a > c+0.1) || (b_neg_integer && b > c+0.1)) {
      /* 2F1 terminates early */
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else {
      /* 2F1 does not terminate early enough, so something survives */
      /* [Abramowitz+Stegun, 15.1.2] */
      gsl_sf_result g1, g2, g3, g4, g5;
      double s1, s2, s3, s4, s5;
      int stat = 0;
      stat += gsl_sf_lngamma_sgn_e(a-c+1, &g1, &s1);
      stat += gsl_sf_lngamma_sgn_e(b-c+1, &g2, &s2);
      stat += gsl_sf_lngamma_sgn_e(a, &g3, &s3);
      stat += gsl_sf_lngamma_sgn_e(b, &g4, &s4);
      stat += gsl_sf_lngamma_sgn_e(-c+2, &g5, &s5);
      if(stat != 0) {
        DOMAIN_ERROR(result);
      }
      else {
        gsl_sf_result F;
        int stat_F = gsl_sf_hyperg_2F1_e(a-c+1, b-c+1, -c+2, x, &F);
        double ln_pre_val = g1.val + g2.val - g3.val - g4.val - g5.val;
        double ln_pre_err = g1.err + g2.err + g3.err + g4.err + g5.err;
        double sg  = s1 * s2 * s3 * s4 * s5;
        int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
                                              sg * F.val, F.err,
                                              result);
        return GSL_ERROR_SELECT_2(stat_e, stat_F);
      }
    }
  }
  else {
    /* generic c */
    gsl_sf_result F;
    gsl_sf_result lng;
    double sgn;
    int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);
    int stat_F = gsl_sf_hyperg_2F1_e(a, b, c, x, &F);
    int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,
                                          sgn*F.val, F.err,
                                          result);
    return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);
  }
}

int
gsl_sf_hyperg_2F1_conj_renorm_e(const double aR, const double aI, const double c,
                                   const double x,
                                   gsl_sf_result * result
                                   )
{
  const double rintc = floor(c  + 0.5);
  const double rinta = floor(aR + 0.5);
  const int a_neg_integer = ( aR < 0.0 && fabs(aR-rinta) < locEPS && aI == 0.0);
  const int c_neg_integer = (  c < 0.0 && fabs(c - rintc) < locEPS );

  if(c_neg_integer) {
    if(a_neg_integer && aR > c+0.1) {
      /* 2F1 terminates early */
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else {
      /* 2F1 does not terminate early enough, so something survives */
      /* [Abramowitz+Stegun, 15.1.2] */
      gsl_sf_result g1, g2;
      gsl_sf_result g3;
      gsl_sf_result a1, a2;
      int stat = 0;
      stat += gsl_sf_lngamma_complex_e(aR-c+1, aI, &g1, &a1);
      stat += gsl_sf_lngamma_complex_e(aR, aI, &g2, &a2);
      stat += gsl_sf_lngamma_e(-c+2.0, &g3);
      if(stat != 0) {
        DOMAIN_ERROR(result);
      }
      else {
        gsl_sf_result F;
        int stat_F = gsl_sf_hyperg_2F1_conj_e(aR-c+1, aI, -c+2, x, &F);
        double ln_pre_val = 2.0*(g1.val - g2.val) - g3.val;
        double ln_pre_err = 2.0 * (g1.err + g2.err) + g3.err;
        int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
                                              F.val, F.err,
                                              result);
        return GSL_ERROR_SELECT_2(stat_e, stat_F);
      }
    }
  }
  else {
    /* generic c */
    gsl_sf_result F;
    gsl_sf_result lng;
    double sgn;
    int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);
    int stat_F = gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &F);
    int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,
                                          sgn*F.val, F.err,
                                          result);
    return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);
  }
}

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

#include "eval.h"

double gsl_sf_hyperg_2F1(double a, double b, double c, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_2F1_e(a, b, c, x, &result));
}

double gsl_sf_hyperg_2F1_conj(double aR, double aI, double c, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &result));
}

double gsl_sf_hyperg_2F1_renorm(double a, double b, double c, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_2F1_renorm_e(a, b, c, x, &result));
}

double gsl_sf_hyperg_2F1_conj_renorm(double aR, double aI, double c, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_2F1_conj_renorm_e(aR, aI, c, x, &result));
}

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
/* monte/vegas.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Michael Booth
 *
 * 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.
 */

/* Author: MJB */
/* Modified by: Brian Gough, 12/2000 */

/* This is an implementation of the adaptive Monte-Carlo algorithm
   of G. P. Lepage, originally described in J. Comp. Phys. 27, 192(1978).
   The current version of the algorithm was described in the Cornell
   preprint CLNS-80/447 of March, 1980.

   This code follows most closely the c version by D.R.Yennie, coded
   in 1984.

   The input coordinates are x[j], with upper and lower limits xu[j]
   and xl[j].  The integration length in the j-th direction is
   delx[j].  Each coordinate x[j] is rescaled to a variable y[j] in
   the range 0 to 1.  The range is divided into bins with boundaries
   xi[i][j], where i=0 corresponds to y=0 and i=bins to y=1.  The grid
   is refined (ie, bins are adjusted) using d[i][j] which is some
   variation on the squared sum.  A third parameter used in defining
   the real coordinate using random numbers is called z.  It ranges
   from 0 to bins.  Its integer part gives the lower index of the bin
   into which a call is to be placed, and the remainder gives the
   location inside the bin.

   When stratified sampling is used the bins are grouped into boxes,
   and the algorithm allocates an equal number of function calls to
   each box.

   The variable alpha controls how "stiff" the rebinning algorithm is.
   alpha = 0 means never change the grid.  Alpha is typically set between
   1 and 2.

   */

/* configuration headers */
#include <config.h>

/* standard headers */
#include <math.h>
#include <stdio.h>

/* gsl headers */
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_monte_vegas.h>

/* lib-specific headers */
#define BINS_MAX 50		/* even integer, will be divided by two */

/* A separable grid with coordinates and values */
#define COORD(s,i,j) ((s)->xi[(i)*(s)->dim + (j)])
#define NEW_COORD(s,i) ((s)->xin[(i)])
#define VALUE(s,i,j) ((s)->d[(i)*(s)->dim + (j)])

/* predeclare functions */

typedef int coord;

static void init_grid (gsl_monte_vegas_state * s, double xl[], double xu[],
		size_t dim);
static void reset_grid_values (gsl_monte_vegas_state * s);
static void init_box_coord (gsl_monte_vegas_state * s, coord box[]);
static int change_box_coord (gsl_monte_vegas_state * s, coord box[]);
static void accumulate_distribution (gsl_monte_vegas_state * s, coord bin[],
                                     double y);
static void random_point (double x[], coord bin[], double *bin_vol,
                          const coord box[],
                          const double xl[], const double xu[],
                          gsl_monte_vegas_state * s, gsl_rng * r);
static void resize_grid (gsl_monte_vegas_state * s, unsigned int bins);
static void refine_grid (gsl_monte_vegas_state * s);

static void print_lim (gsl_monte_vegas_state * state,
                       double xl[], double xu[], unsigned long dim);
static void print_head (gsl_monte_vegas_state * state,
                        unsigned long num_dim, unsigned long calls,
                        unsigned int it_num,
                        unsigned int bins, unsigned int boxes);
static void print_res (gsl_monte_vegas_state * state,
                       unsigned int itr, double res, double err,
                       double cum_res, double cum_err,
                       double chi_sq);
static void print_dist (gsl_monte_vegas_state * state, unsigned long dim);
static void print_grid (gsl_monte_vegas_state * state, unsigned long dim);

int
gsl_monte_vegas_integrate (gsl_monte_function * f,
			   double xl[], double xu[],
			   size_t dim, size_t calls,
			   gsl_rng * r,
			   gsl_monte_vegas_state * state,
			   double *result, double *abserr)
{
  double cum_int, cum_sig;
  size_t i, k, it;

  if (dim != state->dim)
    {
      GSL_ERROR ("number of dimensions must match allocated size", GSL_EINVAL);
    }

  for (i = 0; i < dim; i++)
    {
      if (xu[i] <= xl[i])
	{
	  GSL_ERROR ("xu must be greater than xl", GSL_EINVAL);
	}

      if (xu[i] - xl[i] > GSL_DBL_MAX)
	{
	  GSL_ERROR ("Range of integration is too large, please rescale",
		     GSL_EINVAL);
	}
    }

  if (state->stage == 0)
    {
      init_grid (state, xl, xu, dim);

      if (state->verbose >= 0)
	{
	  print_lim (state, xl, xu, dim);
	}
    }

  if (state->stage <= 1)
    {
      state->wtd_int_sum = 0;
      state->sum_wgts = 0;
      state->chi_sum = 0;
      state->it_num = 1;
      state->samples = 0;
    }

  if (state->stage <= 2)
    {
      unsigned int bins = state->bins_max;
      unsigned int boxes = 1;

      if (state->mode != GSL_VEGAS_MODE_IMPORTANCE_ONLY)
	{
	  /* shooting for 2 calls/box */

	  boxes = floor (pow (calls / 2.0, 1.0 / dim));
	  state->mode = GSL_VEGAS_MODE_IMPORTANCE;

	  if (2 * boxes >= state->bins_max)
	    {
	      /* if bins/box < 2 */
	      int box_per_bin = GSL_MAX (boxes / state->bins_max, 1);

	      bins = GSL_MIN(boxes / box_per_bin, state->bins_max);
	      boxes = box_per_bin * bins;

	      state->mode = GSL_VEGAS_MODE_STRATIFIED;
	    }
	}

      {
	double tot_boxes = pow ((double) boxes, (double) dim);
	state->calls_per_box = GSL_MAX (calls / tot_boxes, 2);
	calls = state->calls_per_box * tot_boxes;
      }

      /* total volume of x-space/(avg num of calls/bin) */
      state->jac = state->vol * pow ((double) bins, (double) dim) / calls;

      state->boxes = boxes;

      /* If the number of bins changes from the previous invocation, bins
         are expanded or contracted accordingly, while preserving bin
         density */

      if (bins != state->bins)
	{
	  resize_grid (state, bins);

	  if (state->verbose > 1)
	    {
	      print_grid (state, dim);
	    }
	}

      if (state->verbose >= 0)
	{
	  print_head (state,
		      dim, calls, state->it_num, state->bins, state->boxes);
	}
    }

  state->it_start = state->it_num;

  cum_int = 0.0;
  cum_sig = 0.0;

  state->chisq = 0.0;

  for (it = 0; it < state->iterations; it++)
    {
      double intgrl = 0.0, intgrl_sq = 0.0;
      double sig = 0.0;
      double wgt;
      size_t calls_per_box = state->calls_per_box;
      double jacbin = state->jac;
      double *x = state->x;
      coord *bin = state->bin;

      state->it_num = state->it_start + it;

      reset_grid_values (state);
      init_box_coord (state, state->box);

      do
	{
          double m = 0, q = 0;
	  double f_sq_sum = 0.0;

	  for (k = 0; k < calls_per_box; k++)
	    {
	      double fval, bin_vol;

	      random_point (x, bin, &bin_vol, state->box, xl, xu, state, r);

	      fval = jacbin * bin_vol * GSL_MONTE_FN_EVAL (f, x);

              /* recurrence for mean and variance */

              {
                double d = fval - m;
                m += d / (k + 1.0);
                q += d * d * (k / (k + 1.0));
              }

	      if (state->mode != GSL_VEGAS_MODE_STRATIFIED)
		{
                  double f_sq = fval * fval;
		  accumulate_distribution (state, bin, f_sq);
		}
	    }

	  intgrl += m * calls_per_box;

	  f_sq_sum = q * calls_per_box ;

	  sig += f_sq_sum ;

	  if (state->mode == GSL_VEGAS_MODE_STRATIFIED)
	    {
	      accumulate_distribution (state, bin, f_sq_sum);
	    }
	}
      while (change_box_coord (state, state->box));

      /* Compute final results for this iteration   */

      sig = sig / (calls_per_box - 1.0)  ;

      if (sig > 0)
        {
          wgt = 1.0 / sig;
        }
      else if (state->sum_wgts > 0)
        {
          wgt = state->sum_wgts / state->samples;
        }
      else
        {
          wgt = 0.0;
        }

     intgrl_sq = intgrl * intgrl;

     state->result = intgrl;
     state->sigma  = sqrt(sig);

     if (wgt > 0.0)
       {
         state->samples++ ;
         state->sum_wgts += wgt;
         state->wtd_int_sum += intgrl * wgt;
         state->chi_sum += intgrl_sq * wgt;

         cum_int = state->wtd_int_sum / state->sum_wgts;
         cum_sig = sqrt (1 / state->sum_wgts);

         if (state->samples > 1)
           {
             state->chisq = (state->chi_sum - state->wtd_int_sum * cum_int) /
               (state->samples - 1.0);
           }
       }
     else
       {
         cum_int += (intgrl - cum_int) / (it + 1.0);
         cum_sig = 0.0;
       }

      if (state->verbose >= 0)
	{
	  print_res (state,
		     state->it_num, intgrl, sqrt (sig), cum_int, cum_sig,
		     state->chisq);
	  if (it + 1 == state->iterations && state->verbose > 0)
	    {
	      print_grid (state, dim);
	    }
	}

      if (state->verbose > 1)
	{
          print_dist (state, dim);
        }

      refine_grid (state);

      if (state->verbose > 1)
	{
	  print_grid (state, dim);
	}

    }

  /* By setting stage to 1 further calls will generate independent
     estimates based on the same grid, although it may be rebinned. */

  state->stage = 1;

  *result = cum_int;
  *abserr = cum_sig;

  return GSL_SUCCESS;
}

gsl_monte_vegas_state *
gsl_monte_vegas_alloc (size_t dim)
{
  gsl_monte_vegas_state *s =
    (gsl_monte_vegas_state *) malloc (sizeof (gsl_monte_vegas_state));

  if (s == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for vegas state struct",
		     GSL_ENOMEM, 0);
    }

  s->delx = (double *) malloc (dim * sizeof (double));

  if (s->delx == 0)
    {
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for delx", GSL_ENOMEM, 0);
    }

  s->d = (double *) malloc (BINS_MAX * dim * sizeof (double));

  if (s->d == 0)
    {
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for d", GSL_ENOMEM, 0);
    }

  s->xi = (double *) malloc ((BINS_MAX + 1) * dim * sizeof (double));

  if (s->xi == 0)
    {
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for xi", GSL_ENOMEM, 0);
    }

  s->xin = (double *) malloc ((BINS_MAX + 1) * sizeof (double));

  if (s->xin == 0)
    {
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0);
    }

  s->weight = (double *) malloc (BINS_MAX * sizeof (double));

  if (s->weight == 0)
    {
      free (s->xin);
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for xin", GSL_ENOMEM, 0);
    }

  s->box = (coord *) malloc (dim * sizeof (coord));

  if (s->box == 0)
    {
      free (s->weight);
      free (s->xin);
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for box", GSL_ENOMEM, 0);
    }

  s->bin = (coord *) malloc (dim * sizeof (coord));

  if (s->bin == 0)
    {
      free (s->box);
      free (s->weight);
      free (s->xin);
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for bin", GSL_ENOMEM, 0);
    }

  s->x = (double *) malloc (dim * sizeof (double));

  if (s->x == 0)
    {
      free (s->bin);
      free (s->box);
      free (s->weight);
      free (s->xin);
      free (s->xi);
      free (s->d);
      free (s->delx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
    }

  s->dim = dim;
  s->bins_max = BINS_MAX;

  gsl_monte_vegas_init (s);

  return s;
}

/* Set some default values and whatever */
int
gsl_monte_vegas_init (gsl_monte_vegas_state * state)
{
  state->stage = 0;
  state->alpha = 1.5;
  state->verbose = -1;
  state->iterations = 5;
  state->mode = GSL_VEGAS_MODE_IMPORTANCE;
  state->chisq = 0;
  state->bins = state->bins_max;
  state->ostream = stdout;

  return GSL_SUCCESS;
}

void
gsl_monte_vegas_free (gsl_monte_vegas_state * s)
{
  free (s->x);
  free (s->delx);
  free (s->d);
  free (s->xi);
  free (s->xin);
  free (s->weight);
  free (s->box);
  free (s->bin);
  free (s);
}

static void
init_box_coord (gsl_monte_vegas_state * s, coord box[])
{
  size_t i;

  size_t dim = s->dim;

  for (i = 0; i < dim; i++)
    {
      box[i] = 0;
    }
}

/* change_box_coord steps through the box coord like
   {0,0}, {0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, ...
*/
static int
change_box_coord (gsl_monte_vegas_state * s, coord box[])
{
  int j = s->dim - 1;

  int ng = s->boxes;

  while (j >= 0)
    {
      box[j] = (box[j] + 1) % ng;

      if (box[j] != 0)
	{
	  return 1;
	}

      j--;
    }

  return 0;
}

static void
init_grid (gsl_monte_vegas_state * s, double xl[], double xu[], size_t dim)
{
  size_t j;

  double vol = 1.0;

  s->bins = 1;

  for (j = 0; j < dim; j++)
    {
      double dx = xu[j] - xl[j];
      s->delx[j] = dx;
      vol *= dx;

      COORD (s, 0, j) = 0.0;
      COORD (s, 1, j) = 1.0;
    }

  s->vol = vol;
}

static void
reset_grid_values (gsl_monte_vegas_state * s)
{
  size_t i, j;

  size_t dim = s->dim;
  size_t bins = s->bins;

  for (i = 0; i < bins; i++)
    {
      for (j = 0; j < dim; j++)
	{
	  VALUE (s, i, j) = 0.0;
	}
    }
}

static void
accumulate_distribution (gsl_monte_vegas_state * s, coord bin[], double y)
{
  size_t j;
  size_t dim = s->dim;

  for (j = 0; j < dim; j++)
    {
      int i = bin[j];
      VALUE (s, i, j) += y;
    }
}

static void
random_point (double x[], coord bin[], double *bin_vol,
	      const coord box[], const double xl[], const double xu[],
	      gsl_monte_vegas_state * s, gsl_rng * r)
{
  /* Use the random number generator r to return a random position x
     in a given box.  The value of bin gives the bin location of the
     random position (there may be several bins within a given box) */

  double vol = 1.0;

  size_t j;

  size_t dim = s->dim;
  size_t bins = s->bins;
  size_t boxes = s->boxes;

  DISCARD_POINTER(xu); /* prevent warning about unused parameter */

  for (j = 0; j < dim; ++j)
    {
      /* box[j] + ran gives the position in the box units, while z
         is the position in bin units.  */

      double z = ((box[j] + gsl_rng_uniform_pos (r)) / boxes) * bins;

      int k = z;

      double y, bin_width;

      bin[j] = k;

      if (k == 0)
	{
	  bin_width = COORD (s, 1, j);
	  y = z * bin_width;
	}
      else
	{
	  bin_width = COORD (s, k + 1, j) - COORD (s, k, j);
	  y = COORD (s, k, j) + (z - k) * bin_width;
	}

      x[j] = xl[j] + y * s->delx[j];

      vol *= bin_width;
    }

  *bin_vol = vol;
}

static void
resize_grid (gsl_monte_vegas_state * s, unsigned int bins)
{
  size_t j, k;
  size_t dim = s->dim;

  /* weight is ratio of bin sizes */

  double pts_per_bin = (double) s->bins / (double) bins;

  for (j = 0; j < dim; j++)
    {
      double xold;
      double xnew = 0;
      double dw = 0;
      int i = 1;

      for (k = 1; k <= s->bins; k++)
	{
	  dw += 1.0;
	  xold = xnew;
	  xnew = COORD (s, k, j);

	  for (; dw > pts_per_bin; i++)
	    {
	      dw -= pts_per_bin;
	      NEW_COORD (s, i) = xnew - (xnew - xold) * dw;
	    }
	}

      for (k = 1 ; k < bins; k++)
        {
          COORD(s, k, j) = NEW_COORD(s, k);
        }

      COORD (s, bins, j) = 1;
    }

  s->bins = bins;
}

static void
refine_grid (gsl_monte_vegas_state * s)
{
  size_t i, j, k;
  size_t dim = s->dim;
  size_t bins = s->bins;

  for (j = 0; j < dim; j++)
    {
      double grid_tot_j, tot_weight;
      double * weight = s->weight;

      double oldg = VALUE (s, 0, j);
      double newg = VALUE (s, 1, j);

      VALUE (s, 0, j) = (oldg + newg) / 2;
      grid_tot_j = VALUE (s, 0, j);

      /* This implements gs[i][j] = (gs[i-1][j]+gs[i][j]+gs[i+1][j])/3 */

      for (i = 1; i < bins - 1; i++)
	{
	  double rc = oldg + newg;
	  oldg = newg;
	  newg = VALUE (s, i + 1, j);
	  VALUE (s, i, j) = (rc + newg) / 3;
	  grid_tot_j += VALUE (s, i, j);
	}
      VALUE (s, bins - 1, j) = (newg + oldg) / 2;

      grid_tot_j += VALUE (s, bins - 1, j);

      tot_weight = 0;

      for (i = 0; i < bins; i++)
	{
	  weight[i] = 0;

	  if (VALUE (s, i, j) > 0)
	    {
	      oldg = grid_tot_j / VALUE (s, i, j);
	      /* damped change */
	      weight[i] = pow (((oldg - 1) / oldg / log (oldg)), s->alpha);
	    }

	  tot_weight += weight[i];

#ifdef DEBUG
          printf("weight[%d] = %g\n", i, weight[i]);
#endif
	}

      {
	double pts_per_bin = tot_weight / bins;

	double xold;
	double xnew = 0;
	double dw = 0;
	i = 1;

	for (k = 0; k < bins; k++)
	  {
	    dw += weight[k];
	    xold = xnew;
	    xnew = COORD (s, k + 1, j);

	    for (; dw > pts_per_bin; i++)
	      {
		dw -= pts_per_bin;
		NEW_COORD (s, i) = xnew - (xnew - xold) * dw / weight[k];
	      }
	  }

        for (k = 1 ; k < bins ; k++)
          {
            COORD(s, k, j) = NEW_COORD(s, k);
          }

	COORD (s, bins, j) = 1;
      }
    }
}

static void
print_lim (gsl_monte_vegas_state * state,
	   double xl[], double xu[], unsigned long dim)
{
  unsigned long j;

  fprintf (state->ostream, "The limits of integration are:\n");
  for (j = 0; j < dim; ++j)
    fprintf (state->ostream, "\nxl[%lu]=%f    xu[%lu]=%f", j, xl[j], j, xu[j]);
  fprintf (state->ostream, "\n");
  fflush (state->ostream);
}

static void
print_head (gsl_monte_vegas_state * state,
	    unsigned long num_dim, unsigned long calls,
	    unsigned int it_num, unsigned int bins, unsigned int boxes)
{
  fprintf (state->ostream,
	   "\nnum_dim=%lu, calls=%lu, it_num=%d, max_it_num=%d ",
	   num_dim, calls, it_num, state->iterations);
  fprintf (state->ostream,
	   "verb=%d, alph=%.2f,\nmode=%d, bins=%d, boxes=%d\n",
	   state->verbose, state->alpha, state->mode, bins, boxes);
  fprintf (state->ostream,
	   "\n       single.......iteration                   ");
  fprintf (state->ostream, "accumulated......results   \n");

  fprintf (state->ostream,
	   "iteration     integral    sigma             integral   ");
  fprintf (state->ostream, "      sigma     chi-sq/it\n\n");
  fflush (state->ostream);

}

static void
print_res (gsl_monte_vegas_state * state,
	   unsigned int itr,
           double res, double err,
           double cum_res, double cum_err,
	   double chi_sq)
{
  fprintf (state->ostream,
	   "%4d        %6.4e %10.2e          %6.4e      %8.2e  %10.2e\n",
	   itr, res, err, cum_res, cum_err, chi_sq);
  fflush (state->ostream);
}

static void
print_dist (gsl_monte_vegas_state * state, unsigned long dim)
{
  unsigned long i, j;
  int p = state->verbose;
  if (p < 1)
    return;

  for (j = 0; j < dim; ++j)
    {
      fprintf (state->ostream, "\n axis %lu \n", j);
      fprintf (state->ostream, "      x   g\n");
      for (i = 0; i < state->bins; i++)
	{
	  fprintf (state->ostream, "weight [%11.2e , %11.2e] = ",
                   COORD (state, i, j), COORD(state,i+1,j));
	  fprintf (state->ostream, " %11.2e\n", VALUE (state, i, j));

	}
      fprintf (state->ostream, "\n");
    }
  fprintf (state->ostream, "\n");
  fflush (state->ostream);

}

static void
print_grid (gsl_monte_vegas_state * state, unsigned long dim)
{
  unsigned long i, j;
  int p = state->verbose;
  if (p < 1)
    return;

  for (j = 0; j < dim; ++j)
    {
      fprintf (state->ostream, "\n axis %lu \n", j);
      fprintf (state->ostream, "      x   \n");
      for (i = 0; i <= state->bins; i++)
	{
	  fprintf (state->ostream, "%11.2e", COORD (state, i, j));
	  if (i % 5 == 4)
	    fprintf (state->ostream, "\n");
	}
      fprintf (state->ostream, "\n");
    }
  fprintf (state->ostream, "\n");
  fflush (state->ostream);

}


text:
/* These tests are based on the NIST Statistical Reference Datasets
   See http://www.nist.gov/itl/div898/strd/index.html for more
   information. */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_fit.h>

#include <gsl/gsl_ieee_utils.h>

size_t norris_n = 36;

double norris_x[] = { 0.2, 337.4, 118.2, 884.6, 10.1, 226.5, 666.3, 996.3,
                      448.6, 777.0, 558.2, 0.4, 0.6, 775.5, 666.9, 338.0,
                      447.5, 11.6, 556.0, 228.1, 995.8, 887.6, 120.2, 0.3,
                      0.3, 556.8, 339.1, 887.2, 999.0, 779.0, 11.1, 118.3,
                      229.2, 669.1, 448.9, 0.5 } ;

double norris_y[] = { 0.1, 338.8, 118.1, 888.0, 9.2, 228.1, 668.5, 998.5,
                      449.1, 778.9, 559.2, 0.3, 0.1, 778.1, 668.8, 339.3,
                      448.9, 10.8, 557.7, 228.3, 998.0, 888.8, 119.6, 0.3,
                      0.6, 557.6, 339.3, 888.0, 998.5, 778.9, 10.2, 117.6,
                      228.9, 668.4, 449.2, 0.2};

size_t noint1_n = 11;
double noint1_x[] = { 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70 };
double noint1_y[] = { 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140};

size_t noint2_n = 3;
double noint2_x[] = { 4, 5, 6 } ;
double noint2_y[] = { 3, 4, 4 } ;

int
main (void)
{

  double x[1000], y[1000], w[1000];

  size_t xstride = 2, wstride = 3, ystride = 5;
  size_t i;

  for (i = 0; i < norris_n; i++)
    {
      x[i*xstride] = norris_x[i];
      w[i*wstride] = 1.0;
      y[i*ystride] = norris_y[i];
    }

  gsl_ieee_env_setup();

  {
    double c0, c1, cov00, cov01, cov11, sumsq;

    double expected_c0 = -0.262323073774029;
    double expected_c1 =  1.00211681802045;
    double expected_cov00 = pow(0.232818234301152, 2.0);
    double expected_cov01 = -7.74327536339570e-05;  /* computed from octave */
    double expected_cov11 = pow(0.429796848199937E-03, 2.0);
    double expected_sumsq = 26.6173985294224;

    gsl_fit_linear (x, xstride, y, ystride, norris_n,
                    &c0, &c1, &cov00, &cov01, &cov11, &sumsq);

    /* gsl_fit_wlinear (x, xstride, w, wstride, y, ystride, norris_n,
                     &c0, &c1, &cov00, &cov01, &cov11, &sumsq); */

    gsl_test_rel (c0, expected_c0, 1e-10, "norris gsl_fit_linear c0") ;
    gsl_test_rel (c1, expected_c1, 1e-10, "norris gsl_fit_linear c1") ;
    gsl_test_rel (cov00, expected_cov00, 1e-10, "norris gsl_fit_linear cov00") ;
    gsl_test_rel (cov01, expected_cov01, 1e-10, "norris gsl_fit_linear cov01") ;
    gsl_test_rel (cov11, expected_cov11, 1e-10, "norris gsl_fit_linear cov11") ;
    gsl_test_rel (sumsq, expected_sumsq, 1e-10, "norris gsl_fit_linear sumsq") ;
  }

  {
    double c0, c1, cov00, cov01, cov11, sumsq;

    double expected_c0 = -0.262323073774029;
    double expected_c1 =  1.00211681802045;
    double expected_cov00 = 6.92384428759429e-02;  /* computed from octave */
    double expected_cov01 = -9.89095016390515e-05; /* computed from octave */
    double expected_cov11 = 2.35960747164148e-07;  /* computed from octave */
    double expected_sumsq = 26.6173985294224;

    gsl_fit_wlinear (x, xstride, w, wstride, y, ystride, norris_n,
                     &c0, &c1, &cov00, &cov01, &cov11, &sumsq);

    gsl_test_rel (c0, expected_c0, 1e-10, "norris gsl_fit_wlinear c0") ;
    gsl_test_rel (c1, expected_c1, 1e-10, "norris gsl_fit_wlinear c1") ;
    gsl_test_rel (cov00, expected_cov00, 1e-10, "norris gsl_fit_wlinear cov00") ;
    gsl_test_rel (cov01, expected_cov01, 1e-10, "norris gsl_fit_wlinear cov01") ;
    gsl_test_rel (cov11, expected_cov11, 1e-10, "norris gsl_fit_wlinear cov11") ;
    gsl_test_rel (sumsq, expected_sumsq, 1e-10, "norris gsl_fit_wlinear sumsq") ;
  }

  for (i = 0; i < noint1_n; i++)
    {
      x[i*xstride] = noint1_x[i];
      w[i*wstride] = 1.0;
      y[i*ystride] = noint1_y[i];
    }

  {
    double c1, cov11, sumsq;

    double expected_c1 = 2.07438016528926;
    double expected_cov11 = pow(0.165289256198347E-01, 2.0);
    double expected_sumsq = 127.272727272727;

    gsl_fit_mul (x, xstride, y, ystride, noint1_n, &c1, &cov11, &sumsq);

    gsl_test_rel (c1, expected_c1, 1e-10, "noint1 gsl_fit_mul c1") ;
    gsl_test_rel (cov11, expected_cov11, 1e-10, "noint1 gsl_fit_mul cov11") ;
    gsl_test_rel (sumsq, expected_sumsq, 1e-10, "noint1 gsl_fit_mul sumsq") ;
  }

  {
    double c1, cov11, sumsq;

    double expected_c1 = 2.07438016528926;
    double expected_cov11 = 2.14661371686165e-05; /* computed from octave */
    double expected_sumsq = 127.272727272727;

    gsl_fit_wmul (x, xstride, w, wstride, y, ystride, noint1_n, &c1, &cov11, &sumsq);

    gsl_test_rel (c1, expected_c1, 1e-10, "noint1 gsl_fit_wmul c1") ;
    gsl_test_rel (cov11, expected_cov11, 1e-10, "noint1 gsl_fit_wmul cov11") ;
    gsl_test_rel (sumsq, expected_sumsq, 1e-10, "noint1 gsl_fit_wmul sumsq") ;
  }

  for (i = 0; i < noint2_n; i++)
    {
      x[i*xstride] = noint2_x[i];
      w[i*wstride] = 1.0;
      y[i*ystride] = noint2_y[i];
    }

  {
    double c1, cov11, sumsq;

    double expected_c1 = 0.727272727272727;
    double expected_cov11 = pow(0.420827318078432E-01, 2.0);
    double expected_sumsq = 0.272727272727273;

    gsl_fit_mul (x, xstride, y, ystride, noint2_n, &c1, &cov11, &sumsq);

    gsl_test_rel (c1, expected_c1, 1e-10, "noint2 gsl_fit_mul c1") ;
    gsl_test_rel (cov11, expected_cov11, 1e-10, "noint2 gsl_fit_mul cov11") ;
    gsl_test_rel (sumsq, expected_sumsq, 1e-10, "noint2 gsl_fit_mul sumsq") ;
  }

  {
    double c1, cov11, sumsq;

    double expected_c1 = 0.727272727272727;
    double expected_cov11 = 1.29870129870130e-02 ; /* computed from octave */
    double expected_sumsq = 0.272727272727273;

    gsl_fit_wmul (x, xstride, w, wstride, y, ystride, noint2_n, &c1, &cov11, &sumsq);

    gsl_test_rel (c1, expected_c1, 1e-10, "noint2 gsl_fit_wmul c1") ;
    gsl_test_rel (cov11, expected_cov11, 1e-10, "noint2 gsl_fit_wmul cov11") ;
    gsl_test_rel (sumsq, expected_sumsq, 1e-10, "noint2 gsl_fit_wmul sumsq") ;
  }

  /* now summarize the results */

  exit (gsl_test_summary ());
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#pragma once

#include <string_view>
#include <gsl/span>
#include <spirv_cross.hpp>

namespace babylon
{
    class ShaderCompiler
    {
    public:
        ShaderCompiler();
        ~ShaderCompiler();

        struct ShaderInfo
        {
            std::unique_ptr<const spirv_cross::Compiler> Compiler;
            gsl::span<uint8_t> Bytes;
        };

        void Compile(std::string_view vertexSource, std::string_view fragmentSource, std::function<void (ShaderInfo, ShaderInfo)> onCompiled);
    };
}

max_stars_repo_name: TrevorDev/BabylonNative

max_issues_repo_name: TrevorDev/BabylonNative

max_forks_repo_name: ryantrem/BabylonNative


text:
/*
 * Licensed to the Apache Software Foundation (ASF) under one
 * or more contributor license agreements.  See the NOTICE file
 * distributed with this work for additional information
 * regarding copyright ownership.  The ASF licenses this file
 * to you under the Apache License, Version 2.0 (the
 * "License"); you may not use this file except in compliance
 * with the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing,
 * software distributed under the License is distributed on an
 * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
 * KIND, either express or implied.  See the License for the
 * specific language governing permissions and limitations
 * under the License.
 */

#ifndef _libmatrixmult_h
#define _libmatrixmult_h

#define MAX(x, y) (((x) > (y)) ? (x) : (y))
#define MIN(x, y) (((x) < (y)) ? (x) : (y))

// *****************************************************************
// We support Intel MKL (recommended) or OpenBLAS.
// These flags are used for conditional compilation with mkl and openblas
// #define USE_INTEL_MKL
// #define USE_GNU_THREADING
// #define USE_OPEN_BLAS
// *****************************************************************

//#ifdef __cplusplus
//extern "C" {
//#endif
//#ifdef __cplusplus
//}
//#endif

// Since we call cblas_dgemm in openmp for loop,
// we call "extension" APIs for setting number of threads of the given API.
// For example: for OpenBLAS we use openblas_set_num_threads and
// for MKL we use mkl_set_num_threads. This avoids performance degradation due to overprovisioning.
#ifdef USE_OPEN_BLAS
#include <cblas.h>
// extern "C" void openblas_set_num_threads(int numThreads);
#elif defined USE_INTEL_MKL
#include <mkl.h>
#include <mkl_service.h>
#endif

void setNumThreadsForBLAS(int numThreads);

// Multiplies two matrices m1Ptr and m2Ptr in row-major format of shape
// (m1rlen, m1clen) and (m1clen, m2clen)
void matmult(double* m1Ptr, double* m2Ptr, double* retPtr, int m1rlen,
             int m1clen, int m2clen, int numThreads);

void tsmm(double* m1Ptr, double* retPtr, int m1rlen, int m1clen, bool isLeftTranspose,  int numThreads);

#endif

max_stars_repo_name: frreiss/fred-systemml

max_issues_repo_name: frreiss/fred-systemml

max_forks_repo_name: frreiss/fred-systemml


text:
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

/* Computation of the integral,

      I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))

   over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
   is Gamma(1/4)^4/(4 pi^3).  This example is taken from
   C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
   Volume 1", Section 1.1, p21, which cites the original
   paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
   1800 (1977) */

/* For simplicity we compute the integral over the region
   (0,0,0) -> (pi,pi,pi) and multiply by 8 */

double exact = 1.3932039296856768591842462603255;

double
g (double *k, size_t dim, void *params)
{
  double A = 1.0 / (M_PI * M_PI * M_PI);
  return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
}

void
display_results (char *title, double result, double error)
{
  printf ("%s ==================\n", title);
  printf ("result = % .6f\n", result);
  printf ("sigma  = % .6f\n", error);
  printf ("exact  = % .6f\n", exact);
  printf ("error  = % .6f = %.2g sigma\n", result - exact,
          fabs (result - exact) / error);
}

int
main (void)
{
  double res, err;

  double xl[3] = { 0, 0, 0 };
  double xu[3] = { M_PI, M_PI, M_PI };

  const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G = { &g, 3, 0 };

  size_t calls = 500000;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  {
    gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
    gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s,
                               &res, &err);
    gsl_monte_plain_free (s);

    display_results ("plain", res, err);
  }

  {
    gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
    gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
                               &res, &err);
    gsl_monte_miser_free (s);

    display_results ("miser", res, err);
  }

  {
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);

    gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
                               &res, &err);
    display_results ("vegas warm-up", res, err);

    printf ("converging...\n");

    do
      {
        gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
                                   &res, &err);
        printf ("result = % .6f sigma = % .6f "
                "chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s));
      }
    while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);

    display_results ("vegas final", res, err);

    gsl_monte_vegas_free (s);
  }

  gsl_rng_free (r);

  return 0;
}


text:
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>

/* number of data points to fit */
#define N 40

struct data {
  size_t n;
  double * y;
};

int
expb_f (const gsl_vector * x, void *data,
        gsl_vector * f)
{
  size_t n = ((struct data *)data)->n;
  double *y = ((struct data *)data)->y;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);
  double b = gsl_vector_get (x, 2);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp (-lambda * t) + b;
      gsl_vector_set (f, i, Yi - y[i]);
    }

  return GSL_SUCCESS;
}

int
expb_df (const gsl_vector * x, void *data,
         gsl_matrix * J)
{
  size_t n = ((struct data *)data)->n;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double e = exp(-lambda * t);
      gsl_matrix_set (J, i, 0, e);
      gsl_matrix_set (J, i, 1, -t * A * e);
      gsl_matrix_set (J, i, 2, 1.0);
    }
  return GSL_SUCCESS;
}

void
callback(const size_t iter, void *params,
         const gsl_multifit_nlinear_workspace *w)
{
  gsl_vector *f = gsl_multifit_nlinear_residual(w);
  gsl_vector *x = gsl_multifit_nlinear_position(w);
  double rcond;

  /* compute reciprocal condition number of J(x) */
  gsl_multifit_nlinear_rcond(&rcond, w);

  fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n",
          iter,
          gsl_vector_get(x, 0),
          gsl_vector_get(x, 1),
          gsl_vector_get(x, 2),
          1.0 / rcond,
          gsl_blas_dnrm2(f));
}

int
main (void)
{
  const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
  gsl_multifit_nlinear_workspace *w;
  gsl_multifit_nlinear_fdf fdf;
  gsl_multifit_nlinear_parameters fdf_params =
    gsl_multifit_nlinear_default_parameters();
  const size_t n = N;
  const size_t p = 3;

  gsl_vector *f;
  gsl_matrix *J;
  gsl_matrix *covar = gsl_matrix_alloc (p, p);
  double y[N], weights[N];
  struct data d = { n, y };
  double x_init[3] = { 1.0, 1.0, 0.0 }; /* starting values */
  gsl_vector_view x = gsl_vector_view_array (x_init, p);
  gsl_vector_view wts = gsl_vector_view_array(weights, n);
  gsl_rng * r;
  double chisq, chisq0;
  int status, info;
  size_t i;

  const double xtol = 1e-8;
  const double gtol = 1e-8;
  const double ftol = 0.0;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* define the function to be minimized */
  fdf.f = expb_f;
  fdf.df = expb_df;   /* set to NULL for finite-difference Jacobian */
  fdf.fvv = NULL;     /* not using geodesic acceleration */
  fdf.n = n;
  fdf.p = p;
  fdf.params = &d;

  /* this is the data to be fitted */
  for (i = 0; i < n; i++)
    {
      double t = i;
      double yi = 1.0 + 5 * exp (-0.1 * t);
      double si = 0.1 * yi;
      double dy = gsl_ran_gaussian(r, si);

      weights[i] = 1.0 / (si * si);
      y[i] = yi + dy;
      printf ("data: %zu %g %g\n", i, y[i], si);
    };

  /* allocate workspace with default parameters */
  w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);

  /* initialize solver with starting point and weights */
  gsl_multifit_nlinear_winit (&x.vector, &wts.vector, &fdf, w);

  /* compute initial cost function */
  f = gsl_multifit_nlinear_residual(w);
  gsl_blas_ddot(f, f, &chisq0);

  /* solve the system with a maximum of 20 iterations */
  status = gsl_multifit_nlinear_driver(20, xtol, gtol, ftol,
                                       callback, NULL, &info, w);

  /* compute covariance of best fit parameters */
  J = gsl_multifit_nlinear_jac(w);
  gsl_multifit_nlinear_covar (J, 0.0, covar);

  /* compute final cost */
  gsl_blas_ddot(f, f, &chisq);

#define FIT(i) gsl_vector_get(w->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

  fprintf(stderr, "summary from method '%s/%s'\n",
          gsl_multifit_nlinear_name(w),
          gsl_multifit_nlinear_trs_name(w));
  fprintf(stderr, "number of iterations: %zu\n",
          gsl_multifit_nlinear_niter(w));
  fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
  fprintf(stderr, "Jacobian evaluations: %zu\n", fdf.nevaldf);
  fprintf(stderr, "reason for stopping: %s\n",
          (info == 1) ? "small step size" : "small gradient");
  fprintf(stderr, "initial |f(x)| = %f\n", sqrt(chisq0));
  fprintf(stderr, "final   |f(x)| = %f\n", sqrt(chisq));

  {
    double dof = n - p;
    double c = GSL_MAX_DBL(1, sqrt(chisq / dof));

    fprintf(stderr, "chisq/dof = %g\n", chisq / dof);

    fprintf (stderr, "A      = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
    fprintf (stderr, "lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
    fprintf (stderr, "b      = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
  }

  fprintf (stderr, "status = %s\n", gsl_strerror (status));

  gsl_multifit_nlinear_free (w);
  gsl_matrix_free (covar);
  gsl_rng_free (r);

  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE