/*
 * gslbug.c
 * compile:	gcc -O2 -Wall -Wextra -o gslbug gslbug.c -lgsl -lgslcblas
 *
 * A small demo for some bugs in the GSL (V1.8)
 * "Trigonometric Restriction" routines
 */

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_errno.h>

void
my_handler(const char *reason,
           const char *file,
           int line,
           int gsl_errno) {
  printf("GSL error %d (%s:%d %s) ignored\n",
    gsl_errno, file, line, reason);
}

int
main(void) {
  int i;
  double val, res;
  gsl_sf_result result;
  gsl_error_handler_t *old_handler;
  const double TwoPi = 2*M_PI;

  printf("\n\tDemonstrate bugs in the GSL(-1.8) routines\n"
    "\t\tgsl_sf_angle_restrict_symm and gsl_sf_angle_restrict_pos\n");
  old_handler = gsl_set_error_handler(my_handler);

/*
 * Part I.A. Missing check for huge negative arguments
 */
  printf("\n\tPart I.A. Missing check for huge negative arguments:\n");
  val = 1.0e+14;
  for (i=0; i<2; i++) {
    printf("gsl_sf_angle_restrict_symm(%g) => ", val);
    printf("%g\n", gsl_sf_angle_restrict_symm(val));
    printf("gsl_sf_angle_restrict_pos(%g) => ", val);
    printf("%g\n", gsl_sf_angle_restrict_pos(val));
    printf("gsl_sf_angle_restrict_symm(%g) => ", -val);
    printf("%g\n", gsl_sf_angle_restrict_symm(-val));
    printf("gsl_sf_angle_restrict_pos(%g) => ", -val);
    printf("%g\n", gsl_sf_angle_restrict_pos(-val));
    val *= 10.0;
  }

/*
 * Part I.B. gsl_sf_angle_restrict_pos(x)<0
 */
  printf("\n\tPart I.B. gsl_sf_angle_restrict_pos(x)<0:\n");
  for (i=-3; i<6; i+=3)
    printf("gsl_sf_angle_restrict_pos(TwoPi%+d*DBL_EPSILON) = %.17g\n",
      i, gsl_sf_angle_restrict_pos(TwoPi + i*DBL_EPSILON));

/*
 * Part II.1. Wrong error estimates
 */
  printf("\n\tPart II.1. Wrong error estimates:\n");
  val = 1.0e-13;
  for (i=0; i<3; i++) {
    gsl_sf_angle_restrict_symm_err_e(val, &result);
    printf("gsl_sf_angle_restrict_symm_err_e(%.15e) => %.15e +- %g\n",
           val, result.val, result.err);
    gsl_sf_angle_restrict_symm_err_e(-val, &result);
    printf("gsl_sf_angle_restrict_symm_err_e(%.15e) => %.15e +- %g\n",
           -val, result.val, result.err);
    gsl_sf_angle_restrict_pos_err_e(val, &result);
    printf("gsl_sf_angle_restrict_pos_err_e(%.15e) => %.15e +- %g\n",
           val, result.val, result.err);
    gsl_sf_angle_restrict_pos_err_e(-val, &result);
    printf("gsl_sf_angle_restrict_pos_err_e(%.15e) => %.15e +- %g\n",
           -val, result.val, result.err);
    val /= 10.0;
  }
  printf("\n\t           More error estimates:\n");
  for (i=-3; i<9; i+=3) {
    val = M_PI + i*DBL_EPSILON;
    gsl_sf_angle_restrict_symm_err_e(val, &result);
    printf("gsl_sf_angle_restrict_symm_err_e(Pi%+d*DBL_EPSILON) => %.16f +- %g\n",
           i, result.val, result.err);
  }
  printf("\n");
  for (i=-3; i<9; i+=3) {
    val = TwoPi + i*DBL_EPSILON;
    gsl_sf_angle_restrict_pos_err_e(val, &result);
    printf("gsl_sf_angle_restrict_pos_err_e(TwoPi%+d*DBL_EPSILON) => %.16f +- %g\n",
           i, result.val, result.err);
  }

/*
 * Part II.2. Huge (positive or negative) results
 */
  printf("\n\tPart II.2. Huge (positive or negative) results:\n");
  val = +1.0e+300;
  do {
    res = gsl_sf_angle_restrict_pos(val);
    printf("gsl_sf_angle_restrict_pos(%g) = %g\n", val, res);
    val = res;
  } while (fabs(val) > 1.0e+200);

  printf("\n");
  return 0;
}