! 例題 8-10  ルジャンドル多項式 ［再帰的関数 ］
RECURSIVE FUNCTION p( n, x ) RESULT ( pnx )

  IMPLICIT NONE
  INTEGER, INTENT ( IN ) :: n
  REAL, INTENT ( IN ) :: x
  REAL :: pnx

  IF ( n == 0 ) THEN
    pnx = 1.0E0
  ELSE IF ( n == 1 ) THEN
    pnx = x
  ELSE
    pnx = ( 2.0E0 - 1.0E0 / n ) * p( n - 1, x ) * x      &
        - ( 1.0E0 - 1.0E0 / n ) * p( n - 2, x )
  END IF

END FUNCTION p

! main
PROGRAM ex8_10

  IMPLICIT NONE
  INTEGER :: n, i
  REAL :: x, y, p
  CHARACTER ( LEN = 1 ) :: c( -25 : 25 )

  PRINT *, '何次の多項式にしますか？'
  READ *, n

  PRINT '( "Legendre polynomial p(", I2, ", x )" )', n

  DO i = 0, 20
    x = 0.05E0 * i
    y = p( n, x )
    CALL presentation()
    PRINT '( F5.2, 2X, F10.6, 1X, 51A1 )', x, y, c
  END DO

CONTAINS

  SUBROUTINE presentation( )          ! 内部手続きのため引数不要

    IMPLICIT NONE
    INTEGER :: j

    IF ( i /= 0 ) THEN
      c = ' '
    ELSE
      c = (/ ('-', j = -25, 25) /)
    END IF

    c( -25 ) = '|'
    c( 0 )   = '|'
    c( 25 )  = '|'
    c( NINT( 25 * y ) ) = '*'

  END SUBROUTINE presentation

END PROGRAM ex8_10