DECLARE SUB LegendreCoeff (n%, A!())
DECLARE SUB HORNER (n%, A!(), X0!, K%, IB%, R!())
DECLARE FUNCTION Eval! (n%, A!(), x!)
'***************************************************
'*  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) =  .2231445                               *
'*                                                 *
'* ----------------------------------------------- *
'*             Basic Release By J-P Moreau, Paris. *
'*                      (www.jpmoreau.fr)          *
'***************************************************
'Program Eval_Leg

DEFINT I-N

OPTION BASE 0

DIM SHARED NMAX AS INTEGER

NMAX = 12

DIM A(NMAX)

DIM SHARED B(NMAX, NMAX)

  CLS
  PRINT
  INPUT " give order: ", n
  PRINT
  INPUT " give argument x: ", x

  LegendreCoeff n, A()

  PRINT
  PRINT " P(x) = "; Eval(n, A(), x)
  PRINT

END 'of main program

FUNCTION Eval (n, A(), x)
' Evaluate P(x) of order n for argument x 
DIM A1(NMAX), R(NMAX)
  FOR i = n TO 0 STEP -1
    A1(n - i + 1) = A(i)
  NEXT i
  HORNER n, A1(), x, 0, 1, R()
  Eval = R(1)
END FUNCTION

SUB HORNER (n, A(), X0, K, IB, 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
'     IB      FLAG
'             <>0  EVALUATION ONLY
'             = 0  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
'---------------------------------------------------------------------
  RR = A(1)
  FOR i = 1 TO K + 1
    R(i) = RR
  NEXT i
  FOR J = 2 TO n + 1
    R(1) = R(1) * X0 + A(J)
    NMJ = n - J + 1
    IF NMJ > K THEN
      L = K
    ELSE
      L = NMJ
    END IF
    FOR i = 2 TO L + 1
      R(i) = R(i) * X0 + R(i - 1)
    NEXT i
  NEXT J
  IF IB <> 0 THEN
    L = 1
    FOR i = 2 TO K + 1
      L = i * (i - 1)
      R(i) = R(i) * L
    NEXT i
  END IF
END SUB

'******************************************************
'* 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.                         *
'******************************************************
SUB LegendreCoeff (n, A())
' Establish p0 and p1 coefficients
  B(0, 0) = 1!: B(1, 0) = 0!: B(1, 1) = 1!
' Return if order is less then 2
  IF n > 1 THEN
    FOR i = 2 TO n
      B(i, 0) = -(i - 1) * B(i - 2, 0) / i
      FOR J = 1 TO i
        'Basic recursion relation
        B(i, J) = (i + i - 1) * B(i - 1, J - 1) - (i - 1) * B(i - 2, J)
        B(i, J) = B(i, J) / i
      NEXT J
    NEXT i
    FOR i = 0 TO n
      A(i) = B(n, i)
    NEXT i
  END IF
END SUB