#include <stdio.h>
#include <stdarg.h>
// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears
#define MPFR_USE_VA_LIST
#include <mpfr.h>
#include <libinum.h>

int func(mpfr_t* out, mpfr_t in, void* extra_args);

int main(int argc, const char* argv[]){
  int ret;
  mpfr_t tolerance, val, lower, upper;
  array_mpfr abcissa, weights;

  // precision of MPFR floats
  mpfr_set_default_prec(53);

  // tolerance for computing weights
  mpfr_init(tolerance);
  mpfr_set_d(tolerance, 1.e-11, MPFR_RNDN);

  // compute weights
  ret=gauss_legendre_weights_mpfr(10, tolerance, 10000, &abcissa, &weights);
  // return codes
  if(ret==LIBINUM_ERROR_MAXITER){
    printf("error: maximum number of iterations reached when computing the integration abcissa\n");
    mpfr_clear(tolerance);
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    printf("error: infinity encountered when computing the integration abcissa\n");
    mpfr_clear(tolerance);
    return(ret);
  }
  else{
    printf("abcissa:\n");
    array_mpfr_print(abcissa);
    printf("\nweights:\n");
    array_mpfr_print(weights);
    printf("\n");
  }

  mpfr_clear(tolerance);

  mpfr_inits(val, lower, upper, NULL);

  mpfr_set_ui(lower, 0, MPFR_RNDN);
  mpfr_set_ui(upper, 2, MPFR_RNDN);

  ret=integrate_gauss_mpfr(&val, &func, lower, upper, abcissa, weights, NULL);
  // return codes
  if(ret==LIBINUM_ERROR_SIZE_MISMATCH){
    printf("error: the mpfr_arrays 'abcissa' and 'weights' must have the same length\n");
    mpfr_clears(val, upper, lower, NULL);
    array_mpfr_free(abcissa);
    array_mpfr_free(weights);
    return(ret);
  }
  else if(ret==LIBINUM_ERROR_NAN){
    printf("error: the integrand is singular\n");
    array_mpfr_free(abcissa);
    array_mpfr_free(weights);
    mpfr_clears(val, upper, lower, NULL);
    return(ret);
  }
  else{
    printf("\\int_0^2 x/sin(x) dx = ");
    fprint_mpfr(stdout, val);
    printf("\n");
  }

  mpfr_clears(val, upper, lower, NULL);
  array_mpfr_free(abcissa);
  array_mpfr_free(weights);

  return(0);
}

/*
// e^x
int func(mpfr_t* out, mpfr_t in, void* extra_args){
  mpfr_exp(*out, in, MPFR_RNDN);
  return(0);
}
*/

// x/sin(x)
int func(mpfr_t* out, mpfr_t in, void* extra_args){
  mpfr_sin(*out, in, MPFR_RNDN);
  mpfr_div(*out, in, *out, MPFR_RNDN);
  return(0);
}