program polpak_prb
!
!*******************************************************************************
!
!! POLPAK_PRB calls the POLPAK test routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POLPAK_PRB'
  write ( *, '(a)' ) '  Tests for POLPAK, which computes the values of'
  write ( *, '(a)' ) '  certain special functions and polynomials.'
  write ( *, '(a)' ) ' '

  call test001
  call test002
  call test003
  call test004
  call test005
  call test006
  call test007
  call test008
  call test009
  call test010

  call test011
  call test012
  call test013
  call test014
  call test015
  call test016
  call test017
  call test018
  call test019
  call test020

  call test021
  call test022
  call test023
  call test024
  call test025
  call test026
  call test027
  call test028
  call test029
  call test030

  call test031
  call test032
  call test0325
  call test0326
  call test033
  call test034
  call test035
  call test036
  call test037
  call test038
  call test039
  call test0395
  call test040

  call test041
  call test042
  call test043
  call test044
  call test045
  call test046
  call test047
  call test048

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POLPAK_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test001
!
!*******************************************************************************
!
!! TEST001 tests ALIGN_ENUM.
!
  implicit none
!
  integer, parameter :: m_max = 10
  integer, parameter :: n_max = 10
!
  integer align_enum
  integer i
  integer j
  integer table(0:m_max,0:n_max)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST001'
  write ( *, '(a)' ) '  ALIGN_ENUM counts the number of possible'
  write ( *, '(a)' ) '  alignments of two biological sequences.'

  do i = 0, m_max
    do j = 0, n_max
      table(i,j) = align_enum ( i, j )
    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Alignment enumeration table:'
  write ( *, '(a)' ) ' '
  write ( *, '(2x,5i5,6i8)' ) ( i, i = 0, n_max )
  write ( *, '(a)' ) ' '
  do i = 0, m_max
    write ( *, '(i2,5i5,6i8)' ) i, table(i,0:n_max)
  end do

  return
end
subroutine test002
!
!*******************************************************************************
!
!! TEST002 tests BELL.
!
  implicit none
!
  integer, parameter :: n = 10
!
  integer b(0:n)
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST002'
  write ( *, '(a)' ) '  BELL computes the Bell numbers.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I,  BELL(I)'
  write ( *, '(a)' ) ' '

  call bell ( n, b )

  do i = 0, n
    write ( *, '(i4,2x,i10)' )  i, b(i)
  end do
 
  return
end
subroutine test003
!
!*******************************************************************************
!
!! TEST003 tests BENFORD.
!
  implicit none
!
  real benford
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST003'
  write ( *, '(a)' ) '  BENFORD(I) is the Benford probability of the'
  write ( *, '(a)' ) '  initial digit sequence I.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I,  BENFORD(I)'
  write ( *, '(a)' ) ' '

  do i = 1, 9
    write ( *, '(i4,2x,g14.6)' )  i, benford(i)
  end do
 
  return
end
subroutine test004
!
!*******************************************************************************
!
!! TEST004 tests D_FACTORIAL.
!
  implicit none
!
  double precision d_factorial
  real fx
  double precision fx2
  integer n
  real x
  integer x_copy
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST004:'
  write ( *, '(a)' ) '  D_FACTORIAL evaluates the factorial function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       D_FACTORIAL(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call r_factorial_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x <= 0.0E+00 ) then
      cycle
    end if

    x_copy = int ( x )
    fx2 = d_factorial ( x_copy )

    write ( *, '(f4.0,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test005
!
!*******************************************************************************
!
!! TEST005 tests BERN;
!! TEST005 tests BERN2;
!! TEST005 tests D_BERN3.
!
  implicit none
!
  integer, parameter :: n = 20
!
  real c1(0:n)
  real c2(0:n)
  double precision d_bern3
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST005'
  write ( *, '(a)' ) '  BERN computes Bernoulli numbers;'
  write ( *, '(a)' ) '  BERN2 computes Bernoulli numbers;'
  write ( *, '(a)' ) '  D_BERN3 computes Bernoulli numbers.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   I      B1               B2                B3'
  write ( *, '(a)' ) ' '
 
  call bern ( n, c1 )
  call bern2 ( n, c2 )
 
  do i = 0, n
    write ( *, '(i6,3g18.10)' ) i, c1(i), c2(i), d_bern3(i)
  end do
 
  return
end
subroutine test006
!
!*******************************************************************************
!
!! TEST006 tests BERN_POLY;
!! TEST006 tests BERN_POLY2.
!
  implicit none
!
  integer, parameter :: n = 15
!
  double precision bern_poly2
  real bx
  double precision bx2
  double precision dx
  integer i
  real x
!
  x = 0.2E+00
  dx = dble ( x )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST006'
  write ( *, '(a)' ) '  BERN_POLY evaluates Bernoulli polynomials;'
  write ( *, '(a)' ) '  BERN_POLY2 evaluates Bernoulli polynomials. '
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I          BX          BX2'
  write ( *, '(a)' ) ' '
 
  do i = 1, n
    call bern_poly ( i, x, bx )
    bx2 = bern_poly2 ( i, dx )
    write ( *, '(i2,2x,2g16.8)' ) i, bx, bx2
  end do
 
  return
end
subroutine test007
!
!*******************************************************************************
!
!! TEST007 tests BETA.
!
  implicit none
!
  real beta
  real fxy
  real fxy2
  integer n
  real x
  real y
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST007:'
  write ( *, '(a)' ) '  BETA evaluates the Beta function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X      Y        Exact F       BETA(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call beta_values ( n, x, y, fxy )

    if ( n == 0 ) then
      exit
    end if

    fxy2 = beta ( x, y )

    write ( *, '(2f8.4,2g14.6)' ) x, y, fxy, fxy2

  end do

  return
end
subroutine test008
!
!*******************************************************************************
!
!! TEST008 tests BP01.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real bern(0:n)
  integer i
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST008'
  write ( *, '(a)' ) '  BP01 evaluates Bernstein polynomials.'
  write ( *, '(a)' ) ' '
 
  x = 0.3E+00
 
  call bp01 ( n, bern, x )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  The Bernstein polynomials of degree ', n
  write ( *, '(a,g14.6)' ) '  at X = ', x
  write ( *, '(a)' ) ' '
 
  do i = 0, n
    write ( *, '(i4,2x,g14.6)' )  i, bern(i)
  end do
 
  return
end
subroutine test009
!
!*******************************************************************************
!
!! TEST009 tests BPAB.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real a
  real b
  real bern(0:n)
  integer i
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST009'
  write ( *, '(a)' ) '  BPAB evaluates Bernstein polynomials.'
  write ( *, '(a)' ) ' '

  x = 0.3E+00
  a = 0.0E+00
  b = 1.0E+00
  call bpab ( n, bern, x, a, b )
 
  write ( *, '(a,i4)' ) '  The Bernstein polynomials of degree ', n
  write ( *, '(a,g14.6)' ) '  based on the interval from ', a
  write ( *, '(a,g14.6)' ) '  to ', b
  write ( *, '(a,g14.6)' ) '  evaluated at X = ', x
  write ( *, '(a)' ) ' '
 
  do i = 0, n
    write ( *, '(i4,2x,g14.6)' )  i, bern(i)
  end do
 
  return
end
subroutine test010
!
!*******************************************************************************
!
!! TEST010 tests CARDAN.
!! TEST010 tests CARDAN_COEF.
!
  implicit none
!
  integer, parameter :: n_max = 10
!
  real c(0:n_max)
  real cx1
  real cx2(0:n_max)
  integer n
  real s
  real x
!
  s = 1.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST010'
  write ( *, '(a)' ) '  CARDAN_COEFF returns the coefficients of a'
  write ( *, '(a)' ) '    Cardan polynomial.'
  write ( *, '(a)' ) '  CARDAN evaluates a Cardan polynomial directly.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  We use the parameter S = ', s
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Table of polynomial coefficients:'
  write ( *, '(a)' ) ' '

  do n = 0, n_max
    call cardan_coeff ( n, s, c )
    write ( *, '(i2,11f7.0)' ) n, c(0:n)
  end do

  s = 0.5E+00
  x = 0.25E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Compare CARDAN_COEFF + RPOLY_VAL_HORNER versus CARDAN.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Evaluate polynomials at X = ', x
  write ( *, '(a,g14.6)' ) '  We use the parameter S = ', s
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Order, Horner, Direct'
  write ( *, '(a)' ) ' '

  call cardan ( n, x, s, cx2 )

  do n = 0, n_max

    call cardan_coeff ( n, s, c )
    call rpoly_val_horner ( n, c, x, cx1 )

    write ( *, '(i2,2g14.6)' ) n, cx1, cx2(n)

  end do

  return
end
subroutine test011
!
!*******************************************************************************
!
!! TEST011 tests CATALAN.
!
  implicit none
!
  integer, parameter :: n = 10
!
  integer i
  integer c(0:n)
  integer, dimension(0:n), parameter :: c2 = &
    (/ 1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796 /)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST011'
  write ( *, '(a)' ) '  CATALAN computes Catalan numbers.'
  write ( *, '(a)' ) ' '
 
  call catalan ( n, c )
 
  write ( *, '(a)' ) '  I, computed Cat(I), correct Cat(I)'
  write ( *, '(a)' ) ' '
 
  do i = 0, n
    write ( *, '(i3,2i10)' ) i, c(i), c2(i)
  end do
 
  return
end
subroutine test012
!
!*******************************************************************************
!
!! TEST012 tests CATALAN_ROW.
!
  implicit none
!
  integer, parameter :: n = 10
!
  integer c(0:n)
  integer i
  integer ido
  integer j
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST012'
  write ( *, '(a)' ) '  CATALAN_ROW computes a row of Catalan''s triangle.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  First, compute row 7:'

  ido = 0
  i = 7
  call catalan_row ( ido, i, c )
  write ( *, '(i2,2x,11i6)' ) i, c(0:i)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Now compute rows one at a time:'
  write ( *, '(a)' ) ' '

  ido = 0
 
  do i = 0, n
    call catalan_row ( ido, i, c )
    ido = 1
    write ( *, '(i2,2x,11i6)' ) i, c(0:i)
  end do
 
  return
end
subroutine test013
!
!*******************************************************************************
!
!! TEST013 tests CHEBY1.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real c(0:n)
  integer i
  real x
!
  x = 0.2E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST013'
  write ( *, '(a)' ) '  CHEBY1 evaluates Cheybyshev polynomials of the '
  write ( *, '(a)' ) '  first kind.'
  write ( *, '(a,g14.6)' ) '  Use X = ', x
  write ( *, '(a)' ) ' '
 
  call cheby1 ( n, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test014
!
!*******************************************************************************
!
!! TEST014 tests CHEBY2.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real c(0:n)
  integer i
  real x
!
  x = 0.2E+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST014'
  write ( *, '(a)' ) '  CHEBY2 evaluates Chebyshev polynomials of the '
  write ( *, '(a)' ) '  second kind.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)' ) ' '
 
  call cheby2 ( n, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test015
!
!*******************************************************************************
!
!! TEST015 tests COMBIN.
!
  implicit none
!
  real cnk
  integer k
  integer n
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST015'
  write ( *, '(a)' ) '  COMBIN evaluates C(N,K).'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   N     K    CNK'
  write ( *, '(a)' ) ' '
 
  do n = 0, 4
    do k = 0, n
      call combin ( n, k, cnk )
      write ( *, '(i6,i6,g14.6)' ) n, k, cnk
    end do
  end do
 
  return
end
subroutine test016
!
!*******************************************************************************
!
!! TEST016 tests COMBIN2.
!
  implicit none
!
  integer cnk
  integer k
  integer n
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST016'
  write ( *, '(a)' ) '  COMBIN2 evaluates C(N,K).'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   N     K    CNK'
  write ( *, '(a)' ) ' '
 
  do n = 0, 4
    do k = 0, n
      call combin2 ( n, k, cnk )
      write ( *, '(i6,i6,g14.6)' ) n, k, cnk
    end do
  end do
 
  return
end
subroutine test017
!
!*******************************************************************************
!
!! TEST017 tests COMB_ROW.
!
  implicit none
!
  integer, parameter :: n = 10
!
  integer c(0:n)
  integer i
  integer ido
  integer j
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST017'
  write ( *, '(a)' ) '  COMB_ROW computes a row of Pascal''s triangle.'
  write ( *, '(a)' ) ' '
 
  ido = 0
 
  do i = 0, n
    call comb_row ( ido, i, c )
    ido = 1
    write ( *, '(i2,2x,11i5)' ) i, c(0:i)
  end do
 
  return
end
subroutine test018
!
!*******************************************************************************
!
!! TEST018 tests EULER;
!! TEST018 tests D_EULER2.
!
  implicit none
!
  integer, parameter :: n = 16
!
  double precision d_euler2
  integer i
  integer e(0:n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST018'
  write ( *, '(a)' ) '  EULER computes Euler numbers;'
  write ( *, '(a)' ) '  D_EULER2 computes Euler numbers.'
  write ( *, '(a)' ) ' '

  call euler ( n, e )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'EULER, D_EULER2:'
  write ( *, '(a)' ) ' '
  do i = 0, n
    write ( *,'(i3,i12,g17.9)') i, e(i), d_euler2(i)
  end do

  return
end
subroutine test019
!
!*******************************************************************************
!
!! TEST019 calls EULERIAN for the Eulerian numbers.
!
  implicit none
!
  integer, parameter :: n = 7
!
  integer e(n,n)
  integer i
  integer j
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST019'
  write ( *, '(a)' ) '  EULERIAN evaluates Eulerian numbers.'
  write ( *, '(a)' ) ' '
 
  call eulerian ( n, e )

  do i = 1, n
    write ( *, '(10i6)' )  e(i,1:n)
  end do
 
  return
end
subroutine test020
!
!*******************************************************************************
!
!! TEST020 tests EULER_POLY.
!
  implicit none
!
  integer, parameter :: n = 15
!
  double precision euler_poly
  integer i
  double precision x
!
  x = 0.5E+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST020'
  write ( *, '(a)' ) '  EULER_POLY evaluates Euler polynomials.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Evaluate at X = ', x
  write ( *, '(a)' ) ' '
 
  do i = 0, n
    write ( *, '(i2,2x,g14.6)' )  i, euler_poly ( i, x )
  end do
 
  return
end
subroutine test021
!
!*******************************************************************************
!
!! TEST021 tests F_HOFSTADTER.
!
  implicit none
!
  integer f
  integer f_hofstadter
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST021'
  write ( *, '(a)' ) '  F_HOFSTADTER evaluates Hofstadter''s recursive'
  write ( *, '(a)' ) '  F function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     N   F(N)'
  write ( *, '(a)' ) ' '

  do i = 0, 30
    f = f_hofstadter ( i )
    write ( *, '(2i6)' ) i, f
  end do

  return
end
subroutine test022
!
!*******************************************************************************
!
!! TEST022 tests LOG_FACTORIAL.
!
  implicit none
!
  real fx
  real fx2
  real log_factorial
  integer n
  real x
  integer x_copy
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST022:'
  write ( *, '(a)' ) '  LOG_FACTORIAL evaluates the logarithm of the '
  write ( *, '(a)' ) '  factorial function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X	   Exact F	 LOG_FACTORIAL(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call r_factorial_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x <= 0.0E+00 ) then
      cycle
    end if

    x_copy = int ( x )
    fx2 = log_factorial ( x_copy )

    write ( *, '(f4.0,2g14.6)' ) x, log ( fx ), fx2

  end do

  return
end
subroutine test023
!
!*******************************************************************************
!
!! TEST023 tests FIBONACCI_DIRECT.
!
  implicit none
!
  integer f
  integer i
  integer n
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST023'
  write ( *, '(a)' ) '  FIBONACCI_DIRECT evalutes a Fibonacci number directly.'
  write ( *, '(a)' ) ' '
 
  n = 20
 
  do i = 1, n
    call fibonacci_direct ( i, f )
    write ( *, '(i6,i10)' ) i, f
  end do
 
  return
end
subroutine test024
!
!*******************************************************************************
!
!! TEST024 tests FIBONACCI_FLOOR.
!
  implicit none
!
  integer f
  integer i
  integer n
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST024'
  write ( *, '(a)' ) '  FIBONACCI_FLOOR computes the largest Fibonacci number'
  write ( *, '(a)' ) '  less than or equal to a given positive integer.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     N  Fibonacci  Index'
  write ( *, '(a)' ) ' ' 

  do n = 1, 20
    call fibonacci_floor ( n, f, i )
    write ( *, '(i6,2x,i6,2x,i6)' ) n, f, i
  end do
 
  return
end
subroutine test025
!
!*******************************************************************************
!
!! TEST025 tests FIBONACCI_RECURSIVE.
!
  implicit none
!
  integer, parameter :: n = 20
!
  integer f(n)
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST025'
  write ( *, '(a)' ) '  FIBONACCI_RECURSIVE computes the Fibonacci sequence.'
  write ( *, '(a)' ) ' '
 
  call fibonacci_recursive ( n, f )
 
  do i = 1, n
    write ( *, '(i6,i10)' ) i, f(i)
  end do
 
  return
end
subroutine test026
!
!*******************************************************************************
!
!! TEST026 tests G_HOFSTADTER.
!
  implicit none
!
  integer g
  integer g_hofstadter
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST026'
  write ( *, '(a)' ) '  G_HOFSTADTER evaluates Hofstadter''s recursive'
  write ( *, '(a)' ) '  G function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     N   G(N)'
  write ( *, '(a)' ) ' '

  do i = 0, 30
    g = g_hofstadter ( i )
    write ( *, '(2i6)' ) i, g
  end do

  return
end
subroutine test027
!
!*******************************************************************************
!
!! TEST027 tests GAMMA.
!
  implicit none
!
  real fx
  real fx2
  real gamma
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST027:'
  write ( *, '(a)' ) '  GAMMA evaluates the Gamma function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       GAMMA(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call gamma_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    fx2 = gamma ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test028
!
!*******************************************************************************
!
!! TEST028 tests GAMMA_LOG.
!
  implicit none
!
  real fx
  real fx2
  real gamma_log
  integer n
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST028:'
  write ( *, '(a)' ) '  GAMMA_LOG evaluates the logarithm of the '
  write ( *, '(a)' ) '  Gamma function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       GAMMA_LOG(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call gamma_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    fx = log ( fx )
    fx2 = gamma_log ( x )

    write ( *, '(f8.4,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test029
!
!*******************************************************************************
!
!! TEST029 tests GEGENBAUER.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real alfa
  real c(0:n)
  integer i
  real x
!
  alfa = 0.5E+00
  x = 1.2E+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST029'
  write ( *, '(a)' ) '  GEGENBAUER evaluates Gegenbauer polynomials.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Alfa = ', alfa 
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)' ) ' '
 
  call gegenbauer ( n, alfa, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test030
!
!*******************************************************************************
!
!! TEST030 tests HAIL.
!
  implicit none
!
  integer hail
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST030'
  write ( *, '(a)' ) '  HAIL(I) computes the length of the hail sequence'
  write ( *, '(a)' ) '  for I, also known as the 3*N+1 sequence.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I,  HAIL(I)'
  write ( *, '(a)' ) ' '

  do i = 1, 20
    write ( *, '(i4,2x,i6)' )  i, hail(i)
  end do
 
  return
end
subroutine test031
!
!*******************************************************************************
!
!! TEST031 tests H_HOFSTADTER.
!
  implicit none
!
  integer h
  integer h_hofstadter
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST031'
  write ( *, '(a)' ) '  H_HOFSTADTER evaluates Hofstadter''s recursive'
  write ( *, '(a)' ) '  H function.'

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     N   H(N)'
  write ( *, '(a)' ) ' '

  do i = 0, 30
    h = h_hofstadter ( i )
    write ( *, '(2i6)' ) i, h
  end do

  return
end
subroutine test032
!
!*******************************************************************************
!
!! TEST032 tests HERMITE.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real c(0:n)
  integer i
  real x
!
  x = 0.5E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST032' 
  write ( *, '(a)' ) '  HERMITE evaluates the Hermite polynomials.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)') ' '

  call hermite ( n, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
  
  return
end
subroutine test0325
!
!*******************************************************************************
!
!! TEST0325 tests I_FACTORIAL.
!
  implicit none
!
  real fx
  integer fx2
  integer i_factorial
  integer n
  real x
  integer x_copy
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST0325:'
  write ( *, '(a)' ) '  I_FACTORIAL evaluates the factorial function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       I_FACTORIAL(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call r_factorial_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x <= 0.0E+00 ) then
      cycle
    end if

    x_copy = int ( x )
    fx2 = i_factorial ( x_copy )

    write ( *, '(f4.0,g14.6,i14)' ) x, fx, fx2

  end do

  return
end
subroutine test0326
!
!*******************************************************************************
!
!! TEST0326 tests I_FACTORIAL2.
!
  implicit none
!
  integer fx
  integer i
  integer i_factorial2
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST0326:'
  write ( *, '(a)' ) '  I_FACTORIAL2 evaluates the factorial function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       I_FACTORIAL2(X)'
  write ( *, '(a)' ) ' '

  do i = 0, 16

    fx = i_factorial2 ( i )

    write ( *, '(i4,i14)' ) i, fx

  end do

  return
end
subroutine test033
!
!*******************************************************************************
!
!! TEST033 tests JACOBI.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real alfa
  real beta
  real c(0:n)
  integer i
  real x
!
  alfa = -0.5E+00
  beta = 0.5E+00
  x = 0.5E+00
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST033'
  write ( *, '(a)' ) '  JACOBI evaluates the Jacobi polynomials.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Alfa = ', alfa
  write ( *, '(a,g14.6)' ) '  Beta = ', beta
  write ( *, '(a,g14.6)' ) '  X =    ', x
  write ( *, '(a)' ) ' '
 
  call jacobi ( n, alfa, beta, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test034
!
!*******************************************************************************
!
!! TEST034 tests LAGUERRE_GEN.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real alfa
  real c(0:n)
  integer i
  real x
!
  x = 0.5E+00
  alfa = 0.1E+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST034'
  write ( *, '(a)' ) '  LAGUERRE_GEN evaluates the generalized Laguerre '
  write ( *, '(a)' ) '  polynomials.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  ALFA = ', alfa
  write ( *, '(a,g14.6)' ) '  X =    ', x
  write ( *, '(a)' ) ' '
 
  call laguerre_gen ( n, alfa, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test035
!
!*******************************************************************************
!
!! TEST035 tests LAGUERRE_LNM.
!
  implicit none
!
  integer, parameter :: n = 5
  integer, parameter :: m = 1
!
  real c(0:n)
  integer i
  real x
!
  x = 0.5E+00
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST035'
  write ( *, '(a)' ) '  LAGUERRE_LNM evaluates the associated Laguerre'
  write ( *, '(a)' ) '  polynomials,'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  M = ', m
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)' ) ' '
 
  call laguerre_lnm ( n, m, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test036
!
!*******************************************************************************
!
!! TEST036 tests LAGUERRE.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real c(0:n)
  integer i
  real x
!
  x = 0.5E+00
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST036'
  write ( *, '(a)' ) '  LAGUERRE evaluates the Laguerre polynomials.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)' ) ' '
 
  call laguerre ( n, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test037
!
!*******************************************************************************
!
!! TEST037 tests LEGENDRE_PN.
!
  implicit none
!
  integer, parameter :: n = 10
!
  real c(0:n)
  real cp(0:n)
  integer i
  real x
!
  x = 0.5E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST037'
  write ( *, '(a)' ) '  LEGENDRE_PN evaluates the Legendre polynomials of'
  write ( *, '(a)' ) '  the first kind, and derivatives.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X =', x
  write ( *, '(a)' ) ' '
 
  call legendre_pn ( n, x, c, cp )
 
  do i = 0, n
    write ( *, '(i6,2g14.6)' ) i, c(i), cp(i)
  end do
  
  return
end
subroutine test038
!
!*******************************************************************************
!
!! TEST038 tests LEGENDRE_PNM.
!
  implicit none
!
  integer a
  integer b
  real cx(0:20)
  real fx
  real fx2
  integer i
  integer n
  real plgndr
  real x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST038:'
  write ( *, '(a)' ) '  LEGENDRE_PNM evaluates associated Legrendre functions.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      N       M    X     Exact F     PNM(X)'
  write ( *, '(a)' ) ' '
  n = 0

  do

    call legendre_function_values ( n, a, b, x, fx )

    if ( n == 0 ) then
      exit
    end if

    call legendre_pnm ( a, b, x, cx )

    fx2 = cx(a)

    write ( *, '(i8,i8,f8.4,2g14.6)' ) a, b, x, fx, fx2

  end do

  return
end
subroutine test039
!
!*******************************************************************************
!
!! TEST039 tests LEGENDRE_QN.
!
  implicit none
!
  integer, parameter :: n = 4
!
  real c(0:n)
  integer i
  real x
!
  x = 0.5E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST039'
  write ( *, '(a)' ) '  LEGENDRE_QN evaluates Legendre polynomials '
  write ( *, '(a)' ) '  of the second kind.' 
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  X = ', x
  write ( *, '(a)' ) ' '
 
  call legendre_qn ( n, x, c )
 
  do i = 0, n
    write ( *, '(i6,g14.6)' ) i, c(i)
  end do
 
  return
end
subroutine test0395
!
!*******************************************************************************
!
!! TEST0395 tests LOCK.
!
  implicit none
!
  integer, parameter :: n = 10
!
  integer a(0:n)
  integer i
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST0395'
  write ( *, '(a)' ) '  LOCK counts the combinations on a button lock.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I,  LOCK(I)'
  write ( *, '(a)' ) ' '

  call lock ( n, a )

  do i = 0, n
    write ( *, '(i4,2x,i10)' )  i, a(i)
  end do
 
  return
end
subroutine test040
!
!*******************************************************************************
!
!! TEST040 tests PENTAGON_NUM.
!
  implicit none
!
  integer n
  integer p
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST040'
  write ( *, '(a)' ) '  PENTAGON_NUM computes the pentagonal numbers.'
  write ( *, '(a)' ) ' '
 
  do n = 1, 10
    call pentagon_num ( n, p )
    write ( *, '(i4,2x,i6)' ) n, p
  end do
 
  return
end
subroutine test041
!
!*******************************************************************************
!
!! TEST041 tests PYRAMID_NUM.
!
  implicit none
!
  integer n
  integer pyramid_num
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST041'
  write ( *, '(a)' ) '  PYRAMID_NUM computes the pyramidal numbers.'
  write ( *, '(a)' ) ' '
 
  do n = 1, 10
    write ( *, '(i4,2x,i6)' ) n, pyramid_num ( n )
  end do
 
  return
end
subroutine test042
!
!*******************************************************************************
!
!! TEST042 tests R_FACTORIAL.
!
  implicit none
!
  real fx
  real fx2
  integer n
  real r_factorial
  real x
  integer x_copy
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST042:'
  write ( *, '(a)' ) '  R_FACTORIAL evaluates the factorial function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     X       Exact F       R_FACTORIAL(X)'
  write ( *, '(a)' ) ' '

  n = 0

  do

    call r_factorial_values ( n, x, fx )

    if ( n == 0 ) then
      exit
    end if

    if ( x <= 0.0E+00 ) then
      cycle
    end if

    x_copy = int ( x )
    fx2 = r_factorial ( x_copy )

    write ( *, '(f4.0,2g14.6)' ) x, fx, fx2

  end do

  return
end
subroutine test043
!
!*******************************************************************************
!
!! TEST043 tests STIRLING1.
!
  implicit none
!
  integer, parameter :: m = 8
  integer, parameter :: n = m
!
  integer i
  integer j
  integer s1(m,n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST043'
  write ( *, '(a)' ) '  STIRLING1: Stirling numbers of first kind.'
  write ( *, '(a,i4)' ) '  Get rows 1 through ', m
  write ( *, '(a)' ) ' '
 
  call stirling1 ( m, n, s1 )
 
  do i = 1, m
    write ( *, '(i6,8i8)' ) i, s1(i,1:n)
  end do
 
  return
end
subroutine test044
!
!*******************************************************************************
!
!! TEST044 tests STIRLING2.
!
  implicit none
!
  integer, parameter :: m = 8
  integer, parameter :: n = m
!
  integer i
  integer j
  integer s2(m,n)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST044'
  write ( *, '(a)' ) '  STIRLING2: Stirling numbers of second kind.'
  write ( *, '(a,i4)' ) '  Get rows 1 through ', m
  write ( *, '(a)' ) ' '
 
  call stirling2 ( m, n, s2 )
 
  do i = 1, m
    write ( *, '(i6,8i8)' ) i, s2(i,1:n)
  end do
 
  return
end
subroutine test045
!
!*******************************************************************************
!
!! TEST045 tests TRIANGLE_NUM.
!
  implicit none
!
  integer n
  integer triangle_num
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST045'
  write ( *, '(a)' ) '  TRIANGLE_NUM computes the triangular numbers.'
  write ( *, '(a)' ) ' '
 
  do n = 1, 10
    write ( *, '(i4,2x,i6)' ) n, triangle_num ( n )
  end do
 
  return
end
subroutine test046
!
!*******************************************************************************
!
!! TEST046 tests V_HOFSTADTER.
!
  implicit none
!
  integer i
  integer v
  integer v_hofstadter
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST046'
  write ( *, '(a)' ) '  V_HOFSTADTER evaluates Hofstadter''s recursive'
  write ( *, '(a)' ) '  V function.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     N   V(N)'
  write ( *, '(a)' ) ' '

  do i = 0, 30
    v = v_hofstadter ( i )
    write ( *, '(2i6)' ) i, v
  end do

  return
end
subroutine test047
!
!*******************************************************************************
!
!! TEST047 tests VIBONACCI.
!
  implicit none
!
  integer, parameter :: n = 20
  integer, parameter :: n_time = 3
!
  integer i
  integer j
  integer v(n,n_time)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST047'
  write ( *, '(a)' ) '  VIBONACCI computes a Vibonacci sequence.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i4)' ) '  Number of times we compute the series: ', n_time
  write ( *, '(a)' ) ' '

  do j = 1, n_time
    call vibonacci ( n, v(1,j) ) 
  end do

  do i = 1, n
    write ( *, '(i6,3i6)' ) i, v(i,1:n_time)
  end do
 
  return
end
subroutine test048
!
!*******************************************************************************
!
!! TEST048 tests ZECKENDORF.
!
  implicit none
!
  integer i
  integer i_list(20)
  integer f_list(20)
  integer f_sum
  integer m
  integer n
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST048'
  write ( *, '(a)' ) '  ZECKENDORF computes the Zeckendorf decomposition of'
  write ( *, '(a)' ) '  an integer into nonconsecutive Fibonacci numbers.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   N Sum M Parts'
  write ( *, '(a)' ) ' '

  do n = 1, 100

    call zeckendorf ( n, m, i_list, f_list )

    write ( *, '(i4,2x,15i4)' ) n, f_list(1:m)

  end do

  return
end