#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_mathieu.h>

int main(void) {

  gsl_sf_result result;

  int n = 100;
  int n1 = 10000;
  int i;

  double q = 2000.0;

  double diff;

  gsl_sf_mathieu_workspace * workspace =  gsl_sf_mathieu_alloc(n1, q);

  gsl_sf_result old, new;

  gsl_sf_mathieu_a (0, q, &old);

  for (i = 1; i < n; ++i) {
    gsl_sf_mathieu_b (i, q, &new);

    diff = new.val - old.val;
    if (diff <= 0) {
      printf("b[%d] - a[%d-1] , %f\n", i, i, diff);
    }
    old = new;

    gsl_sf_mathieu_a (i, q, &new);

    diff = new.val - old.val;
    if (diff <= 0) {
      printf("a[%d] - b[%d], %f\n", i, i, diff);
    }
    old = new;

  }
  gsl_sf_mathieu_free(workspace);

  return 0;
}