program main

c*********************************************************************72
c
cc hyper_2f1_test() tests hyper_2f1().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 December 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'hyper_2f1_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test hyper_2f1().'

      call hyper_2f1_real_test ( )
      call hyper_2f1_complex_test ( )
      call psi_values_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'hyper_2f1_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine hyper_2f1_real_test ( )

c*********************************************************************72
c
cc hyper_2f1_real_test() tests hyper_2f1_real().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 December 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision a
      double precision b
      double precision c
      double precision fx1
      double precision fx2
      integer n_data
      double precision x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'hyper_2f1_real_test():'
      write ( *, '(a)' ) '  Test hyper_2f1_real() for real arguments.'

      n_data = 0

10    continue

        call hyper_2f1_values ( n_data, a, b, c, x, fx1 )

        if ( n_data == 0 ) then
          go to 20
        end if

        call hyper_2f1_real ( a, b, c, x, fx2 )

        write ( *, '(a)' ) ''
        write ( *, '(a,2x,f10.6,2x,f10.6,2x,f10.6,2x,f10.6)' ) 
     &    '  a,b,c,x: ', a, b, c, x
        write ( *, '(a,2x,g14.6)' ) 
     &    '  exact:   ', fx1
        write ( *, '(a,2x,g14.6)' ) 
     &    '  computed:', fx2
        write ( *, '(a,2x,g14.6)' ) 
     &    '  error:   ', abs ( fx1 - fx2 )

      go to 10

20    continue

      return
      end
      subroutine hyper_2f1_values ( n_data, a, b, c, x, fx )

c*********************************************************************72
c
cc hyper_2f1_values() returns some values of the hypergeometric function 2F1.
c
c  Discussion:
c
c    In Mathematica, the function can be evaluated by:
c
c      fx = Hypergeometric2F1 [ a, b, c, x ]
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    09 September 2007
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Cambridge University Press, 1999,
c    ISBN: 0-521-64314-7,
c    LC: QA76.95.W65.
c
c    Shanjie Zhang, Jianming Jin,
c    Computation of Special Functions,
c    Wiley, 1996,
c    ISBN: 0-471-11963-6,
c    LC: QA351.C45
c
c    Daniel Zwillinger, editor,
c    CRC Standard Mathematical Tables and Formulae,
c    30th Edition,
c    CRC Press, 1996,
c    ISBN: 0-8493-2479-3,
c    LC: QA47.M315.
c
c  Input:
c
c    integer N_DATA: the user sets N_DATA to 0 on the first call.
c
c  Output:
c
c    integer N_DATA.  The routine increments N_DATA by 1,
c    and returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    double precision A, B, C: the parameters.
c
c    double precision X: the argument.
c
c    double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 24 )

      double precision a
      double precision a_vec(n_max)
      double precision b
      double precision b_vec(n_max)
      double precision c
      double precision c_vec(n_max)
      double precision fx
      double precision fx_vec(n_max)
      integer n_data
      double precision x
      double precision x_vec(n_max)

      save a_vec
      save b_vec
      save c_vec
      save fx_vec
      save x_vec

      data a_vec /
     &   -2.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    2.5D+00,
     &   -2.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    2.5D+00,
     &   -2.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    2.5D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00 /
      data b_vec /
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00 /
      data c_vec /
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &   -5.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    4.5D+00,
     &   -5.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    4.5D+00,
     &   -5.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    4.5D+00 /
      data fx_vec /
     &    0.72356129348997784913D+00,
     &    0.97911109345277961340D+00,
     &    1.0216578140088564160D+00,
     &    1.4051563200112126405D+00,
     &    0.46961431639821611095D+00,
     &    0.95296194977446325454D+00,
     &    1.0512814213947987916D+00,
     &    2.3999062904777858999D+00,
     &    0.29106095928414718320D+00,
     &    0.92536967910373175753D+00,
     &    1.0865504094806997287D+00,
     &    5.7381565526189046578D+00,
     &    15090.669748704606754D+00,
     &   -104.31170067364349677D+00,
     &    21.175050707768812938D+00,
     &    4.1946915819031922850D+00,
     &    1.0170777974048815592D+10,
     &   -24708.635322489155868D+00,
     &    1372.2304548384989560D+00,
     &    58.092728706394652211D+00,
     &    5.8682087615124176162D+18,
     &   -4.4635010147295996680D+08,
     &    5.3835057561295731310D+06,
     &    20396.913776019659426D+00 /
      data x_vec /
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00 /

      if ( n_data .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. n_data ) then
        n_data = 0
        a = 0.0D+00
        b = 0.0D+00
        c = 0.0D+00
        x = 0.0D+00
        fx = 0.0D+00
      else
        a = a_vec(n_data)
        b = b_vec(n_data)
        c = c_vec(n_data)
        x = x_vec(n_data)
        fx = fx_vec(n_data)
      end if

      return
      end
      subroutine hyper_2f1_complex_test ( )

c*********************************************************************72
c
cc hyper_2f1_complex_test() tests hyper_2f1_complex().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    21 December 2023
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision a
      double precision b
      double precision c
      double complex fz1
      double complex fz2
      integer n_data
      double complex z

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'hyper_2f1_complex_test():'
      write ( *, '(a)' ) 
     &  '  Test hyper_2f1_complex() for complex arguments.'
      write ( *, '(a)' ) ' '

      n_data = 0

10    continue

        call hyper_2f1_complex_values ( n_data, a, b, c, z, fz1 )

        if ( n_data == 0 ) then
          go to 20
        end if

        call hyper_2f1_complex ( a, b, c, z, fz2 )

        write ( *, '(a)' ) ''
        write ( *, '(a,2x,f10.6,2x,f10.6,2x,f10.6,2x,2f10.6)' ) 
     &    '  a,b,c,x: ', a, b, c, z
        write ( *, '(a,2x,2g14.6)' ) 
     &    '  exact:   ', fz1
        write ( *, '(a,2x,2g14.6)' ) 
     &    '  computed:', fz2
        write ( *, '(a,2x,g14.6)' ) 
     &    '  error:   ', cdabs ( fz1 - fz2 )

      go to 10

20    continue

      return
      end
      subroutine hyper_2f1_complex_values ( n_data, a, b, c, z, fz )

c*********************************************************************72
c
cc hyper_2f1_complex_values() returns some values of the hypergeometric function 2F1.
c
c  Discussion:
c
c    In Mathematica, the function can be evaluated by:
c
c      fz = Hypergeometric2F1 [ a, b, c, z ]
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    22 December 2023
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Cambridge University Press, 1999,
c    ISBN: 0-521-64314-7,
c    LC: QA76.95.W65.
c
c    Shanjie Zhang, Jianming Jin,
c    Computation of Special Functions,
c    Wiley, 1996,
c    ISBN: 0-471-11963-6,
c    LC: QA351.C45
c
c    Daniel Zwillinger, editor,
c    CRC Standard Mathematical Tables and Formulae,
c    30th Edition,
c    CRC Press, 1996,
c    ISBN: 0-8493-2479-3,
c    LC: QA47.M315.
c
c  Input:
c
c    integer N_DATA: The user sets N_DATA to 0 before the first call.
c
c  Output:
c
c    integer N_DATA: On each call, the routine increments N_DATA by 1,
c    and returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    double precision A, B, C: the parameters.
c
c    double complex Z: the argument.
c
c    double complex FZ: the value of the function.
c
      implicit none

      integer, parameter :: n_max = 15

      double precision a
      double precision a_vec(n_max)
      double precision b
      double precision b_vec(n_max)
      double precision c
      double precision c_vec(n_max)
      double complex fz
      double complex fz_vec(n_max)
      integer n_data
      double complex z
      double complex z_vec(n_max)

      save a_vec
      save b_vec
      save c_vec
      save fz_vec
      save z_vec

      data a_vec /
     &   3.2,
     &   3.2,
     &  -5.0,
     &   3.3,
     &  -7.0,
     &   4.3,
     &   3.3,
     &   3.5,
     &   3.3,
     &   7.0,
     &   5.0,
     &   3.5,
     &   2.1,
     &   8.7,
     &   8.7 /
      data b_vec /
     &   1.8,
     &  -1.8,
     &   3.3,
     &  -6.0,
     &   3.3,
     &  -8.0,
     &   5.8,
     &  -2.4,
     &   4.3,
     &   5.0,
     &   7.0,
     &   1.2,
     &   5.4,
     &   3.2,
     &   2.7 /
      data c_vec /
     &   6.7,
     &   6.7,
     &   6.7,
     &   3.7,
     &  -3.7,
     &  -3.7,
     &   6.7,
     &   6.7,
     &   6.7,
     &   4.1,
     &   4.1,
     &   9.7,
     &   9.7,
     &   6.7,
     &   6.7 /
      data fz_vec /
     & (  5.468999154361234D+00,    +0.00000000D+00 ), 
     &  (  0.3375063477462785D+00,   +0.00000000D+00 ), 
     &  (  116.8274991533609D+00,    +603.8909562709345D+00 ), 
     &  (  17620.41819334182D+00,    +38293.80901310932D+00 ), 
     &  ( -11772775115.27448D+00,    -14382285977.20268D+00 ), 
     &  (  1316118577866.058D+00,    -101298889382.4362D+00 ), 
     &  (  1.733055678355656D+00,    +0.6340102904953357D+00 ), 
     &  (  0.6476224071999852D+00,   -0.5211050690999773D+00 ), 
     &  ( -1.483008322270093D+00,    +8.374426179451589D+00 ), 
     &  ( -0.004037609523971226D+00, -0.002956632645480181D+00 ), 
     &  ( -0.004037609523971226D+00, -0.002956632645480181D+00 ), 
     &  (  1.034313610729953D+00,    +0.5447389238499308D+00 ), 
     &  (  0.6885043978280027D+00,   +1.227418679098749D+00 ), 
     &  ( -0.9004649679297319D+00,   -1.11988994714304D+00 ), 
     &  ( -0.4608388640599718D+00,   -0.5457569650549665D+00 ) /
      data z_vec /
     &  ( 1.0, +0.0 ),
     &  ( 1.0, +0.0 ),
     &  ( 5.2, +4.8 ),
     &  ( 5.2, -4.8 ),
     &  ( 5.2, -4.8 ),
     &  ( 5.2, +4.8 ),
     &  ( 0.2, +0.1 ),
     &  ( 0.2, +0.5 ),
     &  ( 0.8, +0.3 ),
     &  ( 3.0, -1.0 ),
     &  ( 3.0, -1.0 ),
     &  ( 0.6, +0.9 ),
     &  ( 0.5, +0.7 ),
     &  ( 0.5, +0.7 ),
     &  ( 0.6, +0.9 ) /

      if ( n_data < 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max < n_data ) then
        n_data = 0
        a = 0.0D+00
        b = 0.0D+00
        c = 0.0D+00
        z = 0.0D+00
        fz = ( 0.0D+00, 0.0D+00 )
      else
        a = a_vec(n_data)
        b = b_vec(n_data)
        c = c_vec(n_data)
        z = z_vec(n_data)
        fz = fz_vec(n_data)
      end if

      return
      end
      subroutine psi_values_test ( )

c*********************************************************************72
c
cc psi_values_test() tests psi_values().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    04 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision fx1
      double precision fx2
      integer n_data
      double precision x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'psi_values_test():'
      write ( *, '(a)' ) '  psi_values() returns values of '
      write ( *, '(a)' ) '  the Psi function.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '          X               PSI(X)'
      write ( *, '(a)' ) ' '

      n_data = 0

10    continue

        call psi_values ( n_data, x, fx1 )

        if ( n_data .eq. 0 ) then
          go to 20
        end if

        call psi ( x, fx2 )

        write ( *, '(2x,f12.8,2x,g24.16,2x,g24.16)' ) x, fx1, fx2

      go to 10

20    continue

      return
      end
      subroutine psi_values ( n_data, x, fx )

c*********************************************************************72
c
cc psi_values() returns some values of the Psi or Digamma function for testing.
c
c  Discussion:
c
c    PSI(X) = d LN ( GAMMA ( X ) ) / d X = GAMMA'(X) / GAMMA(X)
c
c    PSI(1) = - Euler's constant.
c
c    PSI(X+1) = PSI(X) + 1 / X.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    04 June 2013
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c  Input:
c
c    integer N_DATA: the user sets N_DATA to 0 before the first call.  
c
c  Output:
c
c    integer N_DATA.  The routine increments N_DATA by 1, and
c    returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    double precision X, the argument of the function.
c
c    double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 20 )

      double precision fx
      double precision fxvec ( n_max )
      integer n_data
      double precision x
      double precision xvec ( n_max )

      data fxvec /
     & -10.42375494041108D+00, 
     &  -5.289039896592188D+00, 
     &  -3.502524222200133D+00, 
     &  -2.561384544585116D+00, 
     &  -1.963510026021423D+00, 
     &  -1.540619213893190D+00, 
     &  -1.220023553697935D+00, 
     &  -0.9650085667061385D+00, 
     &  -0.7549269499470514D+00, 
     &  -0.5772156649015329D+00,
     &  -0.4237549404110768D+00,
     &  -0.2890398965921883D+00,
     &  -0.1691908888667997D+00,
     &  -0.6138454458511615D-01,
     &   0.3648997397857652D-01,
     &   0.1260474527734763D+00,
     &   0.2085478748734940D+00,
     &   0.2849914332938615D+00,
     &   0.3561841611640597D+00,
     &   0.4227843350984671D+00 /

      data xvec /
     &  0.1D+00,
     &  0.2D+00,
     &  0.3D+00,
     &  0.4D+00,
     &  0.5D+00,
     &  0.6D+00,
     &  0.7D+00,
     &  0.8D+00,
     &  0.9D+00,
     &  1.0D+00,
     &  1.1D+00,
     &  1.2D+00,
     &  1.3D+00,
     &  1.4D+00,
     &  1.5D+00,
     &  1.6D+00,
     &  1.7D+00,
     &  1.8D+00,
     &  1.9D+00,
     &  2.0D+00 /

      if ( n_data .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. n_data ) then
        n_data = 0
        x = 0.0D+00
        fx = 0.0D+00
      else
        x = xvec(n_data)
        fx = fxvec(n_data)
      end if

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end