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

/*----------------------------------------------------------------------*
 * "The purpose of computing is insight, not numbers." - R.W. Hamming   *
 * Hermite polynomials, Hermite functions                               *
 * and their respective arbitrary derivatives                           *
 *----------------------------------------------------------------------*/

/* TODO:
 * - array functions for derivatives of Hermite functions
 * - asymptotic approximation for derivatives of Hermite functions
 * - refine existing asymptotic approximations, especially around x=sqrt(2*n+1) or x=sqrt(2*n+1)*sqrt(2), respectively
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_hermite.h>

#include "error.h"
#include "eval.h"

#define pow2(n) (gsl_sf_pow_int(2,n))
#define RND(x)  ((double) ((x >= 0) ? (int) (x + 0.5) : (int) (x - 0.5)))

/* evaluates the probabilists' Hermite polynomial of order n at position x */
int
gsl_sf_hermite_prob_e(const int n, const double x, gsl_sf_result * result)
{
  if (n < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (n == 0)
    {
      result->val = 1.;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else if (n == 1)
    {
      result->val = x;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else if (x == 0.0)
    {
      if (GSL_IS_ODD(n))
        {
          result->val = 0.0;
          result->err = 0.0;
          return GSL_SUCCESS;
        }
      else
        {
          /* for n even, He_n(0) = (-1)^{n/2} (n - 1)!! */

          int status = GSL_SUCCESS;

          if (n - 1 > GSL_SF_DOUBLEFACT_NMAX)  /* test if (n - 1)!! will overflow */
            {
              status = GSL_EOVRFLW;
              result->val = GSL_IS_ODD(n / 2) ? GSL_NEGINF : GSL_POSINF;
              result->err = GSL_POSINF;
            }
          else
            {
              gsl_sf_doublefact_e(n - 1, result);
              if (GSL_IS_ODD(n / 2))
                result->val = -result->val;
            }

          return status;
        }
    }
  else
    {
      /* upward recurrence: He_{n+1} = x He_n - n He_{n-1} */

      int status = GSL_SUCCESS;
      const double abs_x = fabs(x);
      const double thresh1 = abs_x > 1.0 ? 0.9 * GSL_DBL_MAX / abs_x : GSL_DBL_MAX;
      const double thresh2 = 0.9 * GSL_DBL_MAX;

      double p_n0 = 1.0;    /* He_0(x) */
      double p_n1 = x;      /* He_1(x) */
      double p_n = p_n1;

      double e_n0 = GSL_DBL_EPSILON;
      double e_n1 = fabs(x)*GSL_DBL_EPSILON;
      double e_n = e_n1;

      int j;

      for (j = 1; j < n; j++)
        {
          if (fabs(p_n1) > thresh1 ||    /* test if x*p_n1 will overflow */
              fabs(p_n0) > thresh2 / j)  /* test if j*p_n0 will overflow */
            {
              status = GSL_EOVRFLW;
              break;
            }

          p_n  = x*p_n1 - j*p_n0;
          p_n0 = p_n1;
          p_n1 = p_n;

          e_n  = fabs(x)*e_n1+j*e_n0;
          e_n0 = e_n1;
          e_n1 = e_n;
        }

      result->val = p_n;
      result->err = e_n + fabs(result->val)*GSL_DBL_EPSILON;

      return status;
    }
}

double gsl_sf_hermite_prob(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_prob_e(n, x, &result));
}

/* Evaluates the m-th derivative of the probabilists' Hermite polynomial of order n at position x.
 * The direct formula He^{(m)}_n = n!/(n-m)!*He_{n-m}(x) (where He_j(x) is the j-th probabilists' Hermite polynomial and He^{(m)}_j(x) its m-th derivative) is employed. */
int
gsl_sf_hermite_prob_deriv_e(const int m, const int n, const double x, gsl_sf_result * result)
{
  if (n < 0 || m < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (n < m)
    {
      result->val = 0.;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else
    {
      int status;
      double f = gsl_sf_choose(n,m)*gsl_sf_fact(m);
      gsl_sf_result He;

      status = gsl_sf_hermite_prob_e(n - m, x, &He);
      if (status == GSL_SUCCESS)
        {
          result->val = He.val*f;
          result->err = He.err*f + GSL_DBL_EPSILON*fabs(result->val);
        }
      else
        {
          result->val = He.val;
          result->err = GSL_POSINF;
        }

      return status;
    }
}

double
gsl_sf_hermite_prob_deriv(const int m, const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_prob_deriv_e(m, n, x, &result));
}

/* evaluates the physicists' Hermite polynomial of order n at position x */
int
gsl_sf_hermite_e(const int n, const double x, gsl_sf_result * result)
{
  if (n < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (n == 0)
    {
      result->val = 1.;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else if (n == 1)
    {
      result->val = 2.0*x;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else if (x == 0.0)
    {
      if (GSL_IS_ODD(n))
        {
          result->val = 0.0;
          result->err = 0.0;
          return GSL_SUCCESS;
        }
      else
        {
          /* for n even, H_n(0) = (-2)^{n/2} (n - 1)!! */

          int status = GSL_SUCCESS;
          int m = n >> 1;

          if (n - 1 > GSL_SF_DOUBLEFACT_NMAX)  /* test if (n - 1)!! will overflow */
            {
              status = GSL_EOVRFLW;
              result->val = GSL_IS_ODD(m) ? GSL_NEGINF : GSL_POSINF;
              result->err = GSL_POSINF;
            }
          else
            {
              double f = gsl_pow_int(2.0, m);

              gsl_sf_doublefact_e(n - 1, result);

              if (result->val > 0.9 * GSL_DBL_MAX / f) /* test if 2^{n/2} * (n-1)!! will overflow */
                {
                  status = GSL_EOVRFLW;
                  result->val = GSL_IS_ODD(m) ? GSL_NEGINF : GSL_POSINF;
                  result->err = GSL_POSINF;
                }
              else
                {
                  result->val *= f;
                  result->err *= f;
                  if (GSL_IS_ODD(m))
                    result->val = -result->val;
                }
            }

          return status;
        }
    }
  else
    {
      /* upward recurrence: H_{n+1} = 2x H_n - 2n H_{n-1} */

      int status = GSL_SUCCESS;
      const double two_x = 2.0 * x;
      const double abs_two_x = fabs(two_x);
      const double thresh1 = abs_two_x > 1.0 ? 0.9 * GSL_DBL_MAX / abs_two_x : GSL_DBL_MAX;
      const double thresh2 = 0.9 * GSL_DBL_MAX / 2.0;

      double p_n0 = 1.0;    /* H_0(x) */
      double p_n1 = two_x;  /* H_1(x) */
      double p_n = p_n1;

      double e_n0 = GSL_DBL_EPSILON;
      double e_n1 = 2.*fabs(x)*GSL_DBL_EPSILON;
      double e_n = e_n1;

      int j;

      for (j = 1; j <= n - 1; j++)
        {
          if (fabs(p_n1) > thresh1 ||    /* test if 2*x*p_n1 will overflow */
              fabs(p_n0) > thresh2 / j)  /* test if 2*j*p_n0 will overflow */
            {
              status = GSL_EOVRFLW;
              break;
            }

          p_n  = two_x*p_n1 - 2.0*j*p_n0;
          p_n0 = p_n1;
          p_n1 = p_n;

          e_n  = 2.*(fabs(x)*e_n1+j*e_n0);
          e_n0 = e_n1;
          e_n1 = e_n;
        }

      result->val = p_n;
      result->err = e_n + fabs(result->val)*GSL_DBL_EPSILON;

      return status;
    }
}

double
gsl_sf_hermite(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_e(n, x, &result));
}

/* Evaluates the m-th derivative of the physicists' Hermite polynomial of order n at position x.
 * The direct formula H^{(m)}_n = 2**m*n!/(n-m)!*H_{n-m}(x) (where H_j(x) is the j-th physicists' Hermite polynomial and H^{(m)}_j(x) its m-th derivative) is employed. */
int
gsl_sf_hermite_deriv_e(const int m, const int n, const double x, gsl_sf_result * result)
{
  if (n < 0 || m < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (n < m)
    {
      result->val = 0.;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else
    {
      int status;
      double f = gsl_sf_choose(n,m)*gsl_sf_fact(m)*pow2(m);
      gsl_sf_result H;

      status = gsl_sf_hermite_e(n - m, x, &H);
      if (status == GSL_SUCCESS)
        {
          result->val = H.val*f;
          result->err = H.err*f + GSL_DBL_EPSILON*fabs(result->val);
        }
      else
        {
          result->val = H.val;
          result->err = GSL_POSINF;
        }

      return status;
    }
}

double
gsl_sf_hermite_deriv(const int m, const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_deriv_e(m, n, x, &result));
}

/* evaluates the Hermite function of order n at position x */
int
gsl_sf_hermite_func_e(const int n, const double x, gsl_sf_result * result)
{
  if (n < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (x == 0.0)
    {
      if (GSL_IS_ODD(n))
        {
          result->val = 0.;
          result->err = 0.;
          return GSL_SUCCESS;
        }
      else
        {
          double f = (GSL_IS_ODD(n / 2) ? -1.0 : 1.0);
          int j;

          for(j = 1; j < n; j += 2)
            f *= sqrt(j / (j + 1.0));

          result->val = f / sqrt(M_SQRTPI);
          result->err = GSL_DBL_EPSILON * fabs(result->val);
          return GSL_SUCCESS;
        }
    }
  else if (n == 0)
    {
      result->val = exp(-0.5 * x * x) / sqrt(M_SQRTPI);
      result->err = GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
  else if (n == 1)
    {
      result->val = M_SQRT2 * x * exp(-0.5 * x * x) / sqrt(M_SQRTPI);
      result->err = GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
  else
    {
      /*
       * This algorithm is based on the modified recurrence algorithm
       * found in the appendix of:
       *
       * B. Bunck, BIT Numerical Mathematics, 49, 281 (2009)
       *
       * Numerical tests showed that this algorithm is more stable for
       * large x (x > 40) than the standard recurrence relation.
       * Accuracy is comparable to the recurrence relation method
       * for small x and all n.
       *
       * See:
       *
       * https://scicomp.stackexchange.com/questions/30896/generate-high-n-quantum-harmonic-oscillator-states-numerically
       *
       * for further discussion.
       */

      double hi2 = 1.0 / sqrt(M_SQRTPI); /* \hat{h}_0 */
      double hi1 = M_SQRT2 * x * hi2;    /* \hat{h}_1 */
      double hi = 0.0;
      double sum_log_scale = 0.0;
      double abshi;
      int i;

      for (i = 2; i <= n; ++i)
        {
          hi = sqrt(2.0 / i) * x * hi1 - sqrt((i - 1.0) / i) * hi2;
          hi2 = hi1;
          hi1 = hi;

          abshi = fabs(hi);
          if (abshi > 1.0)
            {
              double log_scale = RND(log(abshi));
              double scale = exp(-log_scale);

              hi *= scale;
              hi1 *= scale;
              hi2 *= scale;
              sum_log_scale += log_scale;
            }
        }

      result->val = hi * exp(-0.5 * x * x + sum_log_scale);
      result->err = n * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_SUCCESS;
    }
}

double
gsl_sf_hermite_func(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_func_e(n, x, &result));
}

/*
 * This algorithm is based on the contour integral algorithm of:
 *
 * B. Bunck, BIT Numerical Mathematics, 49, 281 (2009)
 *
 * It has O(sqrt(n)) complexity
 */

int
gsl_sf_hermite_func_fast_e(const int n, const double x, gsl_sf_result * result)
{
  if (n < 1000 || x == 0.0)
    {
      /* for small n, the recurrence method is faster and more accurate */
      return gsl_sf_hermite_func_e(n, x, result);
    }
  else
    {
      size_t j;
      const double k = sqrt(0.5*n);
      const size_t steps = (size_t) ceil(6.211 * sqrt(n));
      const double dt = M_PI/steps;
      const double invn2 = 1.0/(n*n);
      double ex, ex_e, cs, cs_e, sn, sn2, t;
      gsl_sf_result lngamma;

      if (n < 36)
        {
          gsl_sf_lnfact_e(n, &lngamma);
          lngamma.val *= 0.5;
          lngamma.err *= 0.5;
          t = 0.5*n*log(n) + 0.25*M_LNPI;
          cs = 0.5*n;
          lngamma.val += cs - t;
          lngamma.err += (cs + t)*GSL_DBL_EPSILON;
        }
      else
        {
          /* approximate ln(gamma_{n,k}) using Stirling's formula */
          lngamma.val = 0.25*log(2*n);
          lngamma.err = (lngamma.val  + ((((invn2/3360 + 1.0/2520)*invn2 + 1.0/720)*invn2) + 1.0/24)/n)*GSL_DBL_EPSILON;
          lngamma.val -= ((((invn2/3360 - 1.0/2520)*invn2 + 1.0/720)*invn2) - 1.0/24)/n;
        }

      ex = exp(lngamma.val - n - 0.5*x*x - 2*x*k);
      cs = (GSL_IS_ODD(n) ? -1 : 1);
      result->val = 0.5*ex*cs;
      result->err = 0.5*ex*(lngamma.err + (n + 0.5*x*x + fabs(2*x*k) + 1)*GSL_DBL_EPSILON);
      ex = exp(lngamma.val - n - 0.5*x*x + 2*x*k);
      result->val += 0.5*ex;
      result->err += 0.5*ex*(lngamma.err + (n + 0.5*x*x + fabs(2*x*k) + 1)*GSL_DBL_EPSILON);

      for (j = 1; j < steps; j++)
        {
          t = j*dt;
          cs = cos(t);
          ex = exp(lngamma.val - 0.5*x*x + (2*x*k - n*cs)*cs);
          ex_e = ex*(lngamma.err + GSL_DBL_EPSILON*(1 + 0.5*x*x + (fabs(2*x*k) + fabs(n*cs))*fabs(cs)));
          sn = sin(t);
          sn2 = sin(2*t);
          cs = cos(2*x*k*sn - 0.5*n*sn2 - n*t);
          cs_e = GSL_MIN(1.0+fabs(cs), GSL_DBL_EPSILON*(fabs(cs) + (fabs(2*x*k*sn) + fabs(0.5*n*sn2) + n*t)*fabs(sin(2*x*k*sn - 0.5*n*sn2 - n*t))));
          result->val += ex*cs;
          result->err += ex*cs_e + ex_e*fabs(cs) + GSL_DBL_EPSILON*fabs(ex*cs);
        }

      result->val *= M_1_PI*dt;
      result->err = M_1_PI*dt*result->err + GSL_DBL_EPSILON*fabs(result->val);

      return GSL_SUCCESS;
    }
}

double
gsl_sf_hermite_func_fast(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_func_fast_e(n, x, &result));
}

/* Evaluates all probabilists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */
int
gsl_sf_hermite_prob_array(const int nmax, const double x, double * result_array)
{
  if (nmax < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (nmax == 0)
    {
      result_array[0] = 1.0;
      return GSL_SUCCESS;
    }
  else if (nmax == 1)
    {
      result_array[0] = 1.0;
      result_array[1] = x;
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence: He_{n+1} = x He_n - n He_{n-1} */

      int status = GSL_SUCCESS;
      const double abs_x = fabs(x);
      const double thresh1 = abs_x > 1.0 ? 0.9 * GSL_DBL_MAX / abs_x : GSL_DBL_MAX;
      const double thresh2 = 0.9 * GSL_DBL_MAX;

      double p_n0 = 1.0;    /* He_0(x) */
      double p_n1 = x;      /* He_1(x) */
      double p_n = p_n1;

      int j;

      result_array[0] = 1.0;
      result_array[1] = x;

      for (j = 1; j < nmax; j++)
        {
          if (fabs(p_n1) > thresh1 ||    /* test if x*p_n1 will overflow */
              fabs(p_n0) > thresh2 / j)  /* test if j*p_n0 will overflow */
            {
              status = GSL_EOVRFLW;
              break;
            }

          p_n  = x*p_n1 - j*p_n0;
          p_n0 = p_n1;
          p_n1 = p_n;

          result_array[j + 1] = p_n;
        }

      return status;
   }
}

/* Evaluates the m-th derivative of all probabilists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */

int
gsl_sf_hermite_prob_array_deriv(const int m, const int nmax, const double x, double * result_array)
{
  if (nmax < 0 || m < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (m == 0)
    {
      gsl_sf_hermite_prob_array(nmax, x, result_array);
      return GSL_SUCCESS;
    }
  else if (nmax < m)
    {
      int j;

      for (j = 0; j <= nmax; j++)
        result_array[j] = 0.0;

      return GSL_SUCCESS;
    }
  else if (nmax == m)
    {
      int j;

      for (j = 0; j < m; j++)
        result_array[j] = 0.0;

      result_array[nmax] = gsl_sf_fact(m);
      return GSL_SUCCESS;
    }
  else if (nmax == m + 1)
    {
      int j;

      for (j = 0; j < m; j++)
        result_array[j] = 0.0;

      result_array[nmax-1] = gsl_sf_fact(m);
      result_array[nmax] = result_array[nmax-1]*(m+1)*x;
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence: He^{(m)}_{n+1} = (n+1)/(n-m+1)*(x He^{(m)}_n - n He^{(m)}_{n-1}) */

      double p_n0 = gsl_sf_fact(m);    /* He^{(m)}_{m}(x) */
      double p_n1 = p_n0*(m+1)*x;      /* He^{(m)}_{m+1}(x) */
      double p_n = p_n1;
      int j;

      for (j = 0; j < m; j++)
        result_array[j] = 0.0;

      result_array[m] = p_n0;
      result_array[m + 1] = p_n1;

      for (j = m + 1; j <= nmax - 1; j++)
        {
          p_n  = (x*p_n1 - j*p_n0) * (j + 1.0) / (j - m + 1.0);
          p_n0 = p_n1;
          p_n1 = p_n;
          result_array[j + 1] = p_n;
        }

      return GSL_SUCCESS;
    }
}

/* Evaluates all derivatives (starting from 0) up to the mmax-th derivative of the probabilists' Hermite polynomial of order n at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */

int
gsl_sf_hermite_prob_deriv_array(const int mmax, const int n, const double x, double * result_array)
{
  if (n < 0 || mmax < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (n == 0)
    {
      int j;

      result_array[0] = 1.0;
      for (j = 1; j <= mmax; j++)
        result_array[j] = 0.0;

      return GSL_SUCCESS;
    }
  else if (n == 1 && mmax > 0)
    {
      int j;

      result_array[0] = x;
      result_array[1] = 1.0;

      for(j=2; j <= mmax; j++)
        result_array[j] = 0.0;

      return GSL_SUCCESS;
    }
  else if (mmax == 0)
    {
      result_array[0] = gsl_sf_hermite_prob(n,x);
      return GSL_SUCCESS;
    }
  else if (mmax == 1)
    {
      result_array[0] = gsl_sf_hermite_prob(n,x);
      result_array[1] = n*gsl_sf_hermite_prob(n-1,x);
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence */

      int k = GSL_MAX_INT(0, n - mmax);
      double p_n0 = gsl_sf_hermite_prob(k,x);        /* He_k(x) */
      double p_n1 = gsl_sf_hermite_prob(k+1,x);      /* He_{k+1}(x) */
      double p_n  = p_n1;
      int j;

      for(j=n+1; j <= mmax; j++)
        result_array[j] = 0.0;

      result_array[GSL_MIN_INT(n,mmax)] = p_n0;
      result_array[GSL_MIN_INT(n,mmax)-1] = p_n1;

      for (j = GSL_MIN_INT(mmax,n)-1; j > 0; j--)
        {
          k++;
          p_n  = x*p_n1-k*p_n0;
          p_n0 = p_n1;
          p_n1 = p_n;
          result_array[j - 1] = p_n;
        }

      p_n = 1.0;
      for (j = 1; j <= GSL_MIN_INT(n, mmax); j++)
        {
          p_n  = p_n*(n-j+1);
          result_array[j] = p_n*result_array[j];
        }

      return GSL_SUCCESS;
    }
}

/* Evaluates the series sum_{j=0}^n a_j*He_j(x) with He_j being the j-th probabilists' Hermite polynomial.
 * For improved numerical stability the Clenshaw algorithm (Clenshaw, C. W. (July 1955). "A note on the summation of Chebyshev series". Mathematical Tables and other Aids to Computation 9 (51): 118–110.) adapted to probabilists' Hermite polynomials is used. */

int
gsl_sf_hermite_prob_series_e(const int n, const double x, const double * a, gsl_sf_result * result)
{
  if(n < 0) {
    DOMAIN_ERROR(result);
  }
  else if(n == 0) {
    result->val = a[0];
    result->err = 0.;
    return GSL_SUCCESS;
  }
  else if(n == 1) {
    result->val = a[0]+a[1]*x;
    result->err = 2.*GSL_DBL_EPSILON * (fabs(a[0]) + fabs(a[1]*x)) ;
    return GSL_SUCCESS;
  }
  else {
    /* downward recurrence: b_n = a_n + x b_{n+1} - (n+1) b_{n+2} */

    double b0 = 0.;
    double b1 = 0.;
    double btmp = 0.;

    double e0 = 0.;
    double e1 = 0.;
    double etmp = e1;

    int j;

    for(j=n; j >= 0; j--){
      btmp = b0;
      b0  = a[j]+x*b0-(j+1)*b1;
      b1 = btmp;

      etmp = e0;
      e0  = (GSL_DBL_EPSILON*fabs(a[j])+fabs(x)*e0+(j+1)*e1);
      e1 = etmp;
    }

    result->val = b0;
    result->err = e0 + fabs(b0)*GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
}

double
gsl_sf_hermite_prob_series(const int n, const double x, const double * a)
{
  EVAL_RESULT(gsl_sf_hermite_prob_series_e(n, x, a, &result));
}

/* Evaluates all physicists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */
int
gsl_sf_hermite_array(const int nmax, const double x, double * result_array)
{
  if(nmax < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (nmax == 0)
    {
      result_array[0] = 1.0;
      return GSL_SUCCESS;
    }
  else if (nmax == 1)
    {
      result_array[0] = 1.0;
      result_array[1] = 2.0*x;
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence: H_{n+1} = 2x H_n - 2n H_{n-1} */

      int status = GSL_SUCCESS;
      const double two_x = 2.0 * x;
      const double abs_two_x = fabs(two_x);
      const double thresh1 = abs_two_x > 1.0 ? 0.9 * GSL_DBL_MAX / abs_two_x : GSL_DBL_MAX;
      const double thresh2 = 0.9 * GSL_DBL_MAX / 2.0;

      double p_n0 = 1.0;    /* H_0(x) */
      double p_n1 = two_x;  /* H_1(x) */
      double p_n = p_n1;

      int j;

      result_array[0] = 1.0;
      result_array[1] = 2.0*x;

      for (j = 1; j < nmax; j++)
        {
          if (fabs(p_n1) > thresh1 ||    /* test if 2*x*p_n1 will overflow */
              fabs(p_n0) > thresh2 / j)  /* test if 2*j*p_n0 will overflow */
            {
              status = GSL_EOVRFLW;
            }

          p_n  = two_x*p_n1 - 2.0*j*p_n0;
          p_n0 = p_n1;
          p_n1 = p_n;

          result_array[j + 1] = p_n;
        }

      return status;
    }
}

/* Evaluates the m-th derivative of all physicists' Hermite polynomials up to order nmax at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */
int
gsl_sf_hermite_array_deriv(const int m, const int nmax, const double x, double * result_array)
{
  if (nmax < 0 || m < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (m == 0)
    {
      gsl_sf_hermite_array(nmax, x, result_array);
      return GSL_SUCCESS;
    }
  else if (nmax < m)
    {
      int j;

      for(j = 0; j <= nmax; j++)
        result_array[j] = 0.0;

      return GSL_SUCCESS;
    }
  else if (nmax == m)
    {
      int j;

      for(j = 0; j < m; j++)
        result_array[j] = 0.0;

      result_array[nmax] = pow2(m)*gsl_sf_fact(m);

      return GSL_SUCCESS;
    }
  else if (nmax == m + 1)
    {
      int j;

      for(j = 0; j < m; j++)
        result_array[j] = 0.0;

      result_array[nmax-1] = pow2(m)*gsl_sf_fact(m);
      result_array[nmax] = result_array[nmax-1]*2*(m+1)*x;
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence: H^{(m)}_{n+1} = 2(n+1)/(n-m+1)*(x H^{(m)}_n - n H^{(m)}_{n-1}) */

      double p_n0 = pow2(m)*gsl_sf_fact(m);  /* H^{(m)}_{m}(x) */
      double p_n1 = p_n0*2*(m+1)*x;          /* H^{(m)}_{m+1}(x) */
      double p_n;
      int j;

      for(j = 0; j < m; j++)
        result_array[j] = 0.0;

      result_array[m] = p_n0;
      result_array[m+1] = p_n1;

      for (j = m + 1; j < nmax; ++j)
        {
          p_n  = (x*p_n1 - j*p_n0) * 2 * (j + 1.0) / (j - m + 1.0);
          p_n0 = p_n1;
          p_n1 = p_n;
          result_array[j + 1] = p_n;
        }

      return GSL_SUCCESS;
    }
}

/* Evaluates all derivatives (starting from 0) up to the mmax-th derivative of the physicists' Hermite polynomial of order n at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */
int
gsl_sf_hermite_deriv_array(const int mmax, const int n, const double x, double * result_array)
{
  if (n < 0 || mmax < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (n == 0)
    {
      int j;

      result_array[0] = 1.0;
      for(j = 1; j <= mmax; j++)
        result_array[j] = 0.0;

      return GSL_SUCCESS;
    }
  else if (n == 1 && mmax > 0)
    {
      int j;

      result_array[0] = 2*x;
      result_array[1] = 1.0;
      for (j = 2; j <= mmax; j++)
        result_array[j] = 0.0;

      return GSL_SUCCESS;
    }
  else if (mmax == 0)
    {
      result_array[0] = gsl_sf_hermite(n,x);
      return GSL_SUCCESS;
    }
  else if (mmax == 1)
    {
      result_array[0] = gsl_sf_hermite(n,x);
      result_array[1] = 2*n*gsl_sf_hermite(n - 1, x);
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence */

      int k = GSL_MAX_INT(0, n - mmax);
      double p_n0 = gsl_sf_hermite(k, x);        /* H_k(x) */
      double p_n1 = gsl_sf_hermite(k + 1, x);    /* H_{k+1}(x) */
      double p_n  = p_n1;
      int j;

      for (j = n + 1; j <= mmax; j++)
        result_array[j] = 0.0;

      result_array[GSL_MIN_INT(n,mmax)] = p_n0;
      result_array[GSL_MIN_INT(n,mmax)-1] = p_n1;

      for (j = GSL_MIN_INT(mmax, n) - 1; j > 0; j--)
        {
          k++;
          p_n  = 2*x*p_n1 - 2*k*p_n0;
          p_n0 = p_n1;
          p_n1 = p_n;
          result_array[j - 1] = p_n;
        }

      p_n = 1.0;
      for (j = 1; j <= GSL_MIN_INT(n,mmax); j++)
        {
          p_n *= 2.0 * (n - j + 1.0);
          result_array[j] *= p_n;
        }

      return GSL_SUCCESS;
    }
}

/* Evaluates the series sum_{j=0}^n a_j*H_j(x) with H_j being the j-th physicists' Hermite polynomial.
 * For improved numerical stability the Clenshaw algorithm (Clenshaw, C. W. (July 1955). "A note on the summation of Chebyshev series". Mathematical Tables and other Aids to Computation 9 (51): 118–110.) adapted to physicists' Hermite polynomials is used. */
int
gsl_sf_hermite_series_e(const int n, const double x, const double * a, gsl_sf_result * result)
{
  if(n < 0) {
    DOMAIN_ERROR(result);
  }
  else if(n == 0) {
    result->val = a[0];
    result->err = 0.;
    return GSL_SUCCESS;
  }
  else if(n == 1) {
    result->val = a[0]+a[1]*2.*x;
    result->err = 2.*GSL_DBL_EPSILON * (fabs(a[0]) + fabs(a[1]*2.*x)) ;
    return GSL_SUCCESS;
  }
  else {
    /* downward recurrence: b_n = a_n + 2x b_{n+1} - 2(n+1) b_{n+2} */

    double b0 = 0.;
    double b1 = 0.;
    double btmp = 0.;

    double e0 = 0.;
    double e1 = 0.;
    double etmp = e1;

    int j;

    for(j=n; j >= 0; j--){
      btmp = b0;
      b0  = a[j]+2.*x*b0-2.*(j+1)*b1;
      b1 = btmp;

      etmp = e0;
      e0  = (GSL_DBL_EPSILON*fabs(a[j])+fabs(2.*x)*e0+2.*(j+1)*e1);
      e1 = etmp;
    }

    result->val = b0;
    result->err = e0 + fabs(b0)*GSL_DBL_EPSILON;
    return GSL_SUCCESS;
  }
}

double
gsl_sf_hermite_series(const int n, const double x, const double * a)
{
  EVAL_RESULT(gsl_sf_hermite_series_e(n, x, a, &result));
}

/* Evaluates all Hermite functions up to order nmax at position x. The results are stored in result_array.
 * Since all polynomial orders are needed, upward recurrence is employed. */
int
gsl_sf_hermite_func_array(const int nmax, const double x, double * result_array)
{
  if (nmax < 0)
    {
      GSL_ERROR ("domain error", GSL_EDOM);
    }
  else if (nmax == 0)
    {
      result_array[0] = exp(-0.5*x*x)/sqrt(M_SQRTPI);
      return GSL_SUCCESS;
    }
  else if (nmax == 1)
    {
      result_array[0] = exp(-0.5*x*x)/sqrt(M_SQRTPI);
      result_array[1] = result_array[0]*M_SQRT2*x;
      return GSL_SUCCESS;
    }
  else
    {
      /* upward recurrence: Psi_{n+1} = sqrt(2/(n+1))*x Psi_n - sqrt(n/(n+1)) Psi_{n-1} */

      const double arg = -0.5 * x * x;
      double hi2 = 1.0 / sqrt(M_SQRTPI);
      double hi1 = M_SQRT2 * x * hi2;
      double hi = 0.0;
      double sum_log_scale = 0.0;
      double abshi;
      int i;

      result_array[0] = exp(arg) * hi2;
      result_array[1] = result_array[0] * M_SQRT2 * x;

      for (i = 2; i <= nmax; ++i)
        {
          hi = sqrt(2.0 / i) * x * hi1 - sqrt((i - 1.0) / i) * hi2;
          hi2 = hi1;
          hi1 = hi;

          abshi = fabs(hi);
          if (abshi > 1.0)
            {
              double log_scale = RND(log(abshi));
              double scale = exp(-log_scale);

              hi *= scale;
              hi1 *= scale;
              hi2 *= scale;
              sum_log_scale += log_scale;
            }

          result_array[i] = hi * exp(arg + sum_log_scale);
        }

      return GSL_SUCCESS;
    }
}

/* Evaluates the series sum_{j=0}^n a_j*Psi_j(x) with Psi_j being the j-th Hermite function.
 * For improved numerical stability the Clenshaw algorithm (Clenshaw, C. W. (July 1955). "A note on the summation of Chebyshev series". Mathematical Tables and other Aids to Computation 9 (51): 118–110.) adapted to Hermite functions is used. */

int
gsl_sf_hermite_func_series_e(const int n, const double x, const double * a, gsl_sf_result * result)
{
  if (n < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (n == 0)
    {
      result->val = a[0]*exp(-0.5*x*x)/sqrt(M_SQRTPI);
      result->err = GSL_DBL_EPSILON*fabs(result->val);
      return GSL_SUCCESS;
    }
  else if (n == 1)
    {
      result->val = (a[0]+a[1]*M_SQRT2*x)*exp(-0.5*x*x)/sqrt(M_SQRTPI);
      result->err = 2.*GSL_DBL_EPSILON*(fabs(a[0])+fabs(a[1]*M_SQRT2*x))*exp(-0.5*x*x)/sqrt(M_SQRTPI);
      return GSL_SUCCESS;
    }
  else
    {
      /* downward recurrence: b_n = a_n + sqrt(2/(n+1))*x b_{n+1} - sqrt((n+1)/(n+2)) b_{n+2} */

      double b0 = 0.;
      double b1 = 0.;
      double btmp = 0.;

      double e0 = 0.;
      double e1 = 0.;
      double etmp = e1;

      int j;

      for (j = n; j >= 0; j--)
        {
          btmp = b0;
          b0  = a[j]+sqrt(2./(j+1))*x*b0-sqrt((j+1.)/(j+2.))*b1;
          b1 = btmp;

          etmp = e0;
          e0  = (GSL_DBL_EPSILON*fabs(a[j])+sqrt(2./(j+1))*fabs(x)*e0+sqrt((j+1.)/(j+2.))*e1);
          e1 = etmp;
        }

      result->val = b0*exp(-0.5*x*x)/sqrt(M_SQRTPI);
      result->err = e0 + fabs(result->val)*GSL_DBL_EPSILON;

      return GSL_SUCCESS;
    }
}

double
gsl_sf_hermite_func_series(const int n, const double x, const double * a)
{
  EVAL_RESULT(gsl_sf_hermite_func_series_e(n, x, a, &result));
}

/* Evaluates the m-th derivative of the Hermite function of order n at position x.
 * A summation including upward recurrences is used. */

int
gsl_sf_hermite_func_der_e(const int m, const int n, const double x, gsl_sf_result * result)
{
  if(m < 0 || n < 0)
    {
      DOMAIN_ERROR(result);
    }
  else if (m == 0)
    {
      return gsl_sf_hermite_func_e(n, x, result);
    }
  else if (m == 1)
    {
      double hi2 = 1.0 / sqrt(M_SQRTPI);
      double hi1 = M_SQRT2 * x * hi2;
      double hi = 0.0;
      double sum_log_scale = 0.0;
      double abshi;
      int i;

      for (i = 2; i <= n; ++i)
        {
          hi = sqrt(2.0 / i) * x * hi1 - sqrt((i - 1.0) / i) * hi2;
          hi2 = hi1;
          hi1 = hi;

          abshi = fabs(hi);
          if (abshi > 1.0)
            {
              double log_scale = RND(log(abshi));
              double scale = exp(-log_scale);

              hi *= scale;
              hi1 *= scale;
              hi2 *= scale;
              sum_log_scale += log_scale;
            }
        }

      /* psi'_n(x) = sqrt(2 n) psi_{n-1} - x psi_n */
      result->val = (sqrt(2.0*n) * hi2 - x * hi) * exp(-0.5 * x * x + sum_log_scale);
      result->err = n * GSL_DBL_EPSILON * fabs(result->val);

      return GSL_SUCCESS;
    }
  else
    {
      int j;
      double r,er,b;
      double h0 = 1.;
      double h1 = x;
      double eh0 = GSL_DBL_EPSILON;
      double eh1 = GSL_DBL_EPSILON;
      double p0 = 1.;
      double p1 = M_SQRT2*x;
      double ep0 = GSL_DBL_EPSILON;
      double ep1 = M_SQRT2*GSL_DBL_EPSILON;
      double f = 1.;

      for (j=GSL_MAX_INT(1,n-m+1);j<=n;j++)
        f *= sqrt(2.*j);

      if (m > n)
        {
          f = (GSL_IS_ODD(m-n)?-f:f);

          for (j=0;j<GSL_MIN_INT(n,m-n);j++)
            f *= (m-j)/(j+1.);
        }

      for (j=1;j<=m-n;j++)
        {
          b = x*h1-j*h0;
          h0 = h1;
          h1 = b;

          b  = (fabs(x)*eh1+j*eh0);
          eh0 = eh1;
          eh1 = b;
        }

      b = 0.;
      for (j=1;j<=n-m;j++)
        {
          b = (M_SQRT2*x*p1-sqrt(j)*p0)/sqrt(j+1.);
          p0 = p1;
          p1 = b;

          b  = (M_SQRT2*fabs(x)*ep1+sqrt(j)*ep0)/sqrt(j+1.);
          ep0 = ep1;
          ep1 = b;
        }

      b = 0.;
      r = 0.;
      er = 0.;
      for (j=GSL_MAX_INT(0,m-n);j<=m;j++)
        {
          r += f*h0*p0;
          er += eh0*fabs(f*p0)+ep0*fabs(f*h0)+GSL_DBL_EPSILON*fabs(f*h0*p0);

          b = x*h1-(j+1.)*h0;
          h0 = h1;
          h1 = b;

          b  = 0.5*(fabs(x)*eh1+(j+1.)*eh0);
          eh0 = eh1;
          eh1 = b;

          b = (M_SQRT2*x*p1-sqrt(n-m+j+1.)*p0)/sqrt(n-m+j+2.);
          p0 = p1;
          p1 = b;

          b  = 0.5*(M_SQRT2*fabs(x)*ep1+sqrt(n-m+j+1.)*ep0)/sqrt(n-m+j+2.);
          ep0 = ep1;
          ep1 = b;

          f *= -(m-j)/(j+1.)/sqrt(n-m+j+1.)*M_SQRT1_2;
        }

      result->val = r*exp(-0.5*x*x)/sqrt(M_SQRTPI);
      result->err = er*fabs(exp(-0.5*x*x)/sqrt(M_SQRTPI)) + GSL_DBL_EPSILON*fabs(result->val);

      return GSL_SUCCESS;
    }
}

double
gsl_sf_hermite_func_der(const int m, const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_func_der_e(m, n, x, &result));
}

static double
H_zero_init(const int n, const int k)
{
  double p = 1., x = 1., y = 1.;
  if (k == 1 && n > 50) {
    x = (GSL_IS_ODD(n)?1./sqrt((n-1)/6.):1./sqrt(0.5*n));
  }
  else {
    p = -0.7937005259840997373758528196*gsl_sf_airy_zero_Ai(n/2-k+1);
    x = sqrt(2*n+1.);
    y = pow(2*n+1.,1/6.);
    x = x - p/y - 0.1*p*p/(x*y*y) + (9/280. - p*p*p*11/350.)/(x*x*x) + (p*277/12600. - gsl_sf_pow_int(p,4)*823/63000.)/gsl_sf_pow_int(x,4)/y;
  }
  p = acos(x/sqrt(2*n+1.));
  y = M_PI*(-2*(n/2-k)-1.5)/(n+0.5);
  if(gsl_fcmp(y,sin(2.*p)-2*p,GSL_SQRT_DBL_EPSILON)==0) return x; /* initial approx sufficiently accurate */
  if (y > -GSL_DBL_EPSILON) return sqrt(2*n+1.);
  if (p < GSL_DBL_EPSILON) p = GSL_DBL_EPSILON;
  if (p > M_PI_2) p = M_PI_2;
  if (sin(2.*p)-2*p > y){
    x = GSL_MAX((sin(2.*p)-2*p-y)/4.,GSL_SQRT_DBL_EPSILON);
    do{
      x *= 2.;
      p += x;
    } while (sin(2.*p)-2*p > y);
  }
  do {
    x = p;
    p -= (sin(2.*p)-2.*p-y)/(2.*cos(2.*p)-2.);
    if (p<0.||p>M_PI_2) p = M_PI_2;
  } while (gsl_fcmp(x,p,100*GSL_DBL_EPSILON)!=0);
  return sqrt(2*n+1.)*cos(p);
}

/* lookup table for the positive zeros of the probabilists' Hermite polynomials of order 3 through 20 */
static double He_zero_tab[99] = {
  1.73205080756887729352744634151,
  0.741963784302725857648513596726,
  2.33441421833897723931751226721,
  1.35562617997426586583052129087,
  2.85697001387280565416230426401,
  0.616706590192594152193686099399,
  1.88917587775371067550566789858,
  3.32425743355211895236183546247,
  1.154405394739968127239597758838,
  2.36675941073454128861885646856,
  3.75043971772574225630392202571,
  0.539079811351375108072461918694,
  1.63651904243510799922544657297,
  2.80248586128754169911301080618,
  4.14454718612589433206019783917,
  1.023255663789132524828148225810,
  2.07684797867783010652215614374,
  3.20542900285646994336567590292,
  4.51274586339978266756667884317,
  0.484935707515497653046233483105,
  1.46598909439115818325066466416,
  2.48432584163895458087625118368,
  3.58182348355192692277623675546,
  4.85946282833231215015516494660,
  0.928868997381063940144111999584,
  1.87603502015484584534137013967,
  2.86512316064364499771968407254,
  3.93616660712997692868589612142,
  5.18800122437487094818666404539,
  0.444403001944138945299732445510,
  1.34037519715161672153112945211,
  2.25946445100079912386492979448,
  3.22370982877009747166319001956,
  4.27182584793228172295999293076,
  5.50090170446774760081221630899,
  0.856679493519450033897376121795,
  1.72541837958823916151095838741,
  2.62068997343221478063807762201,
  3.56344438028163409162493844661,
  4.59139844893652062705231872720,
  5.80016725238650030586450565322,
  0.412590457954601838167454145167,
  1.24268895548546417895063983219,
  2.08834474570194417097139675101,
  2.96303657983866750254927123447,
  3.88692457505976938384755016476,
  4.89693639734556468372449782879,
  6.08740954690129132226890147034,
  0.799129068324547999424888414207,
  1.60671006902872973652322479373,
  2.43243682700975804116311571682,
  3.28908242439876638890856229770,
  4.19620771126901565957404160583,
  5.19009359130478119946445431715,
  6.36394788882983831771116094427,
  0.386760604500557347721047189801,
  1.16382910055496477419336819907,
  1.95198034571633346449212362880,
  2.76024504763070161684598142269,
  3.60087362417154828824902745506,
  4.49295530252001124266582263095,
  5.47222570594934308841242925805,
  6.63087819839312848022981922233,
  0.751842600703896170737870774614,
  1.50988330779674075905491513417,
  2.28101944025298889535537879396,
  3.07379717532819355851658337833,
  3.90006571719800990903311840097,
  4.77853158962998382710540812497,
  5.74446007865940618125547815768,
  6.88912243989533223256205432938,
  0.365245755507697595916901619097,
  1.09839551809150122773848360538,
  1.83977992150864548966395498992,
  2.59583368891124032910545091458,
  3.37473653577809099529779309480,
  4.18802023162940370448450911428,
  5.05407268544273984538327527397,
  6.00774591135959752029303858752,
  7.13946484914647887560975631213,
  0.712085044042379940413609979021,
  1.42887667607837287134157901452,
  2.15550276131693514033871248449,
  2.89805127651575312007902775275,
  3.66441654745063847665304033851,
  4.46587262683103133615452574019,
  5.32053637733603803162823765939,
  6.26289115651325170419416064557,
  7.38257902403043186766326977122,
  0.346964157081355927973322447164,
  1.04294534880275103146136681143,
  1.74524732081412671493067861704,
  2.45866361117236775131735057433,
  3.18901481655338941485371744116,
  3.94396735065731626033176813604,
  4.73458133404605534390170946748,
  5.57873880589320115268040332802,
  6.51059015701365448636289263918,
  7.61904854167975829138128156060
};

/*
 * Computes the s-th zero the probabilists' Hermite polynomial of order n.
 * A Newton iteration using a continued fraction representation adapted from:
 *
 * [E.T. Whittaker (1914), On the continued fractions which represent the
 * functions of Hermite and other functions defined by differential equations,
 * Proceedings of the Edinburgh Mathematical Society, 32, 65-74]
 *
 * is performed with the initial approximation from
 *
 * [Arpad Elbert and Martin E. Muldoon, Approximations for zeros of Hermite
 * functions, pp. 117-126 in D. Dominici and R. S. Maier, eds, "Special Functions
 * and Orthogonal Polynomials", Contemporary Mathematics, vol 471 (2008)]
 *
 * refined via the bisection method.
 */

int
gsl_sf_hermite_prob_zero_e(const int n, const int s, gsl_sf_result * result)
{
  if (n <= 0 || s < 0 || s > n/2)
    {
      DOMAIN_ERROR(result);
    }
  else if (s == 0)
    {
      if (GSL_IS_ODD(n) == 1)
        {
          result->val = 0.;
          result->err = 0.;
          return GSL_SUCCESS;
        }
      else
        {
          DOMAIN_ERROR(result);
        }
    }
  else if (n == 2)
    {
      result->val = 1.;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else if (n < 21)
    {
      result->val = He_zero_tab[(GSL_IS_ODD(n)?n/2:0)+((n/2)*(n/2-1))+s-2];
      result->err = GSL_DBL_EPSILON*(result->val);
      return GSL_SUCCESS;
    }
  else
    {
      double d = 1., x = 1., x0 = 1.;
      int j;
      x = H_zero_init(n,s) * M_SQRT2;
      do
        {
          x0 = x;
          d = 0.;
          for (j=1; j<n; j++)
            d = j/(x-d);

          x -= (x-d)/n;

          /* gsl_fcmp can be used since the smallest zero approaches 1/sqrt(n) or 1/sqrt((n-1)/3.)
           * for large n and thus all zeros are non-zero (except for the trivial case handled above) */
        }
      while (gsl_fcmp(x, x0, 10*GSL_DBL_EPSILON) != 0);

      result->val = x;
      result->err = 2*GSL_DBL_EPSILON*x + fabs(x-x0);

      return GSL_SUCCESS;
    }
}

double
gsl_sf_hermite_prob_zero(const int n, const int s)
{
  EVAL_RESULT(gsl_sf_hermite_prob_zero_e(n, s, &result));
}

/* lookup table for the positive zeros of the physicists' Hermite polynomials of order 3 through 20 */
static double H_zero_tab[99] = {
  1.22474487139158904909864203735,
  0.524647623275290317884060253835,
  1.65068012388578455588334111112,
  0.958572464613818507112770593893,
  2.02018287045608563292872408814,
  0.436077411927616508679215948251,
  1.335849074013696949714895282970,
  2.35060497367449222283392198706,
  0.816287882858964663038710959027,
  1.67355162876747144503180139830,
  2.65196135683523349244708200652,
  0.381186990207322116854718885584,
  1.157193712446780194720765779063,
  1.98165675669584292585463063977,
  2.93063742025724401922350270524,
  0.723551018752837573322639864579,
  1.46855328921666793166701573925,
  2.26658058453184311180209693284,
  3.19099320178152760723004779538,
  0.342901327223704608789165025557,
  1.03661082978951365417749191676,
  1.75668364929988177345140122011,
  2.53273167423278979640896079775,
  3.43615911883773760332672549432,
  0.656809566882099765024611575383,
  1.32655708449493285594973473558,
  2.02594801582575533516591283121,
  2.78329009978165177083671870152,
  3.66847084655958251845837146485,
  0.314240376254359111276611634095,
  0.947788391240163743704578131060,
  1.59768263515260479670966277090,
  2.27950708050105990018772856942,
  3.02063702512088977171067937518,
  3.88972489786978191927164274724,
  0.605763879171060113080537108602,
  1.22005503659074842622205526637,
  1.85310765160151214200350644316,
  2.51973568567823788343040913628,
  3.24660897837240998812205115236,
  4.10133759617863964117891508007,
  0.291745510672562078446113075799,
  0.878713787329399416114679311861,
  1.47668273114114087058350654421,
  2.09518325850771681573497272630,
  2.74847072498540256862499852415,
  3.46265693360227055020891736115,
  4.30444857047363181262129810037,
  0.565069583255575748526020337198,
  1.13611558521092066631913490556,
  1.71999257518648893241583152515,
  2.32573248617385774545404479449,
  2.96716692790560324848896036355,
  3.66995037340445253472922383312,
  4.49999070730939155366438053053,
  0.273481046138152452158280401965,
  0.822951449144655892582454496734,
  1.38025853919888079637208966969,
  1.95178799091625397743465541496,
  2.54620215784748136215932870545,
  3.17699916197995602681399455926,
  3.86944790486012269871942409801,
  4.68873893930581836468849864875,
  0.531633001342654731349086553718,
  1.06764872574345055363045773799,
  1.61292431422123133311288254454,
  2.17350282666662081927537907149,
  2.75776291570388873092640349574,
  3.37893209114149408338327069289,
  4.06194667587547430689245559698,
  4.87134519367440308834927655662,
  0.258267750519096759258116098711,
  0.776682919267411661316659462284,
  1.30092085838961736566626555439,
  1.83553160426162889225383944409,
  2.38629908916668600026459301424,
  2.96137750553160684477863254906,
  3.57376906848626607950067599377,
  4.24811787356812646302342016090,
  5.04836400887446676837203757885,
  0.503520163423888209373811765050,
  1.01036838713431135136859873726,
  1.52417061939353303183354859367,
  2.04923170985061937575050838669,
  2.59113378979454256492128084112,
  3.15784881834760228184318034120,
  3.76218735196402009751489394104,
  4.42853280660377943723498532226,
  5.22027169053748216460967142500,
  0.245340708300901249903836530634,
  0.737473728545394358705605144252,
  1.23407621539532300788581834696,
  1.73853771211658620678086566214,
  2.25497400208927552308233334473,
  2.78880605842813048052503375640,
  3.34785456738321632691492452300,
  3.94476404011562521037562880052,
  4.60368244955074427307767524898,
  5.38748089001123286201690041068
};

/*
 * Computes the s-th zero the physicists' Hermite polynomial of order n, thus also
 * the s-th zero of the Hermite function of order n. A Newton iteration using a continued
 * fraction representation adapted from:
 *
 * [E.T. Whittaker (1914), On the continued fractions which represent the functions of Hermite
 * and other functions defined by differential equations, Proceedings of the Edinburgh Mathematical
 * Society, 32, 65-74]
 *
 * An initial approximation is used from:
 *
 * [Arpad Elbert and Martin E. Muldoon, Approximations for zeros of Hermite functions,
 * pp. 117-126 in D. Dominici and R. S. Maier, eds, "Special Functions and Orthogonal Polynomials",
 * Contemporary Mathematics, vol 471 (2008)]
 *
 * which is refined via the bisection method.
 */

int
gsl_sf_hermite_zero_e(const int n, const int s, gsl_sf_result * result)
{
  if (n <= 0 || s < 0 || s > n/2)
    {
      DOMAIN_ERROR(result);
    }
  else if (s == 0)
    {
      if (GSL_IS_ODD(n) == 1)
        {
          result->val = 0.;
          result->err = 0.;
          return GSL_SUCCESS;
        }
      else
        {
          DOMAIN_ERROR(result);
        }
    }
  else if (n == 2)
    {
      result->val = M_SQRT1_2;
      result->err = 0.;
      return GSL_SUCCESS;
    }
  else if (n < 21)
    {
      result->val = H_zero_tab[(GSL_IS_ODD(n)?n/2:0)+((n/2)*(n/2-1))+s-2];
      result->err = GSL_DBL_EPSILON*(result->val);
      return GSL_SUCCESS;
    }
  else
    {
      double d = 1., x = 1., x0 = 1.;
      int j;

      x = H_zero_init(n,s);
      do
        {
          x0 = x;
          d = 0.;

          for (j=1; j<n; j++)
            d = 2*j/(2.*x-d);

          x -= (2*x-d)*0.5/n;

        /* gsl_fcmp can be used since the smallest zero approaches 1/sqrt(n) or 1/sqrt((n-1)/3.)
         * for large n and thus all zeros are non-zero (except for the trivial case handled above) */
        }
      while (gsl_fcmp(x, x0, 10*GSL_DBL_EPSILON) != 0);

      result->val = x;
      result->err = 2*GSL_DBL_EPSILON*x + fabs(x-x0);

      return GSL_SUCCESS;
    }
}

double
gsl_sf_hermite_zero(const int n, const int s)
{
  EVAL_RESULT(gsl_sf_hermite_zero_e(n, s, &result));
}

int
gsl_sf_hermite_func_zero_e(const int n, const int s, gsl_sf_result * result)
{
  return gsl_sf_hermite_zero_e(n, s, result);
}

double
gsl_sf_hermite_func_zero(const int n, const int s)
{
  EVAL_RESULT(gsl_sf_hermite_func_zero_e(n, s, &result));
}

#ifndef GSL_DISABLE_DEPRECATED

int
gsl_sf_hermite_phys_e(const int n, const double x, gsl_sf_result * result)
{
  return gsl_sf_hermite_e(n, x, result);
}

double
gsl_sf_hermite_phys(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_phys_e(n, x, &result));
}

int
gsl_sf_hermite_phys_der_e(const int m, const int n, const double x, gsl_sf_result * result)
{
  return gsl_sf_hermite_deriv_e(m, n, x, result);
}

double
gsl_sf_hermite_phys_der(const int m, const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_phys_der_e(m, n, x, &result));
}

int
gsl_sf_hermite_phys_array(const int nmax, const double x, double * result_array)
{
  return gsl_sf_hermite_array(nmax, x, result_array);
}

int
gsl_sf_hermite_phys_series_e(const int n, const double x, const double * a, gsl_sf_result * result)
{
  return gsl_sf_hermite_series_e(n, x, a, result);
}

double
gsl_sf_hermite_phys_series(const int n, const double x, const double * a)
{
  EVAL_RESULT(gsl_sf_hermite_phys_series_e(n, x, a, &result));
}

int
gsl_sf_hermite_phys_array_der(const int m, const int nmax, const double x, double * result_array)
{
  return gsl_sf_hermite_array_deriv(m, nmax, x, result_array);
}

int
gsl_sf_hermite_phys_der_array(const int mmax, const int n, const double x, double * result_array)
{
  return gsl_sf_hermite_deriv_array(mmax, n, x, result_array);
}

int
gsl_sf_hermite_phys_zero_e(const int n, const int s, gsl_sf_result * result)
{
  return gsl_sf_hermite_zero_e(n, s, result);
}

double
gsl_sf_hermite_phys_zero(const int n, const int s)
{
  EVAL_RESULT(gsl_sf_hermite_zero_e(n, s, &result));
}

int
gsl_sf_hermite_prob_array_der(const int m, const int nmax, const double x, double * result_array)
{
  return gsl_sf_hermite_prob_array_deriv(m, nmax, x, result_array);
}

int
gsl_sf_hermite_prob_der_array(const int mmax, const int n, const double x, double * result_array)
{
  return gsl_sf_hermite_prob_deriv_array(mmax, n, x, result_array);
}

int
gsl_sf_hermite_prob_der_e(const int m, const int n, const double x, gsl_sf_result * result)
{
  return gsl_sf_hermite_prob_deriv_e(m, n, x, result);
}

double
gsl_sf_hermite_prob_der(const int m, const int n, const double x)
{
  EVAL_RESULT(gsl_sf_hermite_prob_deriv_e(m, n, x, &result));
}

#endif /* !GSL_DISABLE_DEPRECATED */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 * @file bblas_daccuracy.c
 *
 * @brief BBLAS accuracy tests for double_Complex precision.
 *
 *  BBLAS is a software package provided by Univ. of Manchester,
 *  Univ. of Tennessee.
 *
 * @version 1.0.0
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date    2016-02-20
 *
 **/

#ifndef DOXYGEN_SHOULD_SKIP_THIS
/**
 * Code generation
 * @generated from bblas_zaccuracy.c normal z -> d, Mon Jun  6 09:44:14 2016
 **/
#endif

#include "bblas_common.h"
#include <lapacke.h>
#include <cblas.h>

#define REAL

/**
 * bblas_dgemm_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|B\|_{\infty} +
 *           |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_dgemm_tolerance(bblas_dtest_t *test, double **C_final)
{

    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, Bm, Bn, M, K, N, ldc;
    int max_work_size = max(test->maxK,max(test->maxM, test->maxN));
    double Anorm, Bnorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    double alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    double **C_diff;
    C_diff = (double **) malloc(batch_count*sizeof(double *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_dcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    M = test->M[batch_iter];
	    K   = test->K[batch_iter];
	    N   = test->N[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if (test->transA[batch_iter] == BblasNoTrans) {
		Am = M; An = K;
	    } else {
		Am = K; An = M;
	    }
	    if (test->transB[batch_iter] == BblasNoTrans) {
		Bm = K; Bn = N;
	    } else {
		Bm = N; Bn = K;
	    }

	    /*Compute the error C - C_finial */
	    cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', M, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    result = Rnorm / ((K*cabs(test->alpha[batch_iter])*Anorm*Bnorm + cabs(test->beta[batch_iter])*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	M   = test->M[first_index];
	K   = test->K[first_index];
	ldc = test->ldc[first_index];
	N   = test->N[first_index];

	if (test->transA[first_index] == BblasNoTrans) {
	    Am = M; An = K;
	} else {
	    Am = K; An = M;
	}
	if (test->transB[first_index] == BblasNoTrans) {
	    Bm = K; Bn = N;
	} else {
	    Bm = N; Bn = K;
	}

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
		/*Compute the error C - C_finial */
	    cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);
	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', M, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    /*Compute the relative error */
	    result = Rnorm / ((K*cabs(test->alpha[first_index])*Anorm*Bnorm + cabs(test->beta[first_index])*test->Cinitnorm[batch_iter])*eps);

	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal: Fixed, Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }else
    {
	bblas_error("bblas_dtesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_dstarmm_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|B\|_{\infty} +
 * |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_dstarmm_tolerance(bblas_dtest_t *test, double **C_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, Bm, Bn, N, ldc;
    int max_work_size = max(test->maxM, test->maxN);
    double Anorm, Bnorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    double alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    double **C_diff;
    C_diff = (double **) malloc(batch_count*sizeof(double *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_dcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if(test->side[batch_iter] == BblasLeft )
	    {
		Am = test->M[batch_iter];
	    } else {
		Am = test->N[batch_iter];
	    }
	    Bm = test->ldb[batch_iter];
	    Bn = test->N[batch_iter];
	    /*Compute the error C - C_finial */
	    cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, Am, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    result = Rnorm / ((N*cabs(test->alpha[batch_iter])*Anorm*Bnorm + cabs(test->beta[batch_iter])*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	N = test->N[first_index];
	ldc = test->ldc[first_index];

	if(test->side[first_index] == BblasLeft )
	{
	    Am = test->M[first_index];
	} else {
	    Am = test->N[first_index];
	}
	Bm = test->ldb[first_index];
	Bn = test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    /*Compute the error C - C_finial */
	    cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, Am, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    /*Compute the relative error */
	    result = Rnorm / ((N*cabs(test->alpha[first_index])*Anorm*Bnorm + cabs(test->beta[first_index])*test->Cinitnorm[batch_iter])*eps);

	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal: Fixed, Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }else
    {
	bblas_error("bblas_dtesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_dstark_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|A\|_1 +
 * |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_dstark_tolerance(bblas_dtest_t *test, double **C_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, N, K, ldc;
    char u;
    int max_work_size = max(test->maxK, test->maxN);
    double Anorm, ATnorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    double alpha_norm, beta_norm;
    double alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    double **C_diff;
    C_diff = (double **) malloc(batch_count*sizeof(double *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_dcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    K = test->K[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if (test->trans[batch_iter] == BblasNoTrans) {
		Am = N; An = K;
	    } else {
		Am = K; An = N;
	    }

	    if(test->uplo[batch_iter] == BblasLower )
	      {
		u = 'L';
	      }else
	      {
		u = 'U';
	      }
	    /*Compute the error C - C_finial */
	    cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlansy_work(LAPACK_COL_MAJOR, 'I', u, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the one norm of A */
	    ATnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

	    if (test->routine == BBLAS_HERK)
	    {
		beta_norm   = cabs(test->beta_herk[batch_iter]);
		alpha_norm  = cabs(test->alpha_herk[batch_iter]);
	    }else{
		beta_norm  = cabs(test->beta[batch_iter]);
		alpha_norm  = cabs(test->alpha[batch_iter]);
	    }

	    result = Rnorm / ((K*alpha_norm*Anorm*ATnorm + beta_norm*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	    N = test->N[first_index];
	    K = test->K[first_index];
	    ldc = test->ldc[first_index];

	    if (test->trans[first_index] == BblasNoTrans) {
		Am = N; An = K;
	    } else {
		Am = K; An = N;
	    }
	    if(test->uplo[first_index] == BblasLower )
	      {
		u = 'L';
	      }else
	      {
		u = 'U';
	      }

	    for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	      {
		/*Compute the error C - C_finial */
		cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		Rnorm = LAPACKE_dlansy_work(LAPACK_COL_MAJOR, 'I', u, N, C_diff[batch_iter], ldc, work);

		/*Compute the infinity norm of A */
		Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

		/*Compute the one norm of A */
		ATnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

		if (test->routine == BBLAS_HERK)
		{
		    beta_norm   = cabs(test->beta_herk[first_index]);
		    alpha_norm  = cabs(test->alpha_herk[first_index]);
		}else{
		    beta_norm   = cabs(test->beta[first_index]);
		    alpha_norm  = cabs(test->alpha[first_index]);
		}

		/*Compute the relative error */
		result = Rnorm / ((K*alpha_norm*Anorm*ATnorm +
				   beta_norm*test->Cinitnorm[batch_iter])*eps);

		switch(test->target)
		{
		    case BBLAS_MKL:
			test->mkl_error[batch_iter] = result;
			break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		  test->device_error[batch_iter] = result;
		  break;

		    case BBLAS_OTHER:
			test->other_error[batch_iter] = result;
			break;
		    default:
			printf("In bblas_dcheck_Cfinal: Fixed, Target no defined\n");
			exit(EXIT_FAILURE);
		}
	    }
    }else
    {
	bblas_error("bblas_dtesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_dstar2k_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((K|\alpha|\|A\|_{\infty}\|B\|_1 +
 * K|\alpha|\|B\|_{\infty}\|A\|_1  +
 * |\beta|\|C_{\mathrm{init}}\|_{\infty})\epsilon).@f]
 *
 **/

void bblas_dstar2k_tolerance(bblas_dtest_t *test, double **C_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, Bn, Bm, N, K, ldc;
    int max_work_size = max(test->maxK, test->maxN);
    double Anorm, Bnorm, ATnorm, BTnorm,  Rnorm, result;
    double beta_norm;
    double eps = LAPACKE_dlamch_work('e');
    double alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayC -C_final) */
    double **C_diff;
    C_diff = (double **) malloc(batch_count*sizeof(double *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_dcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    K = test->K[batch_iter];
	    ldc = test->ldc[batch_iter];

	    if (test->trans[batch_iter] == BblasNoTrans) {
		Am = N; An = K;
		Bn = K;
	    } else {
		Am = K; An = N;
		Bn = N;
	    }
	    Bm = test->ldb[batch_iter];

	    /*Compute the error C - C_finial */
	    cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

	    /*Compute the infinity norm of A */
	    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the one norm of A */
	    ATnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

	    /*Compute the infinity norm of B */
	    Bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    /*Compute the one norm of B */
	    BTnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'O', Bm, Bn, test->arrayB[batch_iter], Bm, work);

	    if(test->routine == BBLAS_HER2K)
	    {
		beta_norm = cabs(test->beta_herk[batch_iter]);
	    }else
	    {
		beta_norm = cabs(test->beta[batch_iter]);
	    }

	    result = Rnorm / ((K*cabs(test->alpha[batch_iter])*Anorm*BTnorm +
			       K*cabs(test->alpha[batch_iter])*Bnorm*ATnorm +
			       beta_norm*test->Cinitnorm[batch_iter])*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	N = test->N[first_index];
	K = test->K[first_index];
	ldc = test->ldc[first_index];

	if (test->trans[first_index] == BblasNoTrans) {
		Am = N; An = K;
		Bn = K;
	} else {
	    Am = K; An = N;
	    Bn = N;
	}
	Bm = test->ldb[first_index];

	    for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	    {
		/*Compute the error C - C_finial */
		cblas_daxpy (ldc*N, (alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
		Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldc, N, C_diff[batch_iter], ldc, work);

		/*Compute the infinity norm of A */
		Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Am, An, test->arrayA[batch_iter], Am, work);

		/*Compute the one norm of A */
		ATnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'O', Am, An, test->arrayA[batch_iter], Am, work);

		/*Compute the infinity norm of B */
		Bnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', Bm, Bn, test->arrayB[batch_iter], Bm, work);

		/*Compute the one norm of B */
		BTnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'O', Bm, Bn, test->arrayB[batch_iter], Bm, work);

		if(test->routine == BBLAS_HER2K)
		{
		    beta_norm = cabs(test->beta_herk[first_index]);
		}else
		{
		    beta_norm = cabs(test->beta[first_index]);
		}

		/*Compute the relative error */
		result = Rnorm / ((K*cabs(test->alpha[first_index])*Anorm*BTnorm +
				   K*cabs(test->alpha[first_index])*Bnorm*ATnorm +
				   beta_norm*test->Cinitnorm[batch_iter])*eps);

		switch(test->target)
		{
		case BBLAS_MKL:
		  test->mkl_error[batch_iter] = result;
		  break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		  test->device_error[batch_iter] = result;
		  break;

		case BBLAS_OTHER:
		  test->other_error[batch_iter] = result;
		  break;
		default:
		  printf("In bblas_dcheck_Cfinal: Fixed, Target no defined\n");
		  exit(EXIT_FAILURE);
		}
	    }
    }else
    {
	bblas_error("bblas_dtesting.c", "wrong batch_opts value");
    }

/*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(C_diff[index]);
    }
    free(C_diff);
    free(work);
}

/**
 * bblas_dtrmm_tolerance computes the forward error
 * @f[\|C - C_{\mathrm{seq}}\|_{\infty}/((|\alpha|\|A\|_{\infty}\|B\|_{\infty}N\epsilon).@f]
 *
 **/

void bblas_dtrmm_tolerance(bblas_dtest_t *test, double **B_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, ldb, N, M;
    int max_work_size = max(test->maxK, test->maxN);
    char u, diag;
    double Anorm, Rnorm, result;
    double eps = LAPACKE_dlamch_work('e');
    double alpha =-1;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Temporary buffer to save (test-arrayB -B_final) */
    double **B_diff;
    B_diff = (double **) malloc(batch_count*sizeof(double *));

    /*Make a copy of test->arrayB in B_diff */
    bblas_dcopy_Binit(test, B_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    M = test->M[batch_iter];
	    ldb = test->ldb[batch_iter];

	    Am = test->lda[batch_iter];

	    if(test->side[batch_iter] == BblasLeft )
	    {
		An = M;
	    }else
	    {
		An = N;
	    }

	    if(test->uplo[batch_iter] == BblasLower )
	    {
	       u = 'L';
	    }else
	    {
	       u = 'U';
	    }

	    if(test->diag[batch_iter] == BblasUnit )
	    {
		diag = 'U';
	    }else
	    {
		diag = 'N';
	    }

	    /*Compute the error B - B_finial */
	    cblas_daxpy (ldb*N, (alpha), B_final[batch_iter], 1, B_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, B_diff[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_dlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

		result = Rnorm/(N*cabs(test->alpha[batch_iter])*Anorm*test->Binitnorm[batch_iter]*eps);
	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }

	}
    }else  if( test->batch_opts == BBLAS_FIXED )
    {

	N = test->N[first_index];
	M = test->M[first_index];
	ldb = test->ldb[first_index];

	Am = test->lda[first_index];

	if(test->side[first_index] == BblasLeft )
	{
	    An = M;
	}else
	{
	    An = N;
	}

	if(test->uplo[first_index] == BblasLower )
	{
	    u = 'L';
	}else
	{
	    u = 'U';
	}

	if(test->diag[first_index] == BblasUnit )
	{
	    diag = 'U';
	}else
	{
	    diag = 'N';
	}

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    /*Compute the error B - B_finial */
	    cblas_daxpy (ldb*N, (alpha), B_final[batch_iter], 1, B_diff[batch_iter], 1);

		/*Compute  the infinity norm associated  with the error */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, B_diff[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_dlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

	    result = Rnorm/(N*cabs(test->alpha[first_index])*Anorm*test->Binitnorm[batch_iter]*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }
    else
    {
	bblas_error("bblas_dtesting.c", "wrong batch_opts value");
    }

    /*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(B_diff[index]);
    }
    free(B_diff);
    free(work);
}

/**
 * bblas_dtrsm_tolerance computes the backward error
 * @f[\|\alpha B - AX \|_{\infty}/((\|A\|_{\infty}\|X\| + |\alpha|\|B\|_{\infty})N\epsilon).@f]
 *
 **/

void bblas_dtrsm_tolerance(bblas_dtest_t *test, double **B_final)
{
    /*Local variables */
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int Am, An, ldb, N, M;
    int max_work_size = max(test->maxK, test->maxN);
    char u, diag;
    double Anorm, Rnorm, result, Xnorm;
    double eps = LAPACKE_dlamch_work('e');
    double alpha;
    double *work = (double *)malloc(max_work_size*sizeof(double));

    /*Variable to save the residual alphaB -AX */
    double **residual;
    residual = (double **) malloc(batch_count*sizeof(double *));

    /*Make a copy of test->arrayB (X) in residual */
    bblas_dcopy_Binit(test, residual);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    N = test->N[batch_iter];
	    M = test->M[batch_iter];
	    ldb = test->ldb[batch_iter];

	    Am = test->lda[batch_iter];

	    if(test->side[batch_iter] == BblasLeft )
	    {
		An = M;
	    }else
	    {
		An = N;
	    }

	    if(test->uplo[batch_iter] == BblasLower )
	    {
	       u = 'L';
	    }else
	    {
	       u = 'U';
	    }

	    if(test->diag[batch_iter] == BblasUnit )
	    {
		diag = 'U';
	    }else
	    {
		diag = 'N';
	    }

	    /*Compute residual <- AX  (dtrmm) */
	    alpha = 1;
	    cblas_dtrmm(BblasColMajor, test->side[batch_iter], test->uplo[batch_iter], test->transA[batch_iter],
			test->diag[batch_iter], M, N, (alpha), test->arrayA[batch_iter],
			test->lda[batch_iter], residual[batch_iter], ldb);

	    /*Compute Binit = alpha Binit */
	    for (int index=0;  index < ldb*N; index++)
	    {
		test->Binit[batch_iter][index] = test->alpha[batch_iter]*test->Binit[batch_iter][index];
	    }

	    /*Compute the  residual <- alpha B - Ax  */
	    alpha= -1;
	    cblas_daxpy (ldb*N, (alpha), test->Binit[batch_iter], 1, residual[batch_iter], 1);
		/*Compute  the infinity norm associated  with the residual */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, residual[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_dlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

	    /*Compute  the infinity norm associated  with X (B) */
	    Xnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, test->arrayB[batch_iter], ldb, work);

	    result = Rnorm/(Anorm*Xnorm*N*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }

	}
    }else  if( test->batch_opts == BBLAS_FIXED )
    {

	N = test->N[first_index];
	M = test->M[first_index];
	ldb = test->ldb[first_index];

	Am = test->lda[first_index];

	if(test->side[first_index] == BblasLeft )
	{
	    An = M;
	}else
	{
	    An = N;
	}

	if(test->uplo[first_index] == BblasLower )
	{
	    u = 'L';
	}else
	{
	    u = 'U';
	}

	if(test->diag[first_index] == BblasUnit )
	{
	    diag = 'U';
	}else
	{
	    diag = 'N';
	}

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{

	    /*Compute residual <- AX  (dtrmm) */
	    alpha = 1;
	    cblas_dtrmm(BblasColMajor, test->side[first_index], test->uplo[first_index], test->transA[first_index],
			test->diag[first_index], M, N, (alpha), (const double*) test->arrayA[batch_iter],
			test->lda[first_index], residual[batch_iter], ldb);

	    /*Compute Binit = alpha Binit */
	    for (int index=0;  index < ldb*N; index++)
	    {
		test->Binit[batch_iter][index] = test->alpha[first_index]*test->Binit[batch_iter][index];
	    }

	    /*Compute the  residual <- alpha B - Ax  */
	    alpha = -1;
	    cblas_daxpy (ldb*N, (alpha), test->Binit[batch_iter], 1, residual[batch_iter], 1);

		/*Compute  the infinity norm associated  with the residual */
	    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, residual[batch_iter], ldb, work);

	    /*Infinity norm of the triangular matrix A */
	    Anorm = LAPACKE_dlantr_work(LAPACK_COL_MAJOR, 'I', u, diag, Am,
					An, test->arrayA[batch_iter], Am, work);

	    /*Compute  the infinity norm associated  with X (B) */
	    Xnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'I', ldb, N, test->arrayB[batch_iter], ldb, work);

	    result = Rnorm/(Anorm*Xnorm*N*eps);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = result;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = result;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = result;
		    break;
		default:
		    printf("In bblas_dcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }
    else
    {
	bblas_error("bblas_dtesting.c", "wrong batch_opts value");
    }

    /*Free allocated memory */
    for(int index =0; index < batch_count; index++)
    {
	free(residual[index]);
    }
    free(residual);
    free(work);
}

/**
 * Calls bblas_dstar2k_tolerance.
 **/
#ifdef COMPLEX
void bblas_dsyr2k_tolerance(bblas_dtest_t *test, double **C_final)
{
    bblas_dstar2k_tolerance(test, C_final);

}
#endif
/**
 * Calls bblas_dstar2k_tolerance.
 **/

void bblas_dsyr2k_tolerance(bblas_dtest_t *test, double **C_final)
{
    bblas_dstar2k_tolerance(test, C_final);
}

/**
 * Calls bblas_dstark_tolerance.
 **/

#ifdef COMPLEX
void bblas_dsyrk_tolerance(bblas_dtest_t *test, double **C_final)
{
    bblas_dstark_tolerance(test, C_final);
}
#endif
/**
 * Calls bblas_dstark_tolerance.
 **/

void bblas_dsyrk_tolerance(bblas_dtest_t *test, double **C_final)
{
    bblas_dstark_tolerance(test, C_final);
}

/**
 * Calls bblas_dstarmm_tolerance.
 **/

#ifdef COMPLEX
void bblas_dsymm_tolerance(bblas_dtest_t *test, double **C_final)
{
    bblas_dstarmm_tolerance(test, C_final);
}
#endif
/**
 * Calls bblas_dstamm_tolerance.
 **/

void bblas_dsymm_tolerance(bblas_dtest_t *test, double **C_final)
{
    bblas_dstarmm_tolerance(test, C_final);
}

ext: c

max_stars_repo_name: NLAFET/BBLAS-ref

max_issues_repo_name: NLAFET/BBLAS-ref

max_forks_repo_name: NLAFET/BBLAS-ref


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

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_combination.h>

size_t
gsl_combination_n (const gsl_combination * c)
{
  return c->n ;
}

size_t
gsl_combination_k (const gsl_combination * c)
{
  return c->k ;
}

size_t *
gsl_combination_data (const gsl_combination * c)
{
  return c->data ;
}

#ifndef HIDE_INLINE_STATIC
size_t
gsl_combination_get (const gsl_combination * c, const size_t i)
{
  if (gsl_check_range)
    {
      if (i >= c->k)            /* size_t is unsigned, can't be negative */
        {
          GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
        }
    }

  return c->data[i];
}
#endif

int
gsl_combination_valid (gsl_combination * c)
{
  const size_t n = c->n ;
  const size_t k = c->k ;

  size_t i, j ;

  if( k > n )
    {
      GSL_ERROR("combination has k greater than n", GSL_FAILURE) ;
    }
  for (i = 0; i < k; i++)
    {
      const size_t ci = c->data[i];

      if (ci >= n)
        {
          GSL_ERROR("combination index outside range", GSL_FAILURE) ;
        }

      for (j = 0; j < i; j++)
        {
          if (c->data[j] == ci)
            {
              GSL_ERROR("duplicate combination index", GSL_FAILURE) ;
            }
          if (c->data[j] > ci)
            {
              GSL_ERROR("combination indices not in increasing order",
                        GSL_FAILURE) ;
            }
        }
    }

  return GSL_SUCCESS;
}

int
gsl_combination_next (gsl_combination * c)
{
  /* Replaces c with the next combination (in the standard lexicographical
   * ordering).  Returns GSL_FAILURE if there is no next combination.
   */
  const size_t n = c->n;
  const size_t k = c->k;
  size_t *data = c->data;
  size_t i;

  if(k == 0)
    {
      return GSL_FAILURE;
    }
  i = k - 1;

  while(i > 0 && data[i] == n - k + i)
    {
      i--;
    }
  if(i == 0 && data[i] == n - k)
    {
      return GSL_FAILURE;
    }
  data[i]++;
  for(; i < k - 1; i++)
    {
      data[i + 1] = data[i] + 1;
    }
  return GSL_SUCCESS;
}

int
gsl_combination_prev (gsl_combination * c)
{
  /* Replaces c with the previous combination (in the standard
   * lexicographical ordering).  Returns GSL_FAILURE if there is no
   * previous combination.
   */
  const size_t n = c->n;
  const size_t k = c->k;
  size_t *data = c->data;
  size_t i;

  if(k == 0)
    {
      return GSL_FAILURE;
    }
  i = k - 1;

  while(i > 0 && data[i] == data[i-1] + 1)
    {
      i--;
    }
  if(i == 0 && data[i] == 0)
    {
      return GSL_FAILURE;
    }
  data[i++]--;
  for(; i < k; i++)
    {
      data[i] = n - k + i;
    }
  return GSL_SUCCESS;
}

int
gsl_combination_memcpy (gsl_combination * dest, const gsl_combination * src)
{
   const size_t src_n = src->n;
   const size_t src_k = src->k;
   const size_t dest_n = dest->n;
   const size_t dest_k = dest->k;

   if (src_n != dest_n || src_k != dest_k)
     {
       GSL_ERROR ("combination lengths are not equal", GSL_EBADLEN);
     }

   {
     size_t j;

     for (j = 0; j < src_k; j++)
       {
         dest->data[j] = src->data[j];
       }
   }

   return GSL_SUCCESS;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
 *   1dcol.c
 *
 *   MPI code for performing parallel dense matrix-vector multiplication using 1D
 *   columnwise partitioning.
 *
 *   Written by cetinsamet -*- cetin.samet@metu.edu.tr
 *   April, 2019
 *
 */
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>

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

	int i, j;
    int rank, size;

    double minStart, maxEnd;

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int N 		= atoi(argv[1]);
    int pSize	= N / size;

    /* Initialize result vector */
    double* RESULT;
    RESULT 	= (double*) malloc(N * sizeof(double));

	/* Initialize local result vector */
	double* locResult;
	locResult = (double*) malloc(N * sizeof(double));

    /* Initialize matrix */
    double **matrix	= (double **) malloc(N * sizeof(double*));
	for(i = 0; i < N; i++) {
		matrix[i] = (double *) malloc(N * sizeof(double));
	}
	for (i = 0; i < N; ++i) {
		for (j = 0; j < N; ++j)
			matrix[i][j] = i + j;
	}

	/* Initialize vector */
	double* vector;
	vector = (double*) malloc(N * sizeof(double));
	for (i = 0; i < N; ++i)
		vector[i] = i;

	/* Initialize matrix */
	double **p_matrix	= (double **) malloc(N * sizeof(double*));
	for(i = 0; i < N; i++)
		p_matrix[i] = (double *) malloc(pSize * sizeof(double));

	for (j = 0; j < pSize; ++j) {
		for (i = 0; i < N; ++i)
			p_matrix[i][j] = matrix[i][(rank * pSize) + j];
	}

	/* Initialize local vector */
	double* p_vector;
	p_vector = (double*) malloc(N * sizeof(double));
	for (i = 0; i < pSize; ++i) {
		p_vector[i] = vector[(rank * pSize) + i];
	}

	double pStart = MPI_Wtime(); // <-- Start Time

	/* Calculate local dot product */
	for (i = 0; i < pSize; ++i) {
		for (j = 0; j < N; ++j) {
			locResult[j] += p_matrix[j][i] * p_vector[i];

	for (i = 0; i < N; ++i)
		MPI_Reduce(&locResult[i], &RESULT[i], 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

	double pEnd = MPI_Wtime(); // <-- Stop Time

	MPI_Reduce(&pStart,  &minStart, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
    MPI_Reduce(&pEnd, &maxEnd, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

    /* Print result */
	if (rank == 0) {
		printf("RESULT of MATRIX-VECTOR MULTIPLICATION:\n");
		for (i = 0; i < N; ++i)
			printf("%f ", RESULT[i]);
		printf("\n");
		printf("Vector size: %d\tElapsed time: %f\n", N, maxEnd - minStart);
	}

	MPI_Finalize();
	return 0;
}


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit.h>

int
main()
{
  const size_t n = 1000; /* number of observations */
  const size_t p = 2;    /* number of model parameters */
  size_t i;
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  gsl_matrix *X = gsl_matrix_alloc(n, p);
  gsl_vector *y = gsl_vector_alloc(n);

  for (i = 0; i < n; ++i)
    {
      /* generate first random variable u */
      double ui = 5.0 * gsl_ran_gaussian(r, 1.0);

      /* set v = u + noise */
      double vi = ui + gsl_ran_gaussian(r, 0.001);

      /* set y = u + v + noise */
      double yi = ui + vi + gsl_ran_gaussian(r, 1.0);

      /* since u =~ v, the matrix X is ill-conditioned */
      gsl_matrix_set(X, i, 0, ui);
      gsl_matrix_set(X, i, 1, vi);

      /* rhs vector */
      gsl_vector_set(y, i, yi);
    }

  {
    const size_t npoints = 200;                   /* number of points on L-curve and GCV curve */
    gsl_multifit_linear_workspace *w =
      gsl_multifit_linear_alloc(n, p);
    gsl_vector *c = gsl_vector_alloc(p);          /* OLS solution */
    gsl_vector *c_lcurve = gsl_vector_alloc(p);   /* regularized solution (L-curve) */
    gsl_vector *c_gcv = gsl_vector_alloc(p);      /* regularized solution (GCV) */
    gsl_vector *reg_param = gsl_vector_alloc(npoints);
    gsl_vector *rho = gsl_vector_alloc(npoints);  /* residual norms */
    gsl_vector *eta = gsl_vector_alloc(npoints);  /* solution norms */
    gsl_vector *G = gsl_vector_alloc(npoints);    /* GCV function values */
    double lambda_l;                              /* optimal regularization parameter (L-curve) */
    double lambda_gcv;                            /* optimal regularization parameter (GCV) */
    double G_gcv;                                 /* G(lambda_gcv) */
    size_t reg_idx;                               /* index of optimal lambda */
    double rcond;                                 /* reciprocal condition number of X */
    double chisq, rnorm, snorm;

    /* compute SVD of X */
    gsl_multifit_linear_svd(X, w);

    rcond = gsl_multifit_linear_rcond(w);
    fprintf(stderr, "matrix condition number = %e\n\n", 1.0 / rcond);

    /* unregularized (standard) least squares fit, lambda = 0 */
    gsl_multifit_linear_solve(0.0, X, y, c, &rnorm, &snorm, w);
    chisq = pow(rnorm, 2.0);

    fprintf(stderr, "=== Unregularized fit ===\n");
    fprintf(stderr, "best fit: y = %g u + %g v\n",
      gsl_vector_get(c, 0), gsl_vector_get(c, 1));
    fprintf(stderr, "residual norm = %g\n", rnorm);
    fprintf(stderr, "solution norm = %g\n", snorm);
    fprintf(stderr, "chisq/dof = %g\n", chisq / (n - p));

    /* calculate L-curve and find its corner */
    gsl_multifit_linear_lcurve(y, reg_param, rho, eta, w);
    gsl_multifit_linear_lcorner(rho, eta, &reg_idx);

    /* store optimal regularization parameter */
    lambda_l = gsl_vector_get(reg_param, reg_idx);

    /* regularize with lambda_l */
    gsl_multifit_linear_solve(lambda_l, X, y, c_lcurve, &rnorm, &snorm, w);
    chisq = pow(rnorm, 2.0) + pow(lambda_l * snorm, 2.0);

    fprintf(stderr, "\n=== Regularized fit (L-curve) ===\n");
    fprintf(stderr, "optimal lambda: %g\n", lambda_l);
    fprintf(stderr, "best fit: y = %g u + %g v\n",
            gsl_vector_get(c_lcurve, 0), gsl_vector_get(c_lcurve, 1));
    fprintf(stderr, "residual norm = %g\n", rnorm);
    fprintf(stderr, "solution norm = %g\n", snorm);
    fprintf(stderr, "chisq/dof = %g\n", chisq / (n - p));

    /* calculate GCV curve and find its minimum */
    gsl_multifit_linear_gcv(y, reg_param, G, &lambda_gcv, &G_gcv, w);

    /* regularize with lambda_gcv */
    gsl_multifit_linear_solve(lambda_gcv, X, y, c_gcv, &rnorm, &snorm, w);
    chisq = pow(rnorm, 2.0) + pow(lambda_gcv * snorm, 2.0);

    fprintf(stderr, "\n=== Regularized fit (GCV) ===\n");
    fprintf(stderr, "optimal lambda: %g\n", lambda_gcv);
    fprintf(stderr, "best fit: y = %g u + %g v\n",
            gsl_vector_get(c_gcv, 0), gsl_vector_get(c_gcv, 1));
    fprintf(stderr, "residual norm = %g\n", rnorm);
    fprintf(stderr, "solution norm = %g\n", snorm);
    fprintf(stderr, "chisq/dof = %g\n", chisq / (n - p));

    /* output L-curve and GCV curve */
    for (i = 0; i < npoints; ++i)
      {
        printf("%e %e %e %e\n",
               gsl_vector_get(reg_param, i),
               gsl_vector_get(rho, i),
               gsl_vector_get(eta, i),
               gsl_vector_get(G, i));
      }

    /* output L-curve corner point */
    printf("\n\n%f %f\n",
           gsl_vector_get(rho, reg_idx),
           gsl_vector_get(eta, reg_idx));

    /* output GCV curve corner minimum */
    printf("\n\n%e %e\n",
           lambda_gcv,
           G_gcv);

    gsl_multifit_linear_free(w);
    gsl_vector_free(c);
    gsl_vector_free(c_lcurve);
    gsl_vector_free(reg_param);
    gsl_vector_free(rho);
    gsl_vector_free(eta);
    gsl_vector_free(G);
  }

  gsl_rng_free(r);
  gsl_matrix_free(X);
  gsl_vector_free(y);

  return 0;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE