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

/*************************************************
 * Evaluate a Newton polynomial at a number of
 * points using Horner's rule
 *
 * Compile with
 *     gcc -o blas -Wall -I/usr/include blas.c -lgsl -lgslcblas -lm
 *     ./blas
 *************************************************/

#define N 5

int main() {
        /* The points
             (0.1, 1.2), (0.5, 2.7), (0.7, 3.4), (1.2, 4.7), (1.5, 6.0) */

        double xs[N] = {0.1, 0.5, 0.7, 1.2, 1.5};
        double ys[N] = {1.2, 2.7, 3.8, 4.7, 6.0};

        /* The coefficients of the newton polynomial. */
        double dd[N];
        gsl_poly_dd_init( dd, xs, ys, N );

        double x;

        for ( x = -2.0; x <= 2.0; x += 0.5 ) {
                printf( "(%+2.1lf, %12.8lf)\n",
                        x, gsl_poly_dd_eval( dd, xs, N, x ) );
        }

        return 0;
}