text:
#include "example.h"
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int
print_gsl_sf_bessel_J0 (double x)
{
  double y = gsl_sf_bessel_J0 (x);
  printf ("J0(%g) = %.18e\n", x, y);
  return 0;
}

ext: c

max_stars_repo_name: Dwoosh/cBinder

max_issues_repo_name: Tetrite/cBinder

max_forks_repo_name: Dwoosh/cBinder


text:
// This file is part of playd.
// playd is licensed under the MIT licence: see LICENSE.txt.

/**
 * @file
 * The Ring_buffer class.
 */

#ifndef PLAYD_RING_BUFFER_HPP
#define PLAYD_RING_BUFFER_HPP

#include <atomic>
#include <mutex>
#include <vector>

#undef max
#include <gsl/gsl>

namespace Playd::Audio
{
/**
 * A concurrent ring buffer.
 *
 * This is not particularly efficient, but does the job for playd.
 * It uses two release-acquire-atomic counters to store read and write.
 */
class RingBuffer
{
public:
	/**
	 * Constructs a Ring_buffer.
	 * @param capacity The capacity of the ring buffer, in bytes.
	 */
	explicit RingBuffer(size_t capacity);

	/// Destructs a Ring_buffer.
	~RingBuffer() = default;

	/// Deleted copy constructor.
	RingBuffer(const RingBuffer &) = delete;

	/// Deleted copy-assignment.
	RingBuffer &operator=(const RingBuffer &) = delete;

	/**
	 * The current write capacity.
	 * @return The number of samples this ring buffer has space to store.
	 * @see Write
	 */
	size_t WriteCapacity() const;

	/**
	 * The current read capacity.
	 * @return The number of samples available in this ring buffer.
	 * @see Read
	 */
	size_t ReadCapacity() const;

	/**
	 * Writes samples from a span into the ring buffer.
	 *
	 * * Precondition: @a src is a valid span.
	 * * Postcondition: The ringbuffer has been written to with the contents
	 *     of the first WriteCapacity() bytes of @a src.
	 *
	 * @param src The span of bytes to write into the ring buffer.
	 * @return The number of bytes written.
	 * @see WriteCapacity
	 */
	size_t Write(gsl::span<const std::byte> src);

	/**
	 * Reads samples from the ring buffer into an array.
	 * To read one sample, pass a count of 1 and take a pointer to the
	 * sample variable.
	 *
	 * * Precondition: @a dest is a valid span.
	 * * Postcondition: The first ReadCapacity() bytes of @a dest have been
	 *     filled with the appropriate number of bytes from the front of the
	 *     ring buffer.
	 *
	 * @param dest The span of bytes to fill with bytes read from the ring
	 *  buffer.
	 * @return The number of bytes read.
	 * @see ReadCapacity
	 */
	size_t Read(gsl::span<std::byte> dest);

	/// Empties the ring buffer.
	void Flush();

private:
	/// Empties the ring buffer without acquiring locks.
	void FlushInner();

	std::vector<std::byte> buffer; ///< The array used by the ringbuffer.

	std::vector<std::byte>::const_iterator r_it; ///< The read iterator.
	std::vector<std::byte>::iterator w_it;       ///< The write iterator.

	std::atomic<size_t> count; ///< The current read capacity.
	// Write capacity is the total buffer capacity minus count.

	std::mutex r_lock; ///< The read lock.
	std::mutex w_lock; ///< The write lock.
};

} // namespace Playd::Audio

#endif // PLAYD_RINGBUFFER_HPP

ext: h

max_stars_repo_name: UniversityRadioYork/ury-playd

max_issues_repo_name: UniversityRadioYork/ury-playd

max_forks_repo_name: UniversityRadioYork/ury-playd


text:
#ifndef BIO_GSL_H_
#define BIO_GSL_H_

#include "bio/defs.h"
#include "bio/raii.h"

#include <gsl/gsl_integration.h>

BIO_NS_START

void
gsl_init();

template <>
RAII<gsl_integration_workspace *>::~RAII();

#define BIO_GSL_EXP(x) ((x) < GSL_LOG_DBL_MIN ? 0.0 : gsl_sf_exp(x))

BIO_NS_END

#endif //BIO_GSL_H_

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
/* linalg/choleskyc.c
 *
 * Copyright (C) 2007 Patrick Alken
 * Copyright (C) 2010 Huan Wu (gsl_linalg_complex_cholesky_invert)
 *
 * 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 <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

/*
 * This module contains routines related to the Cholesky decomposition
 * of a complex Hermitian positive definite matrix.
 */

static void cholesky_complex_conj_vector(gsl_vector_complex *v);

/*
gsl_linalg_complex_cholesky_decomp()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix. See Golub & Van Loan, "Matrix Computations" (3rd ed),
algorithm 4.2.2.

Inputs: A - (input/output) complex postive definite matrix

Return: success or error

The lower triangle of A is overwritten with the Cholesky decomposition
*/

int
gsl_linalg_complex_cholesky_decomp(gsl_matrix_complex *A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t i, j;
      gsl_complex z;
      double ajj;

      for (j = 0; j < N; ++j)
        {
          z = gsl_matrix_complex_get(A, j, j);
          ajj = GSL_REAL(z);

          if (j > 0)
            {
              gsl_vector_complex_const_view aj =
                gsl_matrix_complex_const_subrow(A, j, 0, j);

              gsl_blas_zdotc(&aj.vector, &aj.vector, &z);
              ajj -= GSL_REAL(z);
            }

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(A, j, j, z);

          if (j < N - 1)
            {
              gsl_vector_complex_view av =
                gsl_matrix_complex_subcolumn(A, j, j + 1, N - j - 1);

              if (j > 0)
                {
                  gsl_vector_complex_view aj =
                    gsl_matrix_complex_subrow(A, j, 0, j);
                  gsl_matrix_complex_view am =
                    gsl_matrix_complex_submatrix(A, j + 1, 0, N - j - 1, j);

                  cholesky_complex_conj_vector(&aj.vector);

                  gsl_blas_zgemv(CblasNoTrans,
                                 GSL_COMPLEX_NEGONE,
                                 &am.matrix,
                                 &aj.vector,
                                 GSL_COMPLEX_ONE,
                                 &av.vector);

                  cholesky_complex_conj_vector(&aj.vector);
                }

              gsl_blas_zdscal(1.0 / ajj, &av.vector);
            }
        }

      /* Now store L^H in upper triangle */
      for (i = 1; i < N; ++i)
        {
          for (j = 0; j < i; ++j)
            {
              z = gsl_matrix_complex_get(A, i, j);
              gsl_matrix_complex_set(A, j, i, gsl_complex_conjugate(z));
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_decomp() */

/*
gsl_linalg_complex_cholesky_solve()
  Solve A x = b where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                   const gsl_vector_complex * b,
                                   gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_complex_memcpy (x, b);

      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_solve() */

/*
gsl_linalg_complex_cholesky_svx()
  Solve A x = b in place where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                 gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_svx() */

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

gsl_linalg_complex_cholesky_invert()
  Compute the inverse of an Hermitian positive definite matrix in
  Cholesky form.

Inputs: LLT - matrix in cholesky form on input
              A^{-1} = L^{-H} L^{-1} on output

Return: success or error
******************************************************************************/

int
gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * LLT)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      size_t N = LLT->size1;
      size_t i, j;
      gsl_vector_complex_view v1;

      /* invert the lower triangle of LLT */
      for (i = 0; i < N; ++i)
        {
          double ajj;
          gsl_complex z;

          j = N - i - 1;

          {
            gsl_complex z0 = gsl_matrix_complex_get(LLT, j, j);
            ajj = 1.0 / GSL_REAL(z0);
          }

          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(LLT, j, j, z);

          {
            gsl_complex z1 = gsl_matrix_complex_get(LLT, j, j);
            ajj = -GSL_REAL(z1);
          }

          if (j < N - 1)
            {
              gsl_matrix_complex_view m;

              m = gsl_matrix_complex_submatrix(LLT, j + 1, j + 1,
                                       N - j - 1, N - j - 1);
              v1 = gsl_matrix_complex_subcolumn(LLT, j, j + 1, N - j - 1);

              gsl_blas_ztrmv(CblasLower, CblasNoTrans, CblasNonUnit,
                             &m.matrix, &v1.vector);

              gsl_blas_zdscal(ajj, &v1.vector);
            }
        } /* for (i = 0; i < N; ++i) */

      /*
       * The lower triangle of LLT now contains L^{-1}. Now compute
       * A^{-1} = L^{-H} L^{-1}
       *
       * The (ij) element of A^{-1} is column i of conj(L^{-1}) dotted into
       * column j of L^{-1}
       */

      for (i = 0; i < N; ++i)
        {
          gsl_complex sum;
          for (j = i + 1; j < N; ++j)
            {
              gsl_vector_complex_view v2;
              v1 = gsl_matrix_complex_subcolumn(LLT, i, j, N - j);
              v2 = gsl_matrix_complex_subcolumn(LLT, j, j, N - j);

              /* compute Ainv[i,j] = sum_k{conj(Linv[k,i]) * Linv[k,j]} */
              gsl_blas_zdotc(&v1.vector, &v2.vector, &sum);

              /* store in upper triangle */
              gsl_matrix_complex_set(LLT, i, j, sum);
            }

          /* now compute the diagonal element */
          v1 = gsl_matrix_complex_subcolumn(LLT, i, i, N - i);
          gsl_blas_zdotc(&v1.vector, &v1.vector, &sum);
          gsl_matrix_complex_set(LLT, i, i, sum);
        }

      /* copy the Hermitian upper triangle to the lower triangle */

      for (j = 1; j < N; j++)
        {
          for (i = 0; i < j; i++)
            {
              gsl_complex z = gsl_matrix_complex_get(LLT, i, j);
              gsl_matrix_complex_set(LLT, j, i, gsl_complex_conjugate(z));
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_invert() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

static void
cholesky_complex_conj_vector(gsl_vector_complex *v)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(v, i);
      gsl_vector_complex_set(v, i, gsl_complex_conjugate(z));
    }
} /* cholesky_complex_conj_vector() */

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <nlopt.h>
#include <mex.h>

struct mexObjectiveFunction
{
  static double obj_fun(unsigned n, const double *x, double *gradient, void *d_);
  static void precond_fun(unsigned n, const double *x, const double *v, double *vpre, void *f_data);

  mxArray *prhs[2];          // feval mexMatlabCall input arguments for objective function evaluation
  nlopt_opt &opt;
  mxArray *hessmult_args[3];   // Hv = HessMultFcn(x,v)
  mxArray *outfun_args[4];   // feval mexMatlabCall input arguments for OutputFun function evaluation
  mxArray *lasterror;        // trapped MException
  bool stop;                 // true if user issued a stop

  mexObjectiveFunction(nlopt_opt & optin, mxArray * mxFun, mxArray * mxHessMultFcn, mxArray * mxOutputFun);
  ~mexObjectiveFunction();

  double evalFun(unsigned n, const double *x, double *gradient);
  void evalHessMultFcn(unsigned n, const double *x, const double *v, double *vpre);
  bool evalOutputFun(bool init);
  mxArray *evalGrad(mxArray *x);
private:
  static mxArray *create_optimvalues();
  void set_outfun_args(const char *state, mxArray *mxFval, mxArray *mxGradient = NULL);
  bool call_matlab_feval_with_trap(int nlhs, mxArray *plhs[], int nrhs, mxArray *prhs[]);
};

max_stars_repo_name: hokiedsp/matlab-nlopt

max_issues_repo_name: hokiedsp/matlab-nlopt

max_forks_repo_name: hokiedsp/matlab-nlopt


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include "Effective_Potential.h"

#define PI 3.14159265
#define G 4.30091252525*pow(10,-3) //grav constant, in (parsec*km^2)/(Ms*sec^2)
#define c  300000 //speed of light, km/sec

void Geodesic_Orbit(int* time_size, double M, double periapsis, double e)
{
    int dimension = 2;		/* number of differential equations */

    double eps_abs = 1.e-8;	/* absolute error requested */
    double eps_rel = 1.e-10;	/* relative error requested */

    /* define the type of routine for making steps: */
    const gsl_odeiv_step_type *type_ptr = gsl_odeiv_step_rkf45;

    /*
       allocate/initialize the stepper, the control function, and the
       evolution function.
    */
    gsl_odeiv_step *step_ptr = gsl_odeiv_step_alloc (type_ptr, dimension);
    gsl_odeiv_control *control_ptr = gsl_odeiv_control_y_new (eps_abs, eps_rel);
    gsl_odeiv_evolve *evolve_ptr = gsl_odeiv_evolve_alloc (dimension);

    gsl_odeiv_system my_system;	/* structure with the rhs function, etc. */

    double E = E_r(M, periapsis, e);
    double L = L_r(e,periapsis,M);
    double r[*time_size];			/* current solution vector */

    double T, T_next;		/* current and next independent variable */
    double Tmin, Tmax, delta_t;	/* range of t and step size for output */

    double h = 1e-6;		/* starting step size for ode solver */
}

max_stars_repo_name: ladanuzhna/Black-Holes-simulations

max_issues_repo_name: ladanuzhna/Black-Holes-simulations

max_forks_repo_name: ladanuzhna/Black-Holes-simulations


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

/* Author:  G. Jungman */

/* Evaluation of Coulomb wave functions F_L(eta, x), G_L(eta, x),
 * and their derivatives. A combination of Steed's method, asymptotic
 * results, and power series.
 *
 * Steed's method:
 *  [Barnett, CPC 21, 297 (1981)]
 * Power series and other methods:
 *  [Biedenharn et al., PR 97, 542 (1954)]
 *  [Bardin et al., CPC 3, 73 (1972)]
 *  [Abad+Sesma, CPC 71, 110 (1992)]
 */
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_coulomb.h>

#include "error.h"

/* the L=0 normalization constant
 * [Abramowitz+Stegun 14.1.8]
 */
static
double
C0sq(double eta)
{
  double twopieta = 2.0*M_PI*eta;

  if(fabs(eta) < GSL_DBL_EPSILON) {
    return 1.0;
  }
  else if(twopieta > GSL_LOG_DBL_MAX) {
    return 0.0;
  }
  else {
    gsl_sf_result scale;
    gsl_sf_expm1_e(twopieta, &scale);
    return twopieta/scale.val;
  }
}

/* the full definition of C_L(eta) for any valid L and eta
 * [Abramowitz and Stegun 14.1.7]
 * This depends on the complex gamma function. For large
 * arguments the phase of the complex gamma function is not
 * very accurately determined. However the modulus is, and that
 * is all that we need to calculate C_L.
 *
 * This is not valid for L <= -3/2  or  L = -1.
 */
static
int
CLeta(double L, double eta, gsl_sf_result * result)
{
  gsl_sf_result ln1; /* log of numerator Gamma function */
  gsl_sf_result ln2; /* log of denominator Gamma function */
  double sgn = 1.0;
  double arg_val, arg_err;

  if(fabs(eta/(L+1.0)) < GSL_DBL_EPSILON) {
    gsl_sf_lngamma_e(L+1.0, &ln1);
  }
  else {
    gsl_sf_result p1;                 /* phase of numerator Gamma -- not used */
    gsl_sf_lngamma_complex_e(L+1.0, eta, &ln1, &p1); /* should be ok */
  }

  gsl_sf_lngamma_e(2.0*(L+1.0), &ln2);
  if(L < -1.0) sgn = -sgn;

  arg_val  = L*M_LN2 - 0.5*eta*M_PI + ln1.val - ln2.val;
  arg_err  = ln1.err + ln2.err;
  arg_err += GSL_DBL_EPSILON * (fabs(L*M_LN2) + fabs(0.5*eta*M_PI));
  return gsl_sf_exp_err_e(arg_val, arg_err, result);
}

int
gsl_sf_coulomb_CL_e(double lam, double eta, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(lam <= -1.0) {
    DOMAIN_ERROR(result);
  }
  else if(fabs(lam) < GSL_DBL_EPSILON) {
    /* saves a calculation of complex_lngamma(), otherwise not necessary */
    result->val = sqrt(C0sq(eta));
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    return CLeta(lam, eta, result);
  }
}

/* cl[0] .. cl[kmax] = C_{lam_min}(eta) .. C_{lam_min+kmax}(eta)
 */
int
gsl_sf_coulomb_CL_array(double lam_min, int kmax, double eta, double * cl)
{
  int k;
  gsl_sf_result cl_0;
  gsl_sf_coulomb_CL_e(lam_min, eta, &cl_0);
  cl[0] = cl_0.val;

  for(k=1; k<=kmax; k++) {
    double L = lam_min + k;
    cl[k] = cl[k-1] * hypot(L, eta)/(L*(2.0*L+1.0));
  }

  return GSL_SUCCESS;
}

/* Evaluate the series for Phi_L(eta,x) and Phi_L*(eta,x)
 * [Abramowitz+Stegun 14.1.5]
 * [Abramowitz+Stegun 14.1.13]
 *
 * The sequence of coefficients A_k^L is
 * manifestly well-controlled for L >= -1/2
 * and eta < 10.
 *
 * This makes sense since this is the region
 * away from threshold, and you expect
 * the evaluation to become easier as you
 * get farther from threshold.
 *
 * Empirically, this is quite well-behaved for
 *   L >= -1/2
 *   eta < 10
 *   x   < 10
 */
#if 0
static
int
coulomb_Phi_series(const double lam, const double eta, const double x,
                   double * result, double * result_star)
{
  int kmin =   5;
  int kmax = 200;
  int k;
  double Akm2 = 1.0;
  double Akm1 = eta/(lam+1.0);
  double Ak;

  double xpow = x;
  double sum  = Akm2 + Akm1*x;
  double sump = (lam+1.0)*Akm2 + (lam+2.0)*Akm1*x;
  double prev_abs_del   = fabs(Akm1*x);
  double prev_abs_del_p = (lam+2.0) * prev_abs_del;

  for(k=2; k<kmax; k++) {
    double del;
    double del_p;
    double abs_del;
    double abs_del_p;

    Ak = (2.0*eta*Akm1 - Akm2)/(k*(2.0*lam + 1.0 + k));

    xpow *= x;
    del   = Ak*xpow;
    del_p = (k+lam+1.0)*del;
    sum  += del;
    sump += del_p;

    abs_del   = fabs(del);
    abs_del_p = fabs(del_p);

    if(          abs_del/(fabs(sum)+abs_del)          < GSL_DBL_EPSILON
       &&   prev_abs_del/(fabs(sum)+prev_abs_del)     < GSL_DBL_EPSILON
       &&      abs_del_p/(fabs(sump)+abs_del_p)       < GSL_DBL_EPSILON
       && prev_abs_del_p/(fabs(sump)+prev_abs_del_p)  < GSL_DBL_EPSILON
       && k > kmin
       ) break;

    /* We need to keep track of the previous delta because when
     * eta is near zero the odd terms of the sum are very small
     * and this could lead to premature termination.
     */
    prev_abs_del   = abs_del;
    prev_abs_del_p = abs_del_p;

    Akm2 = Akm1;
    Akm1 = Ak;
  }

  *result      = sum;
  *result_star = sump;

  if(k==kmax) {
    GSL_ERROR ("error", GSL_EMAXITER);
  }
  else {
    return GSL_SUCCESS;
  }
}
#endif /* 0 */

/* Determine the connection phase, phi_lambda.
 * See coulomb_FG_series() below. We have
 * to be careful about sin(phi)->0. Note that
 * there is an underflow condition for large
 * positive eta in any case.
 */
static
int
coulomb_connection(const double lam, const double eta,
                   double * cos_phi, double * sin_phi)
{
  if(eta > -GSL_LOG_DBL_MIN/2.0*M_PI-1.0) {
    *cos_phi = 1.0;
    *sin_phi = 0.0;
    GSL_ERROR ("error", GSL_EUNDRFLW);
  }
  else if(eta > -GSL_LOG_DBL_EPSILON/(4.0*M_PI)) {
    const double eps = 2.0 * exp(-2.0*M_PI*eta);
    const double tpl = tan(M_PI * lam);
    const double dth = eps * tpl / (tpl*tpl + 1.0);
    *cos_phi = -1.0 + 0.5 * dth*dth;
    *sin_phi = -dth;
    return GSL_SUCCESS;
  }
  else {
    double X   = tanh(M_PI * eta) / tan(M_PI * lam);
    double phi = -atan(X) - (lam + 0.5) * M_PI;
    *cos_phi = cos(phi);
    *sin_phi = sin(phi);
    return GSL_SUCCESS;
  }
}

/* Evaluate the Frobenius series for F_lam(eta,x) and G_lam(eta,x).
 * Homegrown algebra. Evaluates the series for F_{lam} and
 * F_{-lam-1}, then uses
 *    G_{lam} = (F_{lam} cos(phi) - F_{-lam-1}) / sin(phi)
 * where
 *    phi = Arg[Gamma[1+lam+I eta]] - Arg[Gamma[-lam + I eta]] - (lam+1/2)Pi
 *        = Arg[Sin[Pi(-lam+I eta)] - (lam+1/2)Pi
 *        = atan2(-cos(lam Pi)sinh(eta Pi), -sin(lam Pi)cosh(eta Pi)) - (lam+1/2)Pi
 *
 *        = -atan(X) - (lam+1/2) Pi,  X = tanh(eta Pi)/tan(lam Pi)
 *
 * Not appropriate for lam <= -1/2, lam = 0, or lam >= 1/2.
 */
static
int
coulomb_FG_series(const double lam, const double eta, const double x,
                  gsl_sf_result * F, gsl_sf_result * G)
{
  const int max_iter = 800;
  gsl_sf_result ClamA;
  gsl_sf_result ClamB;
  int stat_A = CLeta(lam, eta, &ClamA);
  int stat_B = CLeta(-lam-1.0, eta, &ClamB);
  const double tlp1 = 2.0*lam + 1.0;
  const double pow_x = pow(x, lam);
  double cos_phi_lam;
  double sin_phi_lam;

  double uA_mm2 = 1.0;                  /* uA sum is for F_{lam} */
  double uA_mm1 = x*eta/(lam+1.0);
  double uA_m;
  double uB_mm2 = 1.0;                  /* uB sum is for F_{-lam-1} */
  double uB_mm1 = -x*eta/lam;
  double uB_m;
  double A_sum = uA_mm2 + uA_mm1;
  double B_sum = uB_mm2 + uB_mm1;
  double A_abs_del_prev = fabs(A_sum);
  double B_abs_del_prev = fabs(B_sum);
  gsl_sf_result FA, FB;
  int m = 2;

  int stat_conn = coulomb_connection(lam, eta, &cos_phi_lam, &sin_phi_lam);

  if(stat_conn == GSL_EUNDRFLW) {
    F->val = 0.0;  /* FIXME: should this be set to Inf too like G? */
    F->err = 0.0;
    OVERFLOW_ERROR(G);
  }

  while(m < max_iter) {
    double abs_dA;
    double abs_dB;
    uA_m = x*(2.0*eta*uA_mm1 - x*uA_mm2)/(m*(m+tlp1));
    uB_m = x*(2.0*eta*uB_mm1 - x*uB_mm2)/(m*(m-tlp1));
    A_sum += uA_m;
    B_sum += uB_m;
    abs_dA = fabs(uA_m);
    abs_dB = fabs(uB_m);
    if(m > 15) {
      /* Don't bother checking until we have gone out a little ways;
       * a minor optimization. Also make sure to check both the
       * current and the previous increment because the odd and even
       * terms of the sum can have very different behaviour, depending
       * on the value of eta.
       */
      double max_abs_dA = GSL_MAX(abs_dA, A_abs_del_prev);
      double max_abs_dB = GSL_MAX(abs_dB, B_abs_del_prev);
      double abs_A = fabs(A_sum);
      double abs_B = fabs(B_sum);
      if(   max_abs_dA/(max_abs_dA + abs_A) < 4.0*GSL_DBL_EPSILON
         && max_abs_dB/(max_abs_dB + abs_B) < 4.0*GSL_DBL_EPSILON
         ) break;
    }
    A_abs_del_prev = abs_dA;
    B_abs_del_prev = abs_dB;
    uA_mm2 = uA_mm1;
    uA_mm1 = uA_m;
    uB_mm2 = uB_mm1;
    uB_mm1 = uB_m;
    m++;
  }

  FA.val = A_sum * ClamA.val * pow_x * x;
  FA.err = fabs(A_sum) * ClamA.err * pow_x * x + 2.0*GSL_DBL_EPSILON*fabs(FA.val);
  FB.val = B_sum * ClamB.val / pow_x;
  FB.err = fabs(B_sum) * ClamB.err / pow_x + 2.0*GSL_DBL_EPSILON*fabs(FB.val);

  F->val = FA.val;
  F->err = FA.err;

  G->val = (FA.val * cos_phi_lam - FB.val)/sin_phi_lam;
  G->err = (FA.err * fabs(cos_phi_lam) + FB.err)/fabs(sin_phi_lam);

  if(m >= max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_ERROR_SELECT_2(stat_A, stat_B);
}

/* Evaluate the Frobenius series for F_0(eta,x) and G_0(eta,x).
 * See [Bardin et al., CPC 3, 73 (1972), (14)-(17)];
 * note the misprint in (17): nu_0=1 is correct, not nu_0=0.
 */
static
int
coulomb_FG0_series(const double eta, const double x,
                   gsl_sf_result * F, gsl_sf_result * G)
{
  const int max_iter = 800;
  const double x2  = x*x;
  const double tex = 2.0*eta*x;
  gsl_sf_result C0;
  int stat_CL = CLeta(0.0, eta, &C0);
  gsl_sf_result r1pie;
  int psi_stat = gsl_sf_psi_1piy_e(eta, &r1pie);
  double u_mm2 = 0.0;  /* u_0 */
  double u_mm1 = x;    /* u_1 */
  double u_m;
  double v_mm2 = 1.0;                               /* nu_0 */
  double v_mm1 = tex*(2.0*M_EULER-1.0+r1pie.val);   /* nu_1 */
  double v_m;
  double u_sum = u_mm2 + u_mm1;
  double v_sum = v_mm2 + v_mm1;
  double u_abs_del_prev = fabs(u_sum);
  double v_abs_del_prev = fabs(v_sum);
  int m = 2;
  double u_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(u_sum);
  double v_sum_err = 2.0 * GSL_DBL_EPSILON * fabs(v_sum);
  double ln2x = log(2.0*x);

  while(m < max_iter) {
    double abs_du;
    double abs_dv;
    double m_mm1 = m*(m-1.0);
    u_m = (tex*u_mm1 - x2*u_mm2)/m_mm1;
    v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*eta*(2*m-1)*u_m)/m_mm1;
    u_sum += u_m;
    v_sum += v_m;
    abs_du = fabs(u_m);
    abs_dv = fabs(v_m);
    u_sum_err += 2.0 * GSL_DBL_EPSILON * abs_du;
    v_sum_err += 2.0 * GSL_DBL_EPSILON * abs_dv;
    if(m > 15) {
      /* Don't bother checking until we have gone out a little ways;
       * a minor optimization. Also make sure to check both the
       * current and the previous increment because the odd and even
       * terms of the sum can have very different behaviour, depending
       * on the value of eta.
       */
      double max_abs_du = GSL_MAX(abs_du, u_abs_del_prev);
      double max_abs_dv = GSL_MAX(abs_dv, v_abs_del_prev);
      double abs_u = fabs(u_sum);
      double abs_v = fabs(v_sum);
      if(   max_abs_du/(max_abs_du + abs_u) < 40.0*GSL_DBL_EPSILON
         && max_abs_dv/(max_abs_dv + abs_v) < 40.0*GSL_DBL_EPSILON
         ) break;
    }
    u_abs_del_prev = abs_du;
    v_abs_del_prev = abs_dv;
    u_mm2 = u_mm1;
    u_mm1 = u_m;
    v_mm2 = v_mm1;
    v_mm1 = v_m;
    m++;
  }

  F->val  = C0.val * u_sum;
  F->err  = C0.err * fabs(u_sum);
  F->err += fabs(C0.val) * u_sum_err;
  F->err += 2.0 * GSL_DBL_EPSILON * fabs(F->val);

  G->val  = (v_sum + 2.0*eta*u_sum * ln2x) / C0.val;
  G->err  = (fabs(v_sum) + fabs(2.0*eta*u_sum * ln2x)) / fabs(C0.val) * fabs(C0.err/C0.val);
  G->err += (v_sum_err + fabs(2.0*eta*u_sum_err*ln2x)) / fabs(C0.val);
  G->err += 2.0 * GSL_DBL_EPSILON * fabs(G->val);

  if(m == max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_ERROR_SELECT_2(psi_stat, stat_CL);
}

/* Evaluate the Frobenius series for F_{-1/2}(eta,x) and G_{-1/2}(eta,x).
 * Homegrown algebra.
 */
static
int
coulomb_FGmhalf_series(const double eta, const double x,
                       gsl_sf_result * F, gsl_sf_result * G)
{
  const int max_iter = 800;
  const double rx  = sqrt(x);
  const double x2  = x*x;
  const double tex = 2.0*eta*x;
  gsl_sf_result Cmhalf;
  int stat_CL = CLeta(-0.5, eta, &Cmhalf);
  double u_mm2 = 1.0;                      /* u_0 */
  double u_mm1 = tex * u_mm2;              /* u_1 */
  double u_m;
  double v_mm2, v_mm1, v_m;
  double f_sum, g_sum;
  double tmp1;
  gsl_sf_result rpsi_1pe;
  gsl_sf_result rpsi_1p2e;
  int m = 2;

  gsl_sf_psi_1piy_e(eta,     &rpsi_1pe);
  gsl_sf_psi_1piy_e(2.0*eta, &rpsi_1p2e);

  v_mm2 = 2.0*M_EULER - M_LN2 - rpsi_1pe.val + 2.0*rpsi_1p2e.val;
  v_mm1 = tex*(v_mm2 - 2.0*u_mm2);

  f_sum = u_mm2 + u_mm1;
  g_sum = v_mm2 + v_mm1;

  while(m < max_iter) {
    double m2 = m*m;
    u_m = (tex*u_mm1 - x2*u_mm2)/m2;
    v_m = (tex*v_mm1 - x2*v_mm2 - 2.0*m*u_m)/m2;
    f_sum += u_m;
    g_sum += v_m;
    if(   f_sum != 0.0
       && g_sum != 0.0
       && (fabs(u_m/f_sum) + fabs(v_m/g_sum) < 10.0*GSL_DBL_EPSILON)) break;
    u_mm2 = u_mm1;
    u_mm1 = u_m;
    v_mm2 = v_mm1;
    v_mm1 = v_m;
    m++;
  }

  F->val = Cmhalf.val * rx * f_sum;
  F->err = Cmhalf.err * fabs(rx * f_sum) + 2.0*GSL_DBL_EPSILON*fabs(F->val);

  tmp1 = f_sum*log(x);
  G->val = -rx*(tmp1 + g_sum)/Cmhalf.val;
  G->err = fabs(rx)*(fabs(tmp1) + fabs(g_sum))/fabs(Cmhalf.val) * fabs(Cmhalf.err/Cmhalf.val);

  if(m == max_iter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return stat_CL;
}

/* Evolve the backwards recurrence for F,F'.
 *
 *    F_{lam-1}  = (S_lam F_lam + F_lam') / R_lam
 *    F_{lam-1}' = (S_lam F_{lam-1} - R_lam F_lam)
 * where
 *    R_lam = sqrt(1 + (eta/lam)^2)
 *    S_lam = lam/x + eta/lam
 *
 */
static
int
coulomb_F_recur(double lam_min, int kmax,
                double eta, double x,
                double F_lam_max, double Fp_lam_max,
                double * F_lam_min, double * Fp_lam_min
                )
{
  double x_inv = 1.0/x;
  double fcl = F_lam_max;
  double fpl = Fp_lam_max;
  double lam_max = lam_min + kmax;
  double lam = lam_max;
  int k;

  for(k=kmax-1; k>=0; k--) {
    double el = eta/lam;
    double rl = hypot(1.0, el);
    double sl = el  + lam*x_inv;
    double fc_lm1;
    fc_lm1 = (fcl*sl + fpl)/rl;
    fpl    =  fc_lm1*sl - fcl*rl;
    fcl    =  fc_lm1;
    lam -= 1.0;
  }

  *F_lam_min  = fcl;
  *Fp_lam_min = fpl;
  return GSL_SUCCESS;
}

/* Evolve the forward recurrence for G,G'.
 *
 *   G_{lam+1}  = (S_lam G_lam - G_lam')/R_lam
 *   G_{lam+1}' = R_{lam+1} G_lam - S_lam G_{lam+1}
 *
 * where S_lam and R_lam are as above in the F recursion.
 */
static
int
coulomb_G_recur(const double lam_min, const int kmax,
                const double eta, const double x,
                const double G_lam_min, const double Gp_lam_min,
                double * G_lam_max, double * Gp_lam_max
                )
{
  double x_inv = 1.0/x;
  double gcl = G_lam_min;
  double gpl = Gp_lam_min;
  double lam = lam_min + 1.0;
  int k;

  for(k=1; k<=kmax; k++) {
    double el = eta/lam;
    double rl = hypot(1.0, el);
    double sl = el + lam*x_inv;
    double gcl1 = (sl*gcl - gpl)/rl;
    gpl   = rl*gcl - sl*gcl1;
    gcl   = gcl1;
    lam += 1.0;
  }

  *G_lam_max  = gcl;
  *Gp_lam_max = gpl;
  return GSL_SUCCESS;
}

/* Evaluate the first continued fraction, giving
 * the ratio F'/F at the upper lambda value.
 * We also determine the sign of F at that point,
 * since it is the sign of the last denominator
 * in the continued fraction.
 */
static
int
coulomb_CF1(double lambda,
            double eta, double x,
            double * fcl_sign,
            double * result,
            int * count
            )
{
  const double CF1_small = 1.e-30;
  const double CF1_abort = 1.0e+05;
  const double CF1_acc   = 2.0*GSL_DBL_EPSILON;
  const double x_inv     = 1.0/x;
  const double px        = lambda + 1.0 + CF1_abort;

  double pk = lambda + 1.0;
  double F  = eta/pk + pk*x_inv;
  double D, C;
  double df;

  *fcl_sign = 1.0;
  *count = 0;

  if(fabs(F) < CF1_small) F = CF1_small;
  D = 0.0;
  C = F;

  do {
    double pk1 = pk + 1.0;
    double ek  = eta / pk;
    double rk2 = 1.0 + ek*ek;
    double tk  = (pk + pk1)*(x_inv + ek/pk1);
    D   =  tk - rk2 * D;
    C   =  tk - rk2 / C;
    if(fabs(C) < CF1_small) C = CF1_small;
    if(fabs(D) < CF1_small) D = CF1_small;
    D = 1.0/D;
    df = D * C;
    F  = F * df;
    if(D < 0.0) {
      /* sign of result depends on sign of denominator */
      *fcl_sign = - *fcl_sign;
    }
    pk = pk1;
    if( pk > px ) {
      *result = F;
      GSL_ERROR ("error", GSL_ERUNAWAY);
    }
    ++(*count);
  }
  while(fabs(df-1.0) > CF1_acc);

  *result = F;
  return GSL_SUCCESS;
}

#if 0
static
int
old_coulomb_CF1(const double lambda,
                double eta, double x,
                double * fcl_sign,
                double * result
                )
{
  const double CF1_abort = 1.e5;
  const double CF1_acc   = 10.0*GSL_DBL_EPSILON;
  const double x_inv     = 1.0/x;
  const double px        = lambda + 1.0 + CF1_abort;

  double pk = lambda + 1.0;

  double D;
  double df;

  double F;
  double p;
  double pk1;
  double ek;

  double fcl = 1.0;

  double tk;

  while(1) {
    ek = eta/pk;
    F = (ek + pk*x_inv)*fcl + (fcl - 1.0)*x_inv;
    pk1 = pk + 1.0;
    if(fabs(eta*x + pk*pk1) > CF1_acc) break;
    fcl = (1.0 + ek*ek)/(1.0 + eta*eta/(pk1*pk1));
    pk = 2.0 + pk;
  }

  D  = 1.0/((pk + pk1)*(x_inv + ek/pk1));
  df = -fcl*(1.0 + ek*ek)*D;

  if(fcl != 1.0) fcl = -1.0;
  if(D    < 0.0) fcl = -fcl;

  F = F + df;

  p = 1.0;
  do {
    pk = pk1;
    pk1 = pk + 1.0;
    ek  = eta / pk;
    tk  = (pk + pk1)*(x_inv + ek/pk1);
    D   =  tk - D*(1.0+ek*ek);
    if(fabs(D) < sqrt(CF1_acc)) {
      p += 1.0;
      if(p > 2.0) {
        printf("HELP............\n");
      }
    }
    D = 1.0/D;
    if(D < 0.0) {
      /* sign of result depends on sign of denominator */
      fcl = -fcl;
    }
    df = df*(D*tk - 1.0);
    F  = F + df;
    if( pk > px ) {
      GSL_ERROR ("error", GSL_ERUNAWAY);
    }
  }
  while(fabs(df) > fabs(F)*CF1_acc);

  *fcl_sign = fcl;
  *result = F;
  return GSL_SUCCESS;
}
#endif /* 0 */

/* Evaluate the second continued fraction to
 * obtain the ratio
 *    (G' + i F')/(G + i F) := P + i Q
 * at the specified lambda value.
 */
static
int
coulomb_CF2(const double lambda, const double eta, const double x,
            double * result_P, double * result_Q, int * count
            )
{
  int status = GSL_SUCCESS;

  const double CF2_acc   = 4.0*GSL_DBL_EPSILON;
  const double CF2_abort = 2.0e+05;

  const double wi    = 2.0*eta;
  const double x_inv = 1.0/x;
  const double e2mm1 = eta*eta + lambda*(lambda + 1.0);

  double ar = -e2mm1;
  double ai =  eta;

  double br =  2.0*(x - eta);
  double bi =  2.0;

  double dr =  br/(br*br + bi*bi);
  double di = -bi/(br*br + bi*bi);

  double dp = -x_inv*(ar*di + ai*dr);
  double dq =  x_inv*(ar*dr - ai*di);

  double A, B, C, D;

  double pk =  0.0;
  double P  =  0.0;
  double Q  =  1.0 - eta*x_inv;

  *count = 0;

  do {
    P += dp;
    Q += dq;
    pk += 2.0;
    ar += pk;
    ai += wi;
    bi += 2.0;
    D  = ar*dr - ai*di + br;
    di = ai*dr + ar*di + bi;
    C  = 1.0/(D*D + di*di);
    dr =  C*D;
    di = -C*di;
    A  = br*dr - bi*di - 1.;
    B  = bi*dr + br*di;
    C  = dp*A  - dq*B;
    dq = dp*B  + dq*A;
    dp = C;
    if(pk > CF2_abort) {
      status = GSL_ERUNAWAY;
      break;
    }
    ++(*count);
  }
  while(fabs(dp)+fabs(dq) > (fabs(P)+fabs(Q))*CF2_acc);

  if(Q < CF2_abort*GSL_DBL_EPSILON*fabs(P)) {
    status = GSL_ELOSS;
  }

  *result_P = P;
  *result_Q = Q;
  return status;
}

/* WKB evaluation of F, G. Assumes  0 < x < turning point.
 * Overflows are trapped, GSL_EOVRFLW is signalled,
 * and an exponent is returned such that:
 *
 *   result_F = fjwkb * exp(-exponent)
 *   result_G = gjwkb * exp( exponent)
 *
 * See [Biedenharn et al. Phys. Rev. 97, 542-554 (1955), Section IV]
 *
 * Unfortunately, this is not very accurate in general. The
 * test cases typically have 3-4 digits of precision. One could
 * argue that this is ok for general use because, for instance,
 * F is exponentially small in this region and so the absolute
 * accuracy is still roughly acceptable. But it would be better
 * to have a systematic method for improving the precision. See
 * the Abad+Sesma method discussion below.
 */
static
int
coulomb_jwkb(const double lam, const double eta, const double x,
             gsl_sf_result * fjwkb, gsl_sf_result * gjwkb,
             double * exponent)
{
  const double llp1      = lam*(lam+1.0) + 6.0/35.0;
  const double llp1_eff  = GSL_MAX(llp1, 0.0);
  const double rho_ghalf = sqrt(x*(2.0*eta - x) + llp1_eff);
  const double sinh_arg  = sqrt(llp1_eff/(eta*eta+llp1_eff)) * rho_ghalf / x;
  const double sinh_inv  = log(sinh_arg + hypot(1.0,sinh_arg));

  const double phi = fabs(rho_ghalf - eta*atan2(rho_ghalf,x-eta) - sqrt(llp1_eff) * sinh_inv);

  const double zeta_half = pow(3.0*phi/2.0, 1.0/3.0);
  const double prefactor = sqrt(M_PI*phi*x/(6.0 * rho_ghalf));

  double F = prefactor * 3.0/zeta_half;
  double G = prefactor * 3.0/zeta_half; /* Note the sqrt(3) from Bi normalization */
  double F_exp;
  double G_exp;

  const double airy_scale_exp = phi;
  gsl_sf_result ai;
  gsl_sf_result bi;
  gsl_sf_airy_Ai_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &ai);
  gsl_sf_airy_Bi_scaled_e(zeta_half*zeta_half, GSL_MODE_DEFAULT, &bi);
  F *= ai.val;
  G *= bi.val;
  F_exp = log(F) - airy_scale_exp;
  G_exp = log(G) + airy_scale_exp;

  if(G_exp >= GSL_LOG_DBL_MAX) {
    fjwkb->val = F;
    gjwkb->val = G;
    fjwkb->err = 1.0e-3 * fabs(F); /* FIXME: real error here ... could be smaller */
    gjwkb->err = 1.0e-3 * fabs(G);
    *exponent = airy_scale_exp;
    GSL_ERROR ("error", GSL_EOVRFLW);
  }
  else {
    fjwkb->val = exp(F_exp);
    gjwkb->val = exp(G_exp);
    fjwkb->err = 1.0e-3 * fabs(fjwkb->val);
    gjwkb->err = 1.0e-3 * fabs(gjwkb->val);
    *exponent = 0.0;
    return GSL_SUCCESS;
  }
}

/* Asymptotic evaluation of F and G below the minimal turning point.
 *
 * This is meant to be a drop-in replacement for coulomb_jwkb().
 * It uses the expressions in [Abad+Sesma]. This requires some
 * work because I am not sure where it is valid. They mumble
 * something about |x| < |lam|^(-1/2) or 8|eta x| > lam when |x| < 1.
 * This seems true, but I thought the result was based on a uniform
 * expansion and could be controlled by simply using more terms.
 */
#if 0
static
int
coulomb_AS_xlt2eta(const double lam, const double eta, const double x,
                   gsl_sf_result * f_AS, gsl_sf_result * g_AS,
                   double * exponent)
{
  /* no time to do this now... */
}
#endif /* 0 */

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

int
gsl_sf_coulomb_wave_FG_e(const double eta, const double x,
                            const double lam_F,
                            const int  k_lam_G,      /* lam_G = lam_F - k_lam_G */
                            gsl_sf_result * F, gsl_sf_result * Fp,
                            gsl_sf_result * G, gsl_sf_result * Gp,
                            double * exp_F, double * exp_G)
{
  const double lam_G = lam_F - k_lam_G;

  if(x < 0.0 || lam_F <= -0.5 || lam_G <= -0.5) {
    GSL_SF_RESULT_SET(F,  0.0, 0.0);
    GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
    GSL_SF_RESULT_SET(G,  0.0, 0.0);
    GSL_SF_RESULT_SET(Gp, 0.0, 0.0);
    *exp_F = 0.0;
    *exp_G = 0.0;
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(x == 0.0) {
    gsl_sf_result C0;
    CLeta(0.0, eta, &C0);
    GSL_SF_RESULT_SET(F,  0.0, 0.0);
    GSL_SF_RESULT_SET(Fp, 0.0, 0.0);
    GSL_SF_RESULT_SET(G,  0.0, 0.0); /* FIXME: should be Inf */
    GSL_SF_RESULT_SET(Gp, 0.0, 0.0); /* FIXME: should be Inf */
    *exp_F = 0.0;
    *exp_G = 0.0;
    if(lam_F == 0.0){
      GSL_SF_RESULT_SET(Fp, C0.val, C0.err);
    }
    if(lam_G == 0.0) {
      GSL_SF_RESULT_SET(Gp, 1.0/C0.val, fabs(C0.err/C0.val)/fabs(C0.val));
    }
    GSL_ERROR ("domain error", GSL_EDOM);
    /* After all, since we are asking for G, this is a domain error... */
  }
  else if(x < 1.2 && 2.0*M_PI*eta < 0.9*(-GSL_LOG_DBL_MIN) && fabs(eta*x) < 10.0) {
    /* Reduce to a small lambda value and use the series
     * representations for F and G. We cannot allow eta to
     * be large and positive because the connection formula
     * for G_lam is badly behaved due to an underflow in sin(phi_lam)
     * [see coulomb_FG_series() and coulomb_connection() above].
     * Note that large negative eta is ok however.
     */
    const double SMALL = GSL_SQRT_DBL_EPSILON;
    const int N    = (int)(lam_F + 0.5);
    const int span = GSL_MAX(k_lam_G, N);
    const double lam_min = lam_F - N;    /* -1/2 <= lam_min < 1/2 */
    double F_lam_F, Fp_lam_F;
    double G_lam_G, Gp_lam_G;
    double F_lam_F_err, Fp_lam_F_err;
    double Fp_over_F_lam_F;
    double F_sign_lam_F;
    double F_lam_min_unnorm, Fp_lam_min_unnorm;
    double Fp_over_F_lam_min;
    gsl_sf_result F_lam_min;
    gsl_sf_result G_lam_min, Gp_lam_min;
    double F_scale;
    double Gerr_frac;
    double F_scale_frac_err;
    double F_unnorm_frac_err;

    /* Determine F'/F at lam_F. */
    int CF1_count;
    int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);

    int stat_ser;
    int stat_Fr;
    int stat_Gr;

    /* Recurse down with unnormalized F,F' values. */
    F_lam_F  = SMALL;
    Fp_lam_F = Fp_over_F_lam_F * F_lam_F;
    if(span != 0) {
      stat_Fr = coulomb_F_recur(lam_min, span, eta, x,
                                F_lam_F, Fp_lam_F,
                                &F_lam_min_unnorm, &Fp_lam_min_unnorm
                                );
    }
    else {
      F_lam_min_unnorm  =  F_lam_F;
      Fp_lam_min_unnorm = Fp_lam_F;
      stat_Fr = GSL_SUCCESS;
    }

    /* Determine F and G at lam_min. */
    if(lam_min == -0.5) {
      stat_ser = coulomb_FGmhalf_series(eta, x, &F_lam_min, &G_lam_min);
    }
    else if(lam_min == 0.0) {
      stat_ser = coulomb_FG0_series(eta, x, &F_lam_min, &G_lam_min);
    }
    else if(lam_min == 0.5) {
      /* This cannot happen. */
      F->val  = F_lam_F;
      F->err  = 2.0 * GSL_DBL_EPSILON * fabs(F->val);
      Fp->val = Fp_lam_F;
      Fp->err = 2.0 * GSL_DBL_EPSILON * fabs(Fp->val);
      G->val  = G_lam_G;
      G->err  = 2.0 * GSL_DBL_EPSILON * fabs(G->val);
      Gp->val = Gp_lam_G;
      Gp->err = 2.0 * GSL_DBL_EPSILON * fabs(Gp->val);
      *exp_F = 0.0;
      *exp_G = 0.0;
      GSL_ERROR ("error", GSL_ESANITY);
    }
    else {
      stat_ser = coulomb_FG_series(lam_min, eta, x, &F_lam_min, &G_lam_min);
    }

    /* Determine remaining quantities. */
    Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;
    Gp_lam_min.val  = Fp_over_F_lam_min*G_lam_min.val - 1.0/F_lam_min.val;
    Gp_lam_min.err  = fabs(Fp_over_F_lam_min)*G_lam_min.err;
    Gp_lam_min.err += fabs(1.0/F_lam_min.val) * fabs(F_lam_min.err/F_lam_min.val);
    F_scale     = F_lam_min.val / F_lam_min_unnorm;

    /* Apply scale to the original F,F' values. */
    F_scale_frac_err  = fabs(F_lam_min.err/F_lam_min.val);
    F_unnorm_frac_err = 2.0*GSL_DBL_EPSILON*(CF1_count+span+1);
    F_lam_F     *= F_scale;
    F_lam_F_err  = fabs(F_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);
    Fp_lam_F    *= F_scale;
    Fp_lam_F_err = fabs(Fp_lam_F) * (F_unnorm_frac_err + F_scale_frac_err);

    /* Recurse up to get the required G,G' values. */
    stat_Gr = coulomb_G_recur(lam_min, GSL_MAX(N-k_lam_G,0), eta, x,
                              G_lam_min.val, Gp_lam_min.val,
                              &G_lam_G, &Gp_lam_G
                              );

    F->val  = F_lam_F;
    F->err  = F_lam_F_err;
    F->err += 2.0 * GSL_DBL_EPSILON * fabs(F_lam_F);

    Fp->val  = Fp_lam_F;
    Fp->err  = Fp_lam_F_err;
    Fp->err += 2.0 * GSL_DBL_EPSILON * fabs(Fp_lam_F);

    Gerr_frac = fabs(G_lam_min.err/G_lam_min.val) + fabs(Gp_lam_min.err/Gp_lam_min.val);

    G->val  = G_lam_G;
    G->err  = Gerr_frac * fabs(G_lam_G);
    G->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(G->val);

    Gp->val  = Gp_lam_G;
    Gp->err  = Gerr_frac * fabs(Gp->val);
    Gp->err += 2.0 * (CF1_count+1) * GSL_DBL_EPSILON * fabs(Gp->val);

    *exp_F = 0.0;
    *exp_G = 0.0;

    return GSL_ERROR_SELECT_4(stat_ser, stat_CF1, stat_Fr, stat_Gr);
  }
  else if(x < 2.0*eta) {
    /* Use WKB approximation to obtain F and G at the two
     * lambda values, and use the Wronskian and the
     * continued fractions for F'/F to obtain F' and G'.
     */
    gsl_sf_result F_lam_F, G_lam_F;
    gsl_sf_result F_lam_G, G_lam_G;
    double exp_lam_F, exp_lam_G;
    int stat_lam_F;
    int stat_lam_G;
    int stat_CF1_lam_F;
    int stat_CF1_lam_G;
    int CF1_count;
    double Fp_over_F_lam_F;
    double Fp_over_F_lam_G;
    double F_sign_lam_F;
    double F_sign_lam_G;

    stat_lam_F = coulomb_jwkb(lam_F, eta, x, &F_lam_F, &G_lam_F, &exp_lam_F);
    if(k_lam_G == 0) {
      stat_lam_G = stat_lam_F;
      F_lam_G = F_lam_F;
      G_lam_G = G_lam_F;
      exp_lam_G = exp_lam_F;
    }
    else {
      stat_lam_G = coulomb_jwkb(lam_G, eta, x, &F_lam_G, &G_lam_G, &exp_lam_G);
    }

    stat_CF1_lam_F = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
    if(k_lam_G == 0) {
      stat_CF1_lam_G  = stat_CF1_lam_F;
      F_sign_lam_G    = F_sign_lam_F;
      Fp_over_F_lam_G = Fp_over_F_lam_F;
    }
    else {
      stat_CF1_lam_G = coulomb_CF1(lam_G, eta, x, &F_sign_lam_G, &Fp_over_F_lam_G, &CF1_count);
    }

    F->val = F_lam_F.val;
    F->err = F_lam_F.err;

    G->val = G_lam_G.val;
    G->err = G_lam_G.err;

    Fp->val  = Fp_over_F_lam_F * F_lam_F.val;
    Fp->err  = fabs(Fp_over_F_lam_F) * F_lam_F.err;
    Fp->err += 2.0*GSL_DBL_EPSILON*fabs(Fp->val);

    Gp->val  = Fp_over_F_lam_G * G_lam_G.val - 1.0/F_lam_G.val;
    Gp->err  = fabs(Fp_over_F_lam_G) * G_lam_G.err;
    Gp->err += fabs(1.0/F_lam_G.val) * fabs(F_lam_G.err/F_lam_G.val);

    *exp_F = exp_lam_F;
    *exp_G = exp_lam_G;

    if(stat_lam_F == GSL_EOVRFLW || stat_lam_G == GSL_EOVRFLW) {
      GSL_ERROR ("overflow", GSL_EOVRFLW);
    }
    else {
      return GSL_ERROR_SELECT_2(stat_lam_F, stat_lam_G);
    }
  }
  else {
    /* x > 2 eta, so we know that we can find a lambda value such
     * that x is above the turning point. We do this, evaluate
     * using Steed's method at that oscillatory point, then
     * use recursion on F and G to obtain the required values.
     *
     * lam_0   = a value of lambda such that x is below the turning point
     * lam_min = minimum of lam_0 and the requested lam_G, since
     *           we must go at least as low as lam_G
     */
    const double SMALL = GSL_SQRT_DBL_EPSILON;
    const double C = sqrt(1.0 + 4.0*x*(x-2.0*eta));
    const int N = ceil(lam_F - C + 0.5);
    const double lam_0   = lam_F - GSL_MAX(N, 0);
    const double lam_min = GSL_MIN(lam_0, lam_G);
    double F_lam_F, Fp_lam_F;
    double G_lam_G, Gp_lam_G;
    double F_lam_min_unnorm, Fp_lam_min_unnorm;
    double F_lam_min, Fp_lam_min;
    double G_lam_min, Gp_lam_min;
    double Fp_over_F_lam_F;
    double Fp_over_F_lam_min;
    double F_sign_lam_F, F_sign_lam_min;
    double P_lam_min, Q_lam_min;
    double alpha;
    double gamma;
    double F_scale;

    int CF1_count;
    int CF2_count;
    int stat_CF1 = coulomb_CF1(lam_F, eta, x, &F_sign_lam_F, &Fp_over_F_lam_F, &CF1_count);
    int stat_CF2;
    int stat_Fr;
    int stat_Gr;

    int F_recur_count;
    int G_recur_count;

    double err_amplify;

    F_lam_F  = F_sign_lam_F * SMALL;  /* unnormalized */
    Fp_lam_F = Fp_over_F_lam_F * F_lam_F;

    /* Backward recurrence to get F,Fp at lam_min */
    F_recur_count = GSL_MAX(k_lam_G, N);
    stat_Fr = coulomb_F_recur(lam_min, F_recur_count, eta, x,
                              F_lam_F, Fp_lam_F,
                              &F_lam_min_unnorm, &Fp_lam_min_unnorm
                              );
    Fp_over_F_lam_min = Fp_lam_min_unnorm / F_lam_min_unnorm;

    /* Steed evaluation to complete evaluation of F,Fp,G,Gp at lam_min */
    stat_CF2 = coulomb_CF2(lam_min, eta, x, &P_lam_min, &Q_lam_min, &CF2_count);
    alpha = Fp_over_F_lam_min - P_lam_min;
    gamma = alpha/Q_lam_min;

    F_sign_lam_min = GSL_SIGN(F_lam_min_unnorm) ;

    F_lam_min  = F_sign_lam_min / sqrt(alpha*alpha/Q_lam_min + Q_lam_min);
    Fp_lam_min = Fp_over_F_lam_min * F_lam_min;
    G_lam_min  = gamma * F_lam_min;
    Gp_lam_min = (P_lam_min * gamma - Q_lam_min) * F_lam_min;

    /* Apply scale to values of F,Fp at lam_F (the top). */
    F_scale = F_lam_min / F_lam_min_unnorm;
    F_lam_F  *= F_scale;
    Fp_lam_F *= F_scale;

    /* Forward recurrence to get G,Gp at lam_G (the top). */
    G_recur_count = GSL_MAX(N-k_lam_G,0);
    stat_Gr = coulomb_G_recur(lam_min, G_recur_count, eta, x,
                              G_lam_min, Gp_lam_min,
                              &G_lam_G, &Gp_lam_G
                              );

    err_amplify = CF1_count + CF2_count + F_recur_count + G_recur_count + 1;

    F->val  = F_lam_F;
    F->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(F->val);

    Fp->val = Fp_lam_F;
    Fp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Fp->val);

    G->val  = G_lam_G;
    G->err  = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(G->val);

    Gp->val = Gp_lam_G;
    Gp->err = 8.0*err_amplify*GSL_DBL_EPSILON * fabs(Gp->val);

    *exp_F = 0.0;
    *exp_G = 0.0;

    return GSL_ERROR_SELECT_4(stat_CF1, stat_CF2, stat_Fr, stat_Gr);
  }
}

int
gsl_sf_coulomb_wave_F_array(double lam_min, int kmax,
                                 double eta, double x,
                                 double * fc_array,
                                 double * F_exp)
{
  if(x == 0.0) {
    int k;
    *F_exp = 0.0;
    for(k=0; k<=kmax; k++) {
      fc_array[k] = 0.0;
    }
    if(lam_min == 0.0){
      gsl_sf_result f_0;
      CLeta(0.0, eta, &f_0);
      fc_array[0] = f_0.val;
    }
    return GSL_SUCCESS;
  }
  else {
    const double x_inv = 1.0/x;
    const double lam_max = lam_min + kmax;
    gsl_sf_result F, Fp;
    gsl_sf_result G, Gp;
    double G_exp;

    int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, 0,
                                              &F, &Fp, &G, &Gp, F_exp, &G_exp);

    double fcl  = F.val;
    double fpl = Fp.val;
    double lam = lam_max;
    int k;

    fc_array[kmax] = F.val;

    for(k=kmax-1; k>=0; k--) {
      double el = eta/lam;
      double rl = hypot(1.0, el);
      double sl = el  + lam*x_inv;
      double fc_lm1 = (fcl*sl + fpl)/rl;
      fc_array[k]   = fc_lm1;
      fpl           =  fc_lm1*sl - fcl*rl;
      fcl           =  fc_lm1;
      lam -= 1.0;
    }

    return stat_FG;
  }
}

int
gsl_sf_coulomb_wave_FG_array(double lam_min, int kmax,
                                  double eta, double x,
                                  double * fc_array, double * gc_array,
                                  double * F_exp, double * G_exp)
{
  const double x_inv = 1.0/x;
  const double lam_max = lam_min + kmax;
  gsl_sf_result F, Fp;
  gsl_sf_result G, Gp;

  int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
                                            &F, &Fp, &G, &Gp, F_exp, G_exp);

  double fcl  = F.val;
  double fpl = Fp.val;
  double lam = lam_max;
  int k;

  double gcl, gpl;

  fc_array[kmax] = F.val;

  for(k=kmax-1; k>=0; k--) {
    double el = eta/lam;
    double rl = hypot(1.0, el);
    double sl = el  + lam*x_inv;
    double fc_lm1;
    fc_lm1 = (fcl*sl + fpl)/rl;
    fc_array[k] = fc_lm1;
    fpl         =  fc_lm1*sl - fcl*rl;
    fcl         =  fc_lm1;
    lam -= 1.0;
  }

  gcl = G.val;
  gpl = Gp.val;
  lam = lam_min + 1.0;

  gc_array[0] = G.val;

  for(k=1; k<=kmax; k++) {
    double el = eta/lam;
    double rl = hypot(1.0, el);
    double sl = el + lam*x_inv;
    double gcl1 = (sl*gcl - gpl)/rl;
    gc_array[k] = gcl1;
    gpl         = rl*gcl - sl*gcl1;
    gcl         = gcl1;
    lam += 1.0;
  }

  return stat_FG;
}

int
gsl_sf_coulomb_wave_FGp_array(double lam_min, int kmax,
                                   double eta, double x,
                                   double * fc_array, double * fcp_array,
                                   double * gc_array, double * gcp_array,
                                   double * F_exp, double * G_exp)

{
  const double x_inv = 1.0/x;
  const double lam_max = lam_min + kmax;
  gsl_sf_result F, Fp;
  gsl_sf_result G, Gp;

  int stat_FG = gsl_sf_coulomb_wave_FG_e(eta, x, lam_max, kmax,
                                            &F, &Fp, &G, &Gp, F_exp, G_exp);

  double fcl  = F.val;
  double fpl = Fp.val;
  double lam = lam_max;
  int k;

  double gcl, gpl;

  fc_array[kmax]  = F.val;
  fcp_array[kmax] = Fp.val;

  for(k=kmax-1; k>=0; k--) {
    double el = eta/lam;
    double rl = hypot(1.0, el);
    double sl = el  + lam*x_inv;
    double fc_lm1;
    fc_lm1 = (fcl*sl + fpl)/rl;
    fc_array[k]  = fc_lm1;
    fpl          = fc_lm1*sl - fcl*rl;
    fcp_array[k] = fpl;
    fcl          =  fc_lm1;
    lam -= 1.0;
  }

  gcl = G.val;
  gpl = Gp.val;
  lam = lam_min + 1.0;

  gc_array[0]  = G.val;
  gcp_array[0] = Gp.val;

  for(k=1; k<=kmax; k++) {
    double el = eta/lam;
    double rl = hypot(1.0, el);
    double sl = el + lam*x_inv;
    double gcl1 = (sl*gcl - gpl)/rl;
    gc_array[k]  = gcl1;
    gpl          = rl*gcl - sl*gcl1;
    gcp_array[k] = gpl;
    gcl          = gcl1;
    lam += 1.0;
  }

  return stat_FG;
}

int
gsl_sf_coulomb_wave_sphF_array(double lam_min, int kmax,
                                    double eta, double x,
                                    double * fc_array,
                                    double * F_exp)
{
  if(x < 0.0 || lam_min < -0.5) {
    GSL_ERROR ("domain error", GSL_EDOM);
  }
  else if(x < 10.0/GSL_DBL_MAX) {
    int k;
    for(k=0; k<=kmax; k++) {
      fc_array[k] = 0.0;
    }
    if(lam_min == 0.0) {
      fc_array[0] = sqrt(C0sq(eta));
    }
    *F_exp = 0.0;
    if(x == 0.0)
      return GSL_SUCCESS;
    else
      GSL_ERROR ("underflow", GSL_EUNDRFLW);
  }
  else {
    int k;
    int stat_F = gsl_sf_coulomb_wave_F_array(lam_min, kmax,
                                                  eta, x,
                                                  fc_array,
                                                  F_exp);

    for(k=0; k<=kmax; k++) {
      fc_array[k] = fc_array[k] / x;
    }
    return stat_F;
  }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

/* vector_bfgs2.c -- Fletcher's implementation of the BFGS method,
   from R.Fletcher, "Practical Method's of Optimization", Second
   Edition, ISBN 0471915475.  Algorithms 2.6.2 and 2.6.4. */

/* Thanks to Alan Irwin irwin@beluga.phys.uvic.ca. for suggesting this
   algorithm and providing sample fortran benchmarks */

#include <config.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_blas.h>

#include "linear_minimize.c"
#include "linear_wrapper.c"

typedef struct
{
  int iter;
  double step;
  double g0norm;
  double pnorm;
  double delta_f;
  double fp0;                   /* f'(0) for f(x-alpha*p) */
  gsl_vector *x0;
  gsl_vector *g0;
  gsl_vector *p;
  /* work space */
  gsl_vector *dx0;
  gsl_vector *dg0;
  gsl_vector *x_alpha;
  gsl_vector *g_alpha;
  /* wrapper function */
  wrapper_t wrap;
  /* minimization parameters */
  double rho;
  double sigma;
  double tau1;
  double tau2;
  double tau3;
  int order;
}
vector_bfgs2_state_t;

static int
vector_bfgs2_alloc (void *vstate, size_t n)
{
  vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;

  state->p = gsl_vector_calloc (n);

  if (state->p == 0)
    {
      GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
    }

  state->x0 = gsl_vector_calloc (n);

  if (state->x0 == 0)
    {
      gsl_vector_free (state->p);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->g0 = gsl_vector_calloc (n);

  if (state->g0 == 0)
    {
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->dx0 = gsl_vector_calloc (n);

  if (state->dx0 == 0)
    {
      gsl_vector_free (state->g0);
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->dg0 = gsl_vector_calloc (n);

  if (state->dg0 == 0)
    {
      gsl_vector_free (state->dx0);
      gsl_vector_free (state->g0);
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->x_alpha = gsl_vector_calloc (n);

  if (state->x_alpha == 0)
    {
      gsl_vector_free (state->dg0);
      gsl_vector_free (state->dx0);
      gsl_vector_free (state->g0);
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  state->g_alpha = gsl_vector_calloc (n);

  if (state->g_alpha == 0)
    {
      gsl_vector_free (state->x_alpha);
      gsl_vector_free (state->dg0);
      gsl_vector_free (state->dx0);
      gsl_vector_free (state->g0);
      gsl_vector_free (state->x0);
      gsl_vector_free (state->p);
      GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
    }

  return GSL_SUCCESS;
}

static int
vector_bfgs2_set (void *vstate, gsl_multimin_function_fdf * fdf,
                  const gsl_vector * x, double *f, gsl_vector * gradient,
                  double step_size, double tol)
{
  vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;

  state->iter = 0;
  state->step = step_size;
  state->delta_f = 0;

  GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);

  /* Use the gradient as the initial direction */

  gsl_vector_memcpy (state->x0, x);
  gsl_vector_memcpy (state->g0, gradient);
  state->g0norm = gsl_blas_dnrm2 (state->g0);

  gsl_vector_memcpy (state->p, gradient);
  gsl_blas_dscal (-1 / state->g0norm, state->p);
  state->pnorm = gsl_blas_dnrm2 (state->p);     /* should be 1 */
  state->fp0 = -state->g0norm;

  /* Prepare the wrapper */

  prepare_wrapper (&state->wrap, fdf,
                   state->x0, *f, state->g0,
                   state->p, state->x_alpha, state->g_alpha);

  /* Prepare 1d minimisation parameters */

  state->rho = 0.01;
  state->sigma = tol;
  state->tau1 = 9;
  state->tau2 = 0.05;
  state->tau3 = 0.5;
  state->order = 3;  /* use cubic interpolation where possible */

  return GSL_SUCCESS;
}

static void
vector_bfgs2_free (void *vstate)
{
  vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;

  gsl_vector_free (state->x_alpha);
  gsl_vector_free (state->g_alpha);
  gsl_vector_free (state->dg0);
  gsl_vector_free (state->dx0);
  gsl_vector_free (state->g0);
  gsl_vector_free (state->x0);
  gsl_vector_free (state->p);
}

static int
vector_bfgs2_restart (void *vstate)
{
  vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;

  state->iter = 0;
  return GSL_SUCCESS;
}

static int
vector_bfgs2_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
                      gsl_vector * x, double *f,
                      gsl_vector * gradient, gsl_vector * dx)
{
  vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;
  double alpha = 0.0, alpha1;
  gsl_vector *x0 = state->x0;
  gsl_vector *g0 = state->g0;
  gsl_vector *p = state->p;

  double g0norm = state->g0norm;
  double pnorm = state->pnorm;
  double delta_f = state->delta_f;
  double pg, dir;
  int status;

  double f0 = *f;

  if (pnorm == 0.0 || g0norm == 0.0 || state->fp0 == 0)
    {
      gsl_vector_set_zero (dx);
      return GSL_ENOPROG;
    }

  if (delta_f < 0)
    {
      double del = GSL_MAX_DBL (-delta_f, 10 * GSL_DBL_EPSILON * fabs(f0));
      alpha1 = GSL_MIN_DBL (1.0, 2.0 * del / (-state->fp0));
    }
  else
    {
      alpha1 = fabs(state->step);
    }

  /* line minimisation, with cubic interpolation (order = 3) */

  status = minimize (&state->wrap.fdf_linear, state->rho, state->sigma,
                     state->tau1, state->tau2, state->tau3, state->order,
                     alpha1,  &alpha);

  if (status != GSL_SUCCESS)
    {
      return status;
    }

  update_position (&(state->wrap), alpha, x, f, gradient);

  state->delta_f = *f - f0;

  /* Choose a new direction for the next step */

  {
    /* This is the BFGS update: */
    /* p' = g1 - A dx - B dg */
    /* A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg */
    /* B = dx.g/dx.dg */

    gsl_vector *dx0 = state->dx0;
    gsl_vector *dg0 = state->dg0;

    double dxg, dgg, dxdg, dgnorm, A, B;

    /* dx0 = x - x0 */
    gsl_vector_memcpy (dx0, x);
    gsl_blas_daxpy (-1.0, x0, dx0);

    gsl_vector_memcpy (dx, dx0);  /* keep a copy */

    /* dg0 = g - g0 */
    gsl_vector_memcpy (dg0, gradient);
    gsl_blas_daxpy (-1.0, g0, dg0);

    gsl_blas_ddot (dx0, gradient, &dxg);
    gsl_blas_ddot (dg0, gradient, &dgg);
    gsl_blas_ddot (dx0, dg0, &dxdg);

    dgnorm = gsl_blas_dnrm2 (dg0);

    if (dxdg != 0)
      {
        B = dxg / dxdg;
        A = -(1.0 + dgnorm * dgnorm / dxdg) * B + dgg / dxdg;
      }
    else
      {
        B = 0;
        A = 0;
      }

    gsl_vector_memcpy (p, gradient);
    gsl_blas_daxpy (-A, dx0, p);
    gsl_blas_daxpy (-B, dg0, p);
  }

  gsl_vector_memcpy (g0, gradient);
  gsl_vector_memcpy (x0, x);
  state->g0norm = gsl_blas_dnrm2 (g0);
  state->pnorm = gsl_blas_dnrm2 (p);

  /* update direction and fp0 */

  gsl_blas_ddot (p, gradient, &pg);
  dir = (pg >= 0.0) ? -1.0 : +1.0;
  gsl_blas_dscal (dir / state->pnorm, p);
  state->pnorm = gsl_blas_dnrm2 (p);
  gsl_blas_ddot (p, g0, &state->fp0);

  change_direction (&state->wrap);

  return GSL_SUCCESS;
}

static const gsl_multimin_fdfminimizer_type vector_bfgs2_type = {
  "vector_bfgs2",               /* name */
  sizeof (vector_bfgs2_state_t),
  &vector_bfgs2_alloc,
  &vector_bfgs2_set,
  &vector_bfgs2_iterate,
  &vector_bfgs2_restart,
  &vector_bfgs2_free
};

const gsl_multimin_fdfminimizer_type
  * gsl_multimin_fdfminimizer_vector_bfgs2 = &vector_bfgs2_type;


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef utilities_9904f6c1_3075_4335_b065_c7099f52e80b_h
#define utilities_9904f6c1_3075_4335_b065_c7099f52e80b_h

#include <ariel/type.h>
#include <gslib/std.h>
#include <gslib/string.h>
#include <gslib/file.h>

__ariel_begin__

/*
 * Text format elements:
 *  [section] { ... }
 *  [section]:[notation] { ... }
 *
 * where sections could be nested
 */

extern bool io_bad_eof(const string& src, int32 curr);
extern int32 io_skip_blank_charactors(const string& src, int32 start);
extern int32 io_read_section_name(const string& src, string& name, int32 start);
extern int32 io_read_notation(const string& src, string& notation, int32 start);
extern int32 io_skip_section(const string& src, int32 start);
extern int32 io_enter_section(const string& src, int32 start, gchar st = _t('{'));
extern int32 io_read_line_of_section(const string& src, string& line, int32 start);

/*
 * Binary format elements:
 *  #[uint32][string]@[uint32]      section about length(in bytes)
 *  $[uint32][string]               notation string(ASCII only)
 */

using gs::file;

class __gs_novtable io_binary_stream abstract
{
public:
    typedef vector<int32> section_stack;

    enum control_type
    {
        ctl_unknown,
        ctl_section,
        ctl_notation,
        ctl_counter,
        ctl_byte_stream_field,
        ctl_word_stream_field,
        ctl_dword_stream_field,
        ctl_qword_stream_field,
    };

public:
    io_binary_stream(int32 size);
    virtual ~io_binary_stream() {}
    control_type read_control_type();
    bool next_byte_valid() const { return next_n_bytes_valid(1); }
    bool next_word_valid() const { return next_n_bytes_valid(2); }
    bool next_dword_valid() const { return next_n_bytes_valid(4); }
    bool next_qword_valid() const { return next_n_bytes_valid(8); }
    bool next_n_bytes_valid(int32 bytes) const;
    void seek_to(int32 bytes);
    void seek_by(int32 bytes);
    float read_float();
    double read_double();
    int32 read_nstring(string& str);
    int32 read_string(string& str, const string& stopch);
    void enter_section(int32 size) { _section_stack.push_back(size); }
    bool exit_section();
    bool skip_current_section();
    bool skip_next_section();

protected:
    int32                   _size;
    int32                   _current;
    section_stack           _section_stack;

public:
    virtual byte read_byte() = 0;
    virtual word read_word() = 0;
    virtual dword read_dword() = 0;
    virtual qword read_qword() = 0;
    virtual bool read_field_to_buf(void* ptr, int32 bytes) = 0;
    virtual void rewind_to(int32 bytes) = 0;
    virtual int32 current_dev_pos() const = 0;

protected:
    void rewind_by(int32 bytes);
    void take_next_n_bytes(int32 n);
    bool section_stack_valid(int32 bytes) const;
};

class io_binary_memory:
    public io_binary_stream
{
public:
    io_binary_memory(const void* ptr, int32 size);
    virtual byte read_byte() override;
    virtual word read_word() override;
    virtual dword read_dword() override;
    virtual qword read_qword() override;
    virtual bool read_field_to_buf(void* ptr, int32 bytes) override;
    virtual void rewind_to(int32 bytes) override {}
    virtual int32 current_dev_pos() const override { return _current; }

protected:
    const byte*             _mem;
};

class io_binary_file:
    public io_binary_stream
{
public:
    io_binary_file(file& pf);
    virtual byte read_byte() override;
    virtual word read_word() override;
    virtual dword read_dword() override;
    virtual qword read_qword() override;
    virtual bool read_field_to_buf(void* ptr, int32 bytes) override;
    virtual void rewind_to(int32 bytes) override { _file.seek(bytes, SEEK_SET); }
    virtual int32 current_dev_pos() const override { return _file.current(); }

protected:
    file&                   _file;
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/*
 * author: Achim Gaedke
 * created: January 2001
 * file: pygsl/src/histogrammodule.c
 * $Id: histogrammodule.c,v 1.3 2008/10/26 17:03:23 schnizer Exp $
 *
 *
 * May 2005
 *    Pierre Schnizer
 *    maintainance: replaced direct python error calls with pygsl calls
 *    added CAST and GET macros to reduce code duplication
 *
 *    Now the warnings are all handled by a separate function.
 */

#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>

enum hist_error{
     NOHIST = 0,
     NOHIST2D,
     ARGNOHIST,
     ARGNOHIST2D,
     HISTP_NULL
};
static const char * filename = __FILE__;
static PyObject * module = NULL;
#include "histogram_doc.ic"

static int
PyGSL_hist_error_helper(const char * function, int line, int myerrno, enum hist_error errtype)
{
     char *tmp;
     switch (errtype){
     case NOHIST:      tmp = "Object was not a histogramm"; break;
     case NOHIST2D:    tmp = "Object was not a 2D histogramm"; break;
     case ARGNOHIST:   tmp = "Argument was not a histogramm"; break;
     case ARGNOHIST2D: tmp = "Argument was not a 2D histogramm"; break;
     case HISTP_NULL:  tmp = "Pointer to GSL histogramm(2d) object was NULL!"; break;
     default:          tmp = "Unknown case in function hist_error_helper"; myerrno = GSL_ESANITY; break;
     }
     PyGSL_add_traceback(module, filename, function, line);
     pygsl_error(tmp, filename, line, myerrno);
     return myerrno;
}

/*
 * Check if the recieved object is of approbriate type and that the histogram is defined.
 * Invokes PyGSL_hist_error_helper if it fails.
 * returns GSL_SUCCESS on success
 *
 * ob      ... the object to check
 * type    ... the Python type object
 * errcode ... the errorcode used to describe the failure type for
 *             PyGSL_hist_error_helper
 */

#define _PyGSL_HIST_CHECK_INT(ob, type, errcode) \
     ((ob->ob_type == &(type)) ?  GSL_SUCCESS :  \
     PyGSL_hist_error_helper(__FUNCTION__, __LINE__, (errcode), GSL_ESANITY))

/*
 * returns the gsl histogram from the object
 */
#define _PyGSL_HIST_CAST(ob, type) (( ((type *)(ob)) ->h ))

/*
 * Check if it the gsl_histogram is not NULL. Will give a descriptive error
 * message if it fails.
 */
#define _PyGSL_HIST_CAST_SAVE(ob, type) \
      (\
          (_PyGSL_HIST_CAST((ob), type)  == NULL ) \
        ? \
          (PyGSL_hist_error_helper(__FUNCTION__, __LINE__, GSL_EFAULT, HISTP_NULL), NULL) \
        : \
          ( _PyGSL_HIST_CAST((ob), type) ) \
      )

/*
 * Try to get the gsl_histogram from the python object.
 * Checks for errors ...
 */
#define _PyGSL_HIST_GET_INT(ob, type, cast, errcode) \
      ( \
          (  _PyGSL_HIST_CHECK_INT((ob), type, (errcode)) == GSL_SUCCESS ) \
        ? \
          _PyGSL_HIST_CAST_SAVE((ob), cast) \
	: \
          NULL \
      )
#define _PyGSL_HIST_GET(ob, type, errcode) _PyGSL_HIST_GET_INT(ob, type ## Type, type ## Object, errcode)

#define _PyGSL_HIST_CHECK(ob, errcode)   _PyGSL_HIST_CHECK_INT((ob), (histogram_histogramType), (errcode))
#define _PyGSL_HIST2D_CHECK(ob, errcode) _PyGSL_HIST_CHECK_INT((ob), (histogram_histogram2dType), (errcode))

#define PyGSL_HIST_CHECK(ob)       _PyGSL_HIST_CHECK((ob), (NOHIST))
#define PyGSL_HIST2D_CHECK(ob)     _PyGSL_HIST2D_CHECK((ob), (NOHIST2D))
#define PyGSL_HIST_ARG_CHECK(ob)   _PyGSL_HIST_CHECK((ob), (ARGNOHIST))
#define PyGSL_HIST2D_ARG_CHECK(ob) _PyGSL_HIST2D_CHECK((ob), (ARGNOHIST2D))

#define PyGSL_HIST_CAST(ob)   _PyGSL_HIST_CAST((ob), histogram_histogramObject)
#define PyGSL_HIST2D_CAST(ob) _PyGSL_HIST_CAST((ob), histogram_histogram2dObject)

#define PyGSL_HIST_GET(ob)        _PyGSL_HIST_GET((ob), histogram_histogram,   (NOHIST))
#define PyGSL_HIST2D_GET(ob)      _PyGSL_HIST_GET((ob), histogram_histogram2d, (NOHIST2D))
#define _PyGSL_HIST_GET_ARG(ob)   _PyGSL_HIST_GET((ob), histogram_histogram,   (ARGNOHIST))
#define _PyGSL_HIST2D_GET_ARG(ob) _PyGSL_HIST_GET((ob), histogram_histogram2d, (ARGNOHIST2D))

#define PyGSL_HIST_GET_ARG(ob)   ( ((ob) == NULL) ? NULL : _PyGSL_HIST_GET_ARG((ob)) )
#define PyGSL_HIST2D_GET_ARG(ob) ( ((ob) == NULL) ? NULL : _PyGSL_HIST2D_GET_ARG((ob)) )

/*
 * Helper function for dealing with warning and errors ...
 */
static int
PyGSL_warn_err(int rcode, int errcode, const char * errdes, const char * file, int line)
{
     int warn_result;
     if (errcode==rcode) {
	  warn_result = PyGSL_warning(errdes, file, line, errcode);
	  if (warn_result==-1)
	       /* exception is raised by PyErr_Warn */
	       return GSL_EFAILED;
     }
     else if (PyGSL_ERROR_FLAG(rcode) != GSL_SUCCESS)
	  return rcode;

     return GSL_SUCCESS;
}

#define PyGSL_WARN_ERR(ob, errcode, errdes) \
   (\
       ((ob) == (GSL_SUCCESS)) \
      ? \
         GSL_SUCCESS \
      : \
          ((ob) == (errcode)) \
        ? PyGSL_warn_err(ob, errcode, errdes, filename, __LINE__) \
        : PyGSL_error_flag((ob)) \
    )

static const char edom_message[] =  "value out of histogram range";
#define PyGSL_HIST_EDOM_WARN(ob) PyGSL_WARN_ERR((ob), GSL_EDOM, edom_message)

typedef int (*hist_op)(void *, const void *);
typedef int (*hist_file)(FILE *, void *);

static PyObject *
histogram_histogram_op(PyObject *self, PyObject * arg, hist_op fptr);
static PyObject *
histogram_histogram2d_op(PyObject *self, PyObject * arg, hist_op fptr);
static PyObject *
histogram_histogram_file(PyObject *self, PyObject * arg, hist_file fptr);
static PyObject *
histogram_histogram2d_file(PyObject *self, PyObject * arg, hist_file fptr);

/*
 *
 * histogram type
 * for 1d histogram data
 *
 */

/* my typedef */

staticforward PyTypeObject histogram_histogramType;
staticforward PyMethodDef histogram_histogram_methods[];

typedef struct {
    PyObject_HEAD
    gsl_histogram* h;
} histogram_histogramObject;

#define _CONCAT2(class, suffix) class ## _ ## suffix

#include "histogram.ic"

#define HISTTYPE gsl_histogram
#define PyGSLHISTTYPE histogram_histogramType
#define PyGSL_HIST_TYPE_GET(ob) PyGSL_HIST_GET((ob))
#define PyGSL_HIST_TYPE_CAST(ob) PyGSL_HIST_CAST((ob))
#define PyGSL_HIST_TYPE_ARG_GET(ob) PyGSL_HIST_GET_ARG((ob))
#define GSLNAME(suffix)  _CONCAT2(gsl_histogram, suffix)
#define FUNCNAME(suffix) _CONCAT2(histogram_histogram, suffix)
#ifdef HISTOGRAM2D
#undef HISTOGRAM2D
#endif
#include "histogram_pdf_common.ic"
#include "histogram_common.ic"

#undef HISTTYPE
#undef PyGSLHISTTYPE
#undef PyGSL_HIST_TYPE_GET
#undef PyGSL_HIST_TYPE_CAST
#undef PyGSL_HIST_TYPE_ARG_GET
#undef FUNCNAME
#undef GSLNAME

static
PyTypeObject histogram_histogramType = {
	PyObject_HEAD_INIT(NULL)
	0,
	"pygsl.histogram.histogram",
	sizeof(histogram_histogramObject),
	0,
	(destructor)histogram_histogram_dealloc, /* tp_dealloc */
	0,                                      /* tp_print */
	histogram_histogram_getattr,            /* tp_getattr */
	0,					/* tp_setattr */
	0,			                /* tp_compare */
	0,                  			/* tp_repr */
	0,					/* tp_as_number */
	0,	                		/* tp_as_sequence */
	&histogram_histogram_as_mapping,	/* tp_as_mapping */
	0,				        /* tp_hash */
	0,					/* tp_call */
	0,					/* tp_str */
	0,		                        /* tp_getattro */
	0,					/* tp_setattro */
	0,					/* tp_as_buffer */
	Py_TPFLAGS_DEFAULT,		        /* tp_flags */
        0,				        /* tp_doc */
	0, 		                        /* tp_traverse */
	0,			                /* tp_clear */
	0,              			/* tp_richcompare */
	0,					/* tp_weaklistoffset */
	0,			                /* tp_iter */
	0,					/* tp_iternext */
	0,				        /* tp_methods */
	0,					/* tp_members */
	0,					/* tp_getset */
	0,					/* tp_base */
	0,					/* tp_dict */
	0,					/* tp_descr_get */
	0,					/* tp_descr_set */
	0,					/* tp_dictoffset */
	(initproc)histogram_histogram_init,	/* tp_init */
	NULL,              			/* tp_alloc */
	NULL,                			/* tp_new */
	NULL         			        /* tp_free */
};

/*
 *
 * here begins the section for the 2d histogram
 *
 */

/* my typedef */

staticforward PyTypeObject histogram_histogram2dType;
staticforward PyMethodDef histogram_histogram2d_methods[];
typedef struct {
    PyObject_HEAD
    gsl_histogram2d* h;
} histogram_histogram2dObject;

static PyObject *
histogram_histogram2d_reset(PyObject *);

#include "histogram2d.ic"
#define HISTOGRAM2D 1
#define HISTTYPE gsl_histogram2d
#define PyGSLHISTTYPE histogram_histogram2dType
#define PyGSL_HIST_TYPE_GET(ob) PyGSL_HIST2D_GET((ob))
#define PyGSL_HIST_TYPE_ARG_GET(ob) PyGSL_HIST2D_GET_ARG((ob))
#define PyGSL_HIST_TYPE_CAST(ob) PyGSL_HIST2D_CAST((ob))
#define GSLNAME(suffix)  _CONCAT2(gsl_histogram2d, suffix)
#define FUNCNAME(suffix) _CONCAT2(histogram_histogram2d, suffix)
#include "histogram_pdf_common.ic"
#include "histogram_common.ic"

static
PyTypeObject histogram_histogram2dType = {
	PyObject_HEAD_INIT(NULL)
	0,
	"pygsl.histogram.histogram2d",
	sizeof(histogram_histogram2dObject),
	0,
	(destructor)histogram_histogram2d_dealloc, /* tp_dealloc */
	0,                                      /* tp_print */
	histogram_histogram2d_getattr,          /* tp_getattr */
	0,					/* tp_setattr */
	0,			                /* tp_compare */
	0,                  			/* tp_repr */
	0,					/* tp_as_number */
	0,	                		/* tp_as_sequence */
	&histogram_histogram2d_as_mapping,      /* tp_as_mapping */
	0,				        /* tp_hash */
	0,					/* tp_call */
	0,					/* tp_str */
	0,		                        /* tp_getattro */
	0,					/* tp_setattro */
	0,					/* tp_as_buffer */
	Py_TPFLAGS_DEFAULT,		        /* tp_flags */
        0,				        /* tp_doc */
	0, 		                        /* tp_traverse */
	0,			                /* tp_clear */
	0,              			/* tp_richcompare */
	0,					/* tp_weaklistoffset */
	0,			                /* tp_iter */
	0,					/* tp_iternext */
	0,				        /* tp_methods */
	0,					/* tp_members */
	0,					/* tp_getset */
	0,					/* tp_base */
	0,					/* tp_dict */
	0,					/* tp_descr_get */
	0,					/* tp_descr_set */
	0,					/* tp_dictoffset */
	(initproc)histogram_histogram2d_init,	/* tp_init */
	NULL,            			/* tp_alloc */
	NULL,			                /* tp_new */
	NULL			                /* tp_free */
};

#include "histogram_pdf.ic"
/*
 *
 * module specific stuff
 *
 */

static PyMethodDef histogramMethods[] = {
  {NULL, NULL, 0, NULL}        /* Sentinel */
};

void
register_type(PyTypeObject *p, char *name)
{
     p->ob_type  = &PyType_Type;
     p->tp_alloc = PyType_GenericAlloc;
     p->tp_new   = PyType_GenericNew;
     p->tp_free  = _PyObject_Del;
     /* install histogram type */
     /* important! must increment histogram type reference counter */
     Py_INCREF((PyObject*)p);
     PyModule_AddObject(module, name, (PyObject*)p);
}

void
inithistogram(void)
{
  PyObject* m;
  m=Py_InitModule("histogram", histogramMethods);
  if(!m)
       return;

  module = m;
  init_pygsl();
  /* init histogram type */
  register_type(&histogram_histogramType, "histogram");
  register_type(&histogram_histogram_pdfType, "histogram_pdf");
  register_type(&histogram_histogram2dType, "histogram2d");
  register_type(&histogram_histogram2d_pdfType, "histogram2d_pdf");
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


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

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_linalg.h>

#include "apply_givens.c"

/* Factorise a general M x N matrix A into
 *
 *   A P = Q R
 *
 * where Q is orthogonal (M x M) and R is upper triangular (M x N).
 * When A is rank deficient, r = rank(A) < n, then the permutation is
 * used to ensure that the lower n - r rows of R are zero and the first
 * r columns of Q form an orthonormal basis for A.
 *
 * Q is stored as a packed set of Householder transformations in the
 * strict lower triangular part of the input matrix.
 *
 * R is stored in the diagonal and upper triangle of the input matrix.
 *
 * P: column j of P is column k of the identity matrix, where k =
 * permutation->data[j]
 *
 * The full matrix for Q can be obtained as the product
 *
 *       Q = Q_k .. Q_2 Q_1
 *
 * where k = MIN(M,N) and
 *
 *       Q_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [1, m(i+1,i), m(i+2,i), ... , m(M,i)]
 *
 * This storage scheme is the same as in LAPACK.  See LAPACK's
 * dgeqpf.f for details.
 *
 */

int
gsl_linalg_QRPT_decomp (gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (p->size != N)
    {
      GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
    }
  else if (norm->size != N)
    {
      GSL_ERROR ("norm size must be N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      *signum = 1;

      gsl_permutation_init (p); /* set to identity */

      /* Compute column norms and store in workspace */

      for (i = 0; i < N; i++)
        {
          gsl_vector_view c = gsl_matrix_column (A, i);
          double x = gsl_blas_dnrm2 (&c.vector);
          gsl_vector_set (norm, i, x);
        }

      for (i = 0; i < GSL_MIN (M, N); i++)
        {
          /* Bring the column of largest norm into the pivot position */

          double max_norm = gsl_vector_get(norm, i);
          size_t j, kmax = i;

          for (j = i + 1; j < N; j++)
            {
              double x = gsl_vector_get (norm, j);

              if (x > max_norm)
                {
                  max_norm = x;
                  kmax = j;
                }
            }

          if (kmax != i)
            {
              gsl_matrix_swap_columns (A, i, kmax);
              gsl_permutation_swap (p, i, kmax);
              gsl_vector_swap_elements(norm,i,kmax);

              (*signum) = -(*signum);
            }

          /* Compute the Householder transformation to reduce the j-th
             column of the matrix to a multiple of the j-th unit vector */

          {
            gsl_vector_view c_full = gsl_matrix_column (A, i);
            gsl_vector_view c = gsl_vector_subvector (&c_full.vector,
                                                      i, M - i);
            double tau_i = gsl_linalg_householder_transform (&c.vector);

            gsl_vector_set (tau, i, tau_i);

            /* Apply the transformation to the remaining columns */

            if (i + 1 < N)
              {
                gsl_matrix_view m = gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i+1));

                gsl_linalg_householder_hm (tau_i, &c.vector, &m.matrix);
              }
          }

          /* Update the norms of the remaining columns too */

          if (i + 1 < M)
            {
              for (j = i + 1; j < N; j++)
                {
                  double x = gsl_vector_get (norm, j);

                  if (x > 0.0)
                    {
                      double y = 0;
                      double temp= gsl_matrix_get (A, i, j) / x;

                      if (fabs (temp) >= 1)
                        y = 0.0;
                      else
                        y = x * sqrt (1 - temp * temp);

                      /* recompute norm to prevent loss of accuracy */

                      if (fabs (y / x) < sqrt (20.0) * GSL_SQRT_DBL_EPSILON)
                        {
                          gsl_vector_view c_full = gsl_matrix_column (A, j);
                          gsl_vector_view c =
                            gsl_vector_subvector(&c_full.vector,
                                                 i+1, M - (i+1));
                          y = gsl_blas_dnrm2 (&c.vector);
                        }

                      gsl_vector_set (norm, j, y);
                    }
                }
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QRPT_decomp2 (const gsl_matrix * A, gsl_matrix * q, gsl_matrix * r, gsl_vector * tau, gsl_permutation * p, int *signum, gsl_vector * norm)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (q->size1 != M || q->size2 !=M)
    {
      GSL_ERROR ("q must be M x M", GSL_EBADLEN);
    }
  else if (r->size1 != M || r->size2 !=N)
    {
      GSL_ERROR ("r must be M x N", GSL_EBADLEN);
    }
  else if (tau->size != GSL_MIN (M, N))
    {
      GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
    }
  else if (p->size != N)
    {
      GSL_ERROR ("permutation size must be N", GSL_EBADLEN);
    }
  else if (norm->size != N)
    {
      GSL_ERROR ("norm size must be N", GSL_EBADLEN);
    }

  gsl_matrix_memcpy (r, A);

  gsl_linalg_QRPT_decomp (r, tau, p, signum, norm);

  /* FIXME:  aliased arguments depends on behavior of unpack routine! */

  gsl_linalg_QR_unpack (r, tau, q, r);

  return GSL_SUCCESS;
}

/* Solves the system A x = b using the Q R P^T factorisation,

   R z = Q^T b

   x = P z;

   to obtain x. Based on SLATEC code. */

int
gsl_linalg_QRPT_solve (const gsl_matrix * QR,
                       const gsl_vector * tau,
                       const gsl_permutation * p,
                       const gsl_vector * b,
                       gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != p->size)
    {
      GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
    }
  else if (QR->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_memcpy (x, b);

      gsl_linalg_QRPT_svx (QR, tau, p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QRPT_svx (const gsl_matrix * QR,
                     const gsl_vector * tau,
                     const gsl_permutation * p,
                     gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != p->size)
    {
      GSL_ERROR ("matrix size must match permutation size", GSL_EBADLEN);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* compute sol = Q^T b */

      gsl_linalg_QR_QTvec (QR, tau, x);

      /* Solve R x = sol, storing x inplace in sol */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

/* Find the least squares solution to the overdetermined system
 *
 *   A x = b
 *
 * for M >= N using the QRPT factorization A P = Q R. Assumes
 * A has full column rank.
 */

int
gsl_linalg_QRPT_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p,
                         const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
{
  const size_t N = QR->size2;
  int status = gsl_linalg_QRPT_lssolve2(QR, tau, p, b, N, x, residual);
  return status;
}

/* Find the least squares solution to the overdetermined system
 *
 *   A x = b
 *
 * for M >= N using the QRPT factorization A P = Q R, where A
 * is assumed rank deficient with a given rank.
 */

int
gsl_linalg_QRPT_lssolve2 (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p,
                          const gsl_vector * b, const size_t rank, gsl_vector * x, gsl_vector * residual)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("QR matrix must have M>=N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (rank == 0 || rank > N)
    {
      GSL_ERROR ("rank must have 0 < rank <= N", GSL_EBADLEN);
    }
  else if (N != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (M != residual->size)
    {
      GSL_ERROR ("matrix size must match residual size", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view R11 = gsl_matrix_const_submatrix (QR, 0, 0, rank, rank);
      gsl_vector_view QTb1 = gsl_vector_subvector(residual, 0, rank);
      gsl_vector_view x1 = gsl_vector_subvector(x, 0, rank);
      size_t i;

      /* compute work = Q^T b */
      gsl_vector_memcpy(residual, b);
      gsl_linalg_QR_QTvec (QR, tau, residual);

      /* solve R_{11} x(1:r) = [Q^T b](1:r) */
      gsl_vector_memcpy(&(x1.vector), &(QTb1.vector));
      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, &(R11.matrix), &(x1.vector));

      /* x(r+1:N) = 0 */
      for (i = rank; i < N; ++i)
        gsl_vector_set(x, i, 0.0);

      /* compute x = P y */
      gsl_permute_vector_inverse (p, x);

      /* compute residual = b - A x = Q (Q^T b - R x) */
      gsl_vector_set_zero(&(QTb1.vector));
      gsl_linalg_QR_Qvec(QR, tau, residual);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q, const gsl_matrix * R,
                         const gsl_permutation * p,
                         const gsl_vector * b,
                         gsl_vector * x)
{
  if (Q->size1 != Q->size2 || R->size1 != R->size2)
    {
      return GSL_ENOTSQR;
    }
  else if (Q->size1 != p->size || Q->size1 != R->size1
           || Q->size1 != b->size)
    {
      return GSL_EBADLEN;
    }
  else
    {
      /* compute b' = Q^T b */

      gsl_blas_dgemv (CblasTrans, 1.0, Q, b, 0.0, x);

      /* Solve R x = b', storing x inplace */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, R, x);

      /* Apply permutation to solution in place */

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR,
                        const gsl_permutation * p,
                        const gsl_vector * b,
                        gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
    }
  else if (p->size != x->size)
    {
      GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
    }
  else
    {
      /* Copy x <- b */

      gsl_vector_memcpy (x, b);

      /* Solve R x = b, storing x inplace */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR,
                      const gsl_permutation * p,
                      gsl_vector * x)
{
  if (QR->size1 != QR->size2)
    {
      GSL_ERROR ("QR matrix must be square", GSL_ENOTSQR);
    }
  else if (QR->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match x size", GSL_EBADLEN);
    }
  else if (p->size != x->size)
    {
      GSL_ERROR ("permutation size must match x size", GSL_EBADLEN);
    }
  else
    {
      /* Solve R x = b, storing x inplace */

      gsl_blas_dtrsv (CblasUpper, CblasNoTrans, CblasNonUnit, QR, x);

      gsl_permute_vector_inverse (p, x);

      return GSL_SUCCESS;
    }
}

/* Update a Q R P^T factorisation for A P= Q R ,  A' = A + u v^T,

   Q' R' P^-1 = QR P^-1 + u v^T
   = Q (R + Q^T u v^T P ) P^-1
   = Q (R + w v^T P) P^-1

   where w = Q^T u.

   Algorithm from Golub and Van Loan, "Matrix Computations", Section
   12.5 (Updating Matrix Factorizations, Rank-One Changes)  */

int
gsl_linalg_QRPT_update (gsl_matrix * Q, gsl_matrix * R,
                        const gsl_permutation * p,
                        gsl_vector * w, const gsl_vector * v)
{
  const size_t M = R->size1;
  const size_t N = R->size2;

  if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q matrix must be M x M if R is M x N", GSL_ENOTSQR);
    }
  else if (w->size != M)
    {
      GSL_ERROR ("w must be length M if R is M x N", GSL_EBADLEN);
    }
  else if (v->size != N)
    {
      GSL_ERROR ("v must be length N if R is M x N", GSL_EBADLEN);
    }
  else
    {
      size_t j, k;
      double w0;

      /* Apply Given's rotations to reduce w to (|w|, 0, 0, ... , 0)

         J_1^T .... J_(n-1)^T w = +/- |w| e_1

         simultaneously applied to R,  H = J_1^T ... J^T_(n-1) R
         so that H is upper Hessenberg.  (12.5.2) */

      for (k = M - 1; k > 0; k--)
        {
          double c, s;
          double wk = gsl_vector_get (w, k);
          double wkm1 = gsl_vector_get (w, k - 1);

          gsl_linalg_givens (wkm1, wk, &c, &s);
          gsl_linalg_givens_gv (w, k - 1, k, c, s);
          apply_givens_qr (M, N, Q, R, k - 1, k, c, s);
        }

      w0 = gsl_vector_get (w, 0);

      /* Add in w v^T  (Equation 12.5.3) */

      for (j = 0; j < N; j++)
        {
          double r0j = gsl_matrix_get (R, 0, j);
          size_t p_j = gsl_permutation_get (p, j);
          double vj = gsl_vector_get (v, p_j);
          gsl_matrix_set (R, 0, j, r0j + w0 * vj);
        }

      /* Apply Givens transformations R' = G_(n-1)^T ... G_1^T H
         Equation 12.5.4 */

     for (k = 1; k < GSL_MIN(M,N+1); k++)
        {
          double c, s;
          double diag = gsl_matrix_get (R, k - 1, k - 1);
          double offdiag = gsl_matrix_get (R, k, k - 1);

          gsl_linalg_givens (diag, offdiag, &c, &s);
          apply_givens_qr (M, N, Q, R, k - 1, k, c, s);

          gsl_matrix_set (R, k, k - 1, 0.0);    /* exact zero of G^T */
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_QRPT_rank()
  Estimate rank of triangular matrix from QRPT decomposition

Inputs: QR  - QRPT decomposed matrix
        tol - tolerance for rank determination; if < 0,
              a default value is used:
              20 * (M + N) * eps(max(|diag(R)|))

Return: rank estimate
*/

size_t
gsl_linalg_QRPT_rank (const gsl_matrix * QR, const double tol)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;
  gsl_vector_const_view diag = gsl_matrix_const_diagonal(QR);
  double eps;
  size_t i;
  size_t r = 0;

  if (tol < 0.0)
    {
      double min, max, absmax;
      int ee;

      gsl_vector_minmax(&diag.vector, &min, &max);
      absmax = GSL_MAX(fabs(min), fabs(max));
      ee = (int) (log(absmax) / M_LN2); /* can't use log2 since its not ANSI */

      eps = 20.0 * (M + N) * pow(2.0, (double) ee) * GSL_DBL_EPSILON;
    }
  else
    eps = tol;

  /* count number of diagonal elements with |di| > eps */
  for (i = 0; i < GSL_MIN(M, N); ++i)
    {
      double di = gsl_vector_get(&diag.vector, i);
      if (fabs(di) > eps)
        ++r;
    }

  return r;
}

int
gsl_linalg_QRPT_rcond(const gsl_matrix * QR, double * rcond, gsl_vector * work)
{
  const size_t M = QR->size1;
  const size_t N = QR->size2;

  if (M < N)
    {
      GSL_ERROR ("M must be >= N", GSL_EBADLEN);
    }
  else if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_const_view R = gsl_matrix_const_submatrix (QR, 0, 0, N, N);
      int status;

      status = gsl_linalg_tri_upper_rcond(&R.matrix, rcond, work);

      return status;
    }
}