/*
* Polynomial Interpolation *
* This program demonstrates a function that performs polynomial
* interpolation. The function is taken from "Numerical Recipes
* in C", Second Edition, by William H. Press, Saul A. Teukolsky,
* William T. Vetterling, and Brian P. Flannery.
*
*/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <omp.h>

#define N 20
#define X 14.5

#define NUM_THREADS 7

/* Number of function sample points */
/* Interpolate at this value of x */
/* Function 'vector' is used to allocate vectors with subscript
range v[nl..nh] */
double *vector (long nl, long nh)
{
  double *v;
  v = (double *) malloc(((nh-nl+2)*sizeof(double)));
  return v-nl+1;
}
/* Function 'free_vector' is used to free up memory allocated
with function 'vector' */
void free_vector(double *v, long nl, long nh)
{
  free ((char *) (v+nl-1));
}
/* Function 'polint' performs a polynomial interpolation */
void polint (double xa[], double ya[], int n, double x, double *y, double *dy) {
  int i, m, ns=1;
  double den,dif,dift,ho,hp,w;
  double *c, *d;
  dif = fabs(x-xa[1]);
  c = vector(1,n);
  d = vector(1,n);
  #pragma omp parallel

  #pragma omp parallel for private (i)
  for (i=1; i<= n; i++) {
    #pragma omp critical

    dift = fabs (x - xa[i]);
    if (dift<dif) {
    ns = i;
    dif = dift;
    }
    c[i] = ya[i];
    d[i] = ya[i];
  }
  *y = ya[ns--];

    #pragma omp parallel
    #pragma omp parallel for private (i)

  for (m = 1; m < n; m++) {
  for (i = 1; i<= n-m; i++) {
       #pragma omp critical

  ho = xa[i] - x;
  hp = xa[i+m] - x;
  w = c[i+1] - d[i]; den = ho - hp; den = w / den; d[i] = hp * den; c[i] = ho * den; }
  *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--])); }
  free_vector (d, 1, n);
  free_vector (c, 1, n);
}

/* Functions 'sign' and 'init' are used to initialize the
x and y vectors holding known values of the function.
*/
int sign (int j)
{
  if (j % 2 == 0) return 1;
  else return -1;
}

void init (int i, double *x, double *y)
{
  // int j;
  *x = (double) i;
  *y = sin(i);
}
/* Function 'main' demonstrates the polynomial interpolation function
by generating some test points and then calling 'polint' with a
value of x between two of the test points. */
int main (int argc, char *argv[])
{
  double x, y, dy;
  double *xa, *ya;
  int i;

  double start_time = omp_get_wtime();
  omp_set_num_threads(NUM_THREADS); //define el numero de hilos

  xa = vector (1, N);
  ya = vector (1, N);
  /* Initialize xa's and ya's */

  #pragma omp parallel
  {
    #pragma omp for
    for (i = 1; i<= N; i++) {
     #pragma omp critical

      init (i, &xa[i], &ya[i]);
      printf ("f(%4.2f) = %13.11f\n", xa[i], ya[i]);
    }
  }

  /* Interpolate polynomial at X */
  polint (xa, ya, N, X, &y, &dy);
  printf ("\nf(%6.3f) = %13.11f with error bound %13.11f\n", X, y, fabs(dy)); free_vector (xa, 1, N);
  free_vector (ya, 1, N);
  printf("Time: \t %f \n", omp_get_wtime()-start_time);
  return 0;
}