C mpsolve library example

#include <gmp.h>
#include <stdio.h>
#include <mps/mps.h>


void
solve_poly(long int degree, mpq_t *coeffs_re, mpq_t *coeffs_im,
           long output_prec)
{
  mpc_t *roots = NULL;
  rdpe_t *rds = NULL;
  static mps_context *ctx;
  mps_monomial_poly *mn_pl;

  ctx = mps_context_new();
  mn_pl = mps_monomial_poly_new (ctx, degree);
  for (int m = 0; m < degree + 1; m++) {
    mps_monomial_poly_set_coefficient_q(ctx, mn_pl, m, coeffs_re[m], coeffs_im[m]);
  }
  mps_context_set_input_poly(ctx, MPS_POLYNOMIAL(mn_pl));
  mps_context_select_algorithm (ctx, MPS_ALGORITHM_SECULAR_GA);
  mps_context_set_output_goal(ctx, MPS_OUTPUT_GOAL_APPROXIMATE);
  mps_context_set_output_prec(ctx, output_prec);
  printf("Enter solver\n");
  mps_mpsolve(ctx);
  printf("Leave solver\n");
  mps_context_get_roots_m(ctx, &roots, &rds);
  for (int m = 0; m < degree; m++) {
    fprintf(stdout, "%.17e + I*%.17e\n", mpf_get_d(mpc_Re(roots[m])),
            mpf_get_d(mpc_Im(roots[m])));
  }
  mpc_vclear(roots, degree);
  free(roots);
  free(rds);
  mps_monomial_poly_free(ctx, MPS_POLYNOMIAL(mn_pl));
  mps_context_free(ctx);
}

void
test_mps(void)
{
  int M = 16;
  long output_prec = 100;
  mpq_t *coeffs_re = NULL, *coeffs_im = NULL;

  coeffs_re = malloc(sizeof(mpq_t)*(M + 1));
  coeffs_im = malloc(sizeof(mpq_t)*(M + 1));
  for (int m = 0; m < M+1; m++) {
    mpq_init(coeffs_re[m]);
    mpq_init(coeffs_im[m]);
  }
  mpq_set_si(coeffs_re[0], 1, 1);
  mpq_set_si(coeffs_re[M], -1, 1);
  solve_poly(M, coeffs_re, coeffs_im, output_prec);
  mpq_vclear(coeffs_re, (unsigned)M + 1);
  mpq_vclear(coeffs_im, (unsigned)M + 1);
  free(coeffs_re);
  free(coeffs_im);
}