#include <stdlib.h> //< exit
#include <assert.h> //< assert
#include <stdio.h>  //< printf
#include <math.h>   //< cos, fabs, M_PI

// COPIEZ-COLLEZ LA FONCTION SUIVANTE DANS VOTRE CODE POUR L'UTILISER.

/**
 * \brief Construit une formule de quadrature de Gauss-Legendre à n points sur l'intervalle [-1, 1].
 *
 * \param[out] points Stocke les n points de quadrature.
 * \param[out] weights Stocke les n poids de quadrature.
 * \param[in] n Nombre de points (doit être un entier pair strictement positif)
 *
 * En entrée, les tableaux points et weights doivent être déjà alloués à la taille n.
 * En sortie :
 *  - les indices [0 : (n-1)/2] contiennent les points < 0
 *  - les indices [(n-1)/2+1 : (n-1)] contiennent les points > 0
 *
 * Les points de quadrature sont les racines du polynôme de Legendre de degré n, calculées avec une 
 * méthode de Newton. 
 * Les poids de quadratures sont calculés avec la formule  w_j = 2 / ((1 - x_j^2) * [P_n'(x_j)^2]).
 * Pour plus de détails, voir "Numerical Recipes, 3rd ed." par William H. Press (p. 179-184).
 *
 */
void create_gauss_legendre_quadrature(double* points, double* weights, int n) {
  assert(points != NULL);
  assert(weights != NULL);

  if (n <= 0 || (n % 2) != 0) {
    printf("ERROR : create_gauss_legendre_quadrature : n doit être > 0 et pair.\n");
    exit(EXIT_FAILURE);
  }

  // On calcule seulement les racines positives car la formule est symétrique
  int const m = (n + 1) / 2;

  for (int i = 0; i < m; ++i) {

    // Valeur initiale pour la methode de Newton
    double z = cos(M_PI * (i + 0.75) / (n + 0.5));
    double z1 = 0;
    double pp = 0;

    // Boucle principale de la méthode de Newton
    do {
      double p1 = 1.0;
      double p2 = 0.;

      // Calcule le polynôme de Legendre de degré n au point z en utilisant la formule de récurrence
      for (int j = 0; j < n; ++j) {
        double const p3 = p2;
        p2 = p1;
        p1 = ((2. * j + 1.) * z * p2 - j * p3) / (j + 1);
      }

      // Dérivée du polynôme de Legendre au point z
      pp = n * (z * p1 - p2) / (z * z - 1);
      z1 = z;

      // Formule d'ajustement de la méthode de Newton
      z = z1 - p1 / pp;
    } while (fabs(z - z1) > 1e-12);

    points[i] = -z;
    points[n - i - 1] = z;
    weights[i] = 2. / ((1. - z * z) * pp * pp);
    weights[n - i - 1] = weights[i];
  }
}