#include <gsl/gsl_sf_erf.h>

#include "cctk.h"



/* Complementary Error Function
 * erfc(x) := 2/Sqrt[Pi] Integrate[Exp[-t^2], {t,x,Infinity}]
 *
 * exceptions: none
 */
int CCTK_FCALL
CCTK_FNAME(gsl_sf_erfc_e) (double const * restrict const x,
                           gsl_sf_result * restrict const result)
{
  return gsl_sf_erfc_e (* x, result);
}

double CCTK_FCALL
CCTK_FNAME(gsl_sf_erfc) (double const * restrict const x)
{
  return gsl_sf_erfc (* x);
}


/* Log Complementary Error Function
 *
 * exceptions: none
 */
int CCTK_FCALL
CCTK_FNAME(gsl_sf_log_erfc_e) (double const * restrict const x,
                               gsl_sf_result * restrict const result)
{
  return gsl_sf_log_erfc_e (* x, result);
}

double CCTK_FCALL
CCTK_FNAME(gsl_sf_log_erfc) (double const * restrict const x)
{
  return gsl_sf_log_erfc (* x);
}


/* Error Function
 * erf(x) := 2/Sqrt[Pi] Integrate[Exp[-t^2], {t,0,x}]
 *
 * exceptions: none
 */
int CCTK_FCALL
CCTK_FNAME(gsl_sf_erf_e) (double const * restrict const x,
                          gsl_sf_result * restrict const result)
{
  return gsl_sf_erf_e (* x, result);
}

double CCTK_FCALL
CCTK_FNAME(gsl_sf_erf) (double const * restrict const x)
{
  return gsl_sf_erf (* x);
}


/* Probability functions:
 * Z(x) :  Abramowitz+Stegun 26.2.1
 * Q(x) :  Abramowitz+Stegun 26.2.3
 *
 * exceptions: none
 */
int CCTK_FCALL
CCTK_FNAME(gsl_sf_erf_Z_e) (double const * restrict const x,
                            gsl_sf_result * restrict const result)
{
  return gsl_sf_erf_Z_e (* x, result);
}

int CCTK_FCALL
CCTK_FNAME(gsl_sf_erf_Q_e) (double const * restrict const x,
                            gsl_sf_result * restrict const result)
{
  return gsl_sf_erf_Q_e (* x, result);
}

double CCTK_FCALL
CCTK_FNAME(gsl_sf_erf_Z) (double const * restrict const x)
{
  return gsl_sf_erf_Z (* x);
}

double CCTK_FCALL
CCTK_FNAME(gsl_sf_erf_Q) (double const * restrict const x)
{
  return gsl_sf_erf_Q (* x);
}


#if 0
/* Does not exist in older versions of GSL*/

/* Hazard function, also known as the inverse Mill's ratio.
 *
 *   H(x) := Z(x)/Q(x)
 *         = Sqrt[2/Pi] Exp[-x^2 / 2] / Erfc[x/Sqrt[2]]
 *
 * exceptions: GSL_EUNDRFLW
 */
int CCTK_FCALL
CCTK_FNAME(gsl_sf_hazard_e) (double const * restrict const x,
                             gsl_sf_result * restrict const result)
{
  return gsl_sf_hazard_e (* x, result);
}

double CCTK_FCALL
CCTK_FNAME(gsl_sf_hazard) (double const * restrict const x)
{
  return gsl_sf_hazard (* x);
}
#endif