#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(n, q);

  gsl_sf_result old, new;

  double diff_array;
  double *old_array = malloc(n*sizeof(double));
  double *new_array = malloc(n*sizeof(double));
  gsl_sf_mathieu_a_array(0, n, q, workspace, old_array);
  gsl_sf_mathieu_b_array(1, n, q, workspace, new_array);

  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;
    diff_array = new_array[i-1] - old_array[i-1];
    if (diff <= 0) {
        printf("b[%d] - a[%d] , %f, %f\n", i, i-1, diff, diff_array);
    }
    old = new;

    gsl_sf_mathieu_a (i, q, &new);

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

  return 0;
}