/**************************************************
*  Evaluate a Legendre Polynomial P(x) of Order n *
*  for argument x by using Horner's Rule          *
* ----------------------------------------------- *
* SAMPLE RUN:                                     *
*                                                 *
*  give order: 7                                  *
*                                                 *
*  give argument x: 0.5                           *
*                                                 *
*  P(x) =  0.223145                               *
*                                                 *
* ----------------------------------------------- *
*               C++ Release By J-P Moreau, Paris. *
*                      (www.jpmoreau.fr)          *
**************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  12
#define  TRUE  1


/*****************************************************
* Legendre series coefficients evaluation subroutine *
* by means of recursion relation. The order of the   *
* polynomial is n. The coefficients are returned in  *
* A(i) by increasing powers.                         *
*****************************************************/
void Legendre_Coeff(int n, double *A) {
  int i,j;
  double B[NMAX][NMAX];
// Establish p0 and p1 coefficients
  B[0][0]=1.0; B[1][0]=0.0; B[1][1]=1.0;
// Return if order is less then 2
  if (n > 1) {
	for (i=2; i<=n; i++) {
      B[i][0] = -(i-1)*B[i-2][0]/i;
      for (j=1; j<=i; j++) {
        // Basic recursion relation
        B[i][j]=(i+i-1)*B[i-1][j-1]-(i-1)*B[i-2][j];
        B[i][j] /= i;
      }
    }
    for (i=0; i<=n; i++)  A[i]=B[n][i];
  }
}

void HORNER (int N, double *A, double X0, int K, bool B, double *R) {
/*--------------------------------------------------------------------
!     EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER'S METHOD
!
!     INPUTS:
!     N       ORDER OF POLYNOMIAL
!     A       VECTOR OF SIZE N+1 STORING THE COEFFICIENTS OF
!             POLYNOMIAL IN DECREASING ORDERS OF POWERS
!             I.E. P(X) = A(1)*X**N+A(2)*X**(N-1)+...+A(N)*X+A(N+1)
!     X0      GIVEN ARGUMENT
!     K       MAXIMUM ORDER OF DERIVATIVES ASKED FOR
!     B       FLAG
!             = .TRUE.  EVALUATION ONLY
!             = .FALSE. DETERMINATION OF POLYNOMIAL COEFFICIENTS
!                NEAR X0
!     R       VECTOR OF SIZE K+1 CONTAINS:
!             -  VALUES P(X0), P'(X0), P''(X0),..PK(X0),
!                IF B IS TRUE AND IF K < N
!             -  THE N+1 COEFFICIENTS OF POLYNOMIAL
!                Q(X-X0) = R(1)+R(2)*(X-X0)+...+R(N+1)*(X-X0)**(N),
!                IF B IS FALSE AND IF K = N.
!
!     REFERENCE:
!     ALGORITHM 337, COLLECTED ALGORITHMS FROM CACM, W.PANKIEWICKZ
!-------------------------------------------------------------------*/
  double RR;
  int I,J,L,NMJ;

  RR = A[1];
  for (I=1; I<=K+1; I++) R[I] = RR;
  for (J=2; J<=N+1; J++) {
    R[1] = R[1]*X0 + A[J];
    NMJ = N-J+1;
    if (NMJ > K)
      L = K;
    else
      L = NMJ;
    for (I=2; I<=L+1; I++) R[I] = R[I]*X0 + R[I-1];
  }
  if (B) {
    L = 1;
    for (I=2; I<=K+1; I++) {
      L *= (I-1);
      R[I] *= L;
    }
  }
}

// Evaluate P(x) of order n for argument x
double Eval(int n, double *A, double x) {
  int i;
  double A1[NMAX], R[NMAX];
  for (i=n; i>=0; i--)  A1[n-i+1] = A[i];
  HORNER(n,A1,x,0,TRUE,R);
  return R[1];
}


void main() {

  double A[NMAX];
  int n;
  double x;

  printf("\n give order: "); scanf("%d", &n);
  printf("\n give argument x: "); scanf("%lf", &x);

  Legendre_Coeff(n, A);

  printf("\n P(x) = %f\n\n", Eval(n,A,x));

}

// end of file eval_leg.cpp