#! /usr/bin/env python
#
def pn_polynomial_coefficients ( n ):

#*****************************************************************************80
#
## PN_POLYNOMIAL_COEFFICIENTS: coefficients of normalized Legendre Pn(n,x).
#
#  Discussion:
#
#    Pn(n,x) = P(n,x) * sqrt ( (2n+1)/2 )
#
#          1       x       x^2     x^3     x^4      x^5    x^6     x^7
#
#    0   0.707
#    1   0.000   1.224
#    2  -0.790   0.000   2.371
#    3   0.000  -2.806   0.000   4.677
#    4   0.795   0.000  -7.954   0.000   9.280
#    5   0.000   4.397   0.000 -20.520   0.000   18.468
#    6  -0.796   0.000  16.731   0.000 -50.193    0.000  36.808
#    7   0.000  -5.990   0.000  53.916   0.000 -118.616   0.000  73.429 
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    17 March 2016
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Milton Abramowitz and Irene Stegun,
#    Handbook of Mathematical Functions,
#    US Department of Commerce, 1964.
#
#    Daniel Zwillinger, editor,
#    CRC Standard Mathematical Tables and Formulae,
#    30th Edition,
#    CRC Press, 1996.
#
#  Parameters:
#
#    Input, integer N, the highest order polynomial to evaluate.
#    Note that polynomials 0 through N will be evaluated.
#
#    Output, real C(1:N+1,1:N+1), the coefficients of the normalized Legendre 
#    polynomials of degree 0 through N.  Each polynomial is stored as a row.
#
  import numpy as np

  c = np.zeros ( [ n + 1, n + 1 ] )
#
#  Compute P(n,x) coefficients.
#
  c[0,0] = 1.0

  if ( 0 < n ):
    c[1,1] = 1.0

  for i in range ( 2, n + 1 ):
    for j in range ( 0, i ):
      c[i,j] = (   - i + 1 ) * c[i-2,j] / ( i )
    for j in range ( 1, i + 1 ):
      c[i,j] = c[i,j] + ( i + i - 1 ) * c[i-1,j-1] / ( i )
#
#  Normalize them.
#
  for i in range ( 0, n + 1 ):
    t = np.sqrt ( ( 2 * i + 1 ) / 2.0 )
    for j in range ( 0, i + 1 ):
      c[i,j] = c[i,j] * t

  return c

def pn_polynomial_coefficients_test ( ):

#*****************************************************************************80
#
## PN_POLYNOMIAL_COEFFICIENTS_TEST tests PN_POLYNOMIAL_COEFFICIENTS.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    17 March 2016
#
#  Author:
#
#    John Burkardt
#
  import platform

  n = 10

  print ( '' )
  print ( 'PN_POLYNOMIAL_COEFFICIENTS_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  PN_POLYNOMIAL_COEFFICIENTS: polynomial coefficients of Pn(n,x).' )

  c = pn_polynomial_coefficients ( n )
 
  for i in range ( 0, n + 1 ):
    print ( '' )
    print ( '  P(%d,x) = ' % ( i ) )
    print ( '' )
    for j in range ( i, -1, -1 ):
      if ( c[i,j] == 0.0 ):
        continue
      elif ( j == 0 ):
        print ( '  %g' % ( c[i,j] ) )
      elif ( j == 1 ):
        print ( '  %g * x' % ( c[i,j] ) )
      else:
        print ( '  %g * x^%d' % ( c[i,j], j ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'PN_POLYNOMIAL_COEFFICIENTS_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( ) 
  pn_polynomial_coefficients_test ( )
  timestamp ( )