/*
 * Command line numerical polynomial equation solver using the GNU Scientific Library (GSL).
 * GSL is available from: "http://www.gnu.org/software/gsl/".
 *
 * To compile:
 *
 * cc -O -Wall -o psolve psolve.c -lgsl -lgslcblas
 */
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_poly.h>

int
main(int argc, char **argv)
{
  int i, highest_power;
  char *cp;
  double *a, *z;
  gsl_poly_complex_workspace *w;

  if (argc <= 1) {
    fprintf(stderr, "Solve polynomial = 0 by giving all real coefficients.\n");
    fprintf(stderr, "Usage: %s highest-power-coefficient ... constant-term\n", argv[0]);
    exit(1);
  }

  highest_power = argc - 2;
  a = calloc(highest_power + 1, sizeof(double)); /* allocate real double input array */
  z = calloc(2 * highest_power, sizeof(double)); /* allocate complex double output array */

  for (i = 0; i < argc - 1; i++) {
    a[i] = strtod(argv[argc-i-1], &cp);
    if (argv[argc-i-1] == cp) {
      fprintf(stderr, "Invalid floating point number.\n");
      exit(1);
    }
  }

  printf("The %d floating point solutions of:\n", highest_power);
  for (i = highest_power; i >= 0; i--) {
    if (a[i]) {
      printf("%+.14g", a[i]);
      if (i) {
        printf("*x");
        if (i > 1) {
          printf("^%d", i);
        }
      }
      printf(" ");
    }
  }
  printf("= 0\nare:\n");

  w = gsl_poly_complex_workspace_alloc(highest_power + 1);

  gsl_poly_complex_solve(a, highest_power + 1, w, z);

  gsl_poly_complex_workspace_free(w);

  for (i = 0; i < highest_power; i++) {
    printf("x = %+.14g", z[2*i]);
    if (z[2*i+1])
      printf(" %+.14g*i", z[2*i+1]);
    printf("\n");
  }

  exit(0);
}