GSL: creating complex numbers

In a GNU Scientific Library C program I am testing for internal consistency of roots of polynomials. In the sample program I am cooking up a random polynomial, find its roots and evaluate the polynomial at those roots. In other words, I am expecting zeroes in the last block of output:

/**
 * Compilation: gcc polynomialevaluation.c -Wall -g `gsl-config --cflags --libs` -o polynomialevaluation
 * Usage: ./polynomialevaluation
 *
 */

#include <stdio.h>
#include <time.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_poly.h>

#define CARDINALITY 3

int i = 0;                     // index of coefficient
double a[CARDINALITY];         // coefficients of polynomial
double z[(CARDINALITY-1)*2];   // roots of polynomial, twice as large as cardinality
gsl_complex zp;

int main (void) {
  srand(time(0));
  // generate coefficients and display them
  printf("Polynomial with random integer coefficients:\n");
  for (i = 0; i < CARDINALITY+1; i++) { a[i] = (rand()%7)-2; printf ("%d: %+.0f\n", i, a[i]); }
  if (a[CARDINALITY] == 0) a[CARDINALITY] = 1;
  printf("\n");
  // a workspace contains parameters used for iterative rootfinding of general polynomials
  gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (CARDINALITY+1);
  gsl_poly_complex_solve (a, CARDINALITY+1, w, z);
  gsl_poly_complex_workspace_free (w);

  printf("Roots:\n");
  for (i = 0; i < CARDINALITY; i++) printf ("z[%d] = %+.5f %+.5f\n", i, z[2*i], z[2*i + 1]);
  printf("\n");
  printf("Evaluation at roots:\n");
  // Attempt 0: show real part only
  for (i = 0; i < CARDINALITY; i++) printf ("z%d = %+.5f\n", i, gsl_poly_eval(a, CARDINALITY+1, z[2*i]));
  // Attempt 1: use gsl_complex_rect
  //for (i = 0; i < CARDINALITY; i++) printf ("z%d = %+.5f\n", i, gsl_poly_eval(a, CARDINALITY+1, gsl_complex_rect(z[2*i], z[2*i + 1])));
  // Attempt 2: use GSL_SET_COMPLEX
  //for (i = 0; i < CARDINALITY; i++) printf ("z%d = %+.5f\n", i, gsl_poly_eval(a, CARDINALITY+1, GSL_SET_COMPLEX(&zp, z[2*i], z[2*i + 1])));
  printf("\n");

  return 0;
}

Attempt 0 is obviously not good enough, as it only evaluates the polynomial at the real part. I am looking for a way to forge complex numbers out of the z-array, ie the array of roots. In the commented out parts two approaches were taken from the documentation, but leading to compilation errors:
Attempt 1: warning: implicit declaration of function ‘gsl_complex_rect’ [-Wimplicit-function-declaration]
Attempt 2: error: expected expression before ‘do’

Also tried with several combinations: with or without ampersand (&), with or without asterisk (*), with or without header file (gsl/gsl_complex_math.h). Can someone tell me how I can create complex numbers out of z to use the gsl_poly_eval function?

Sample output:

$ ./polynomialevaluation 
Polynomial with random integer coefficients:
0: -2
1: +4
2: -1
3: +3

Roots:
z[0] = +0.47581 +0.00000
z[1] = -0.07124 +1.18155
z[2] = -0.07124 -1.18155

Evaluation at roots:
z0 = -0.00000
z1 = -2.29111
z2 = -2.29111

Attempt 1: you may need an additional #include <gsl_complex_math.h>

From https://www.gnu.org/software/gsl/doc/html/complex.html:

The complex types are defined in the header file gsl_complex.h, while the corresponding complex functions and arithmetic operations are defined in gsl_complex_math.h.

Attempt 2: This looks like a bad macro definition (or possible misuse) in one of the include files: do ... while (0) is a common construct for Xwin and complex number code, because it acts as a block wrapper and fixes a difficult syntax issue with C. (Google do while(0) macro in c). The compiler should have given you a filename and line number for this, so you should be able to find it.

Hi @technossomy,

there are some inaccuracies:

• #include <gsl/gsl_complex_math.h> // for gsl_complex_rect(), if used
• double a[CARDINALITY+1] // N+1 coeffs, not N
• double z[CARDINALITY*2] // N roots, so 2*N real/imag values
• GSL_SET_COMPLEX is a macro that returns nothing, but sets values

So this should work:

gsl_complex zp;
for (i = 0; i < CARDINALITY; i++) {
    GSL_SET_COMPLEX(&zp, z[2*i], z[2*i + 1]);
    printf ("z%d = %+.5f\n", i, gsl_poly_complex_eval(a, CARDINALITY+1, zp));
}

gsl_complex_rect() returns a gsl_complex type, so gsl_poly_eval(..., gsl_complex_rect()) won't work.

I have settled for the following, given that gsl_poly_complex_eval does not return a float:

gsl_complex zp;
for (i = 0; i < CARDINALITY; i++) {
    GSL_SET_COMPLEX(&zp, z[2*i], z[2*i + 1]);
    printf ("z%d = %+.5f %+.5f\n", i, GSL_REAL(gsl_poly_complex_eval(a, CARDINALITY+1, zp)), GSL_IMAG(gsl_poly_complex_eval(a, CARDINALITY+1, zp)));
  }

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.