function value = r8_psi ( x )

%*****************************************************************************80
%
%% R8_PSI evaluates the psi function of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the psi function of X.
%
  persistent apsics
  persistent dxrel
  persistent ntapsi
  persistent ntpsi
  persistent psics
  persistent xbig

  if ( isempty ( ntpsi ) )

    apsics = [ ...
      -0.832710791069290760174456932269E-03, ...
      -0.416251842192739352821627121990E-03, ...
      +0.103431560978741291174463193961E-06, ...
      -0.121468184135904152987299556365E-09, ...
      +0.311369431998356155521240278178E-12, ...
      -0.136461337193177041776516100945E-14, ...
      +0.902051751315416565130837974000E-17, ...
      -0.831542997421591464829933635466E-19, ...
      +0.101224257073907254188479482666E-20, ...
      -0.156270249435622507620478933333E-22, ...
      +0.296542716808903896133226666666E-24, ...
      -0.674686886765702163741866666666E-26, ...
      +0.180345311697189904213333333333E-27, ...
      -0.556901618245983607466666666666E-29, ...
      +0.195867922607736251733333333333E-30, ...
      -0.775195892523335680000000000000E-32 ]';
    psics = [ ...
      -0.38057080835217921520437677667039E-01, ...
      +0.49141539302938712748204699654277, ...
      -0.56815747821244730242892064734081E-01, ...
      +0.83578212259143131362775650747862E-02, ...
      -0.13332328579943425998079274172393E-02, ...
      +0.22031328706930824892872397979521E-03, ...
      -0.37040238178456883592889086949229E-04, ...
      +0.62837936548549898933651418717690E-05, ...
      -0.10712639085061849855283541747074E-05, ...
      +0.18312839465484165805731589810378E-06, ...
      -0.31353509361808509869005779796885E-07, ...
      +0.53728087762007766260471919143615E-08, ...
      -0.92116814159784275717880632624730E-09, ...
      +0.15798126521481822782252884032823E-09, ...
      -0.27098646132380443065440589409707E-10, ...
      +0.46487228599096834872947319529549E-11, ...
      -0.79752725638303689726504797772737E-12, ...
      +0.13682723857476992249251053892838E-12, ...
      -0.23475156060658972717320677980719E-13, ...
      +0.40276307155603541107907925006281E-14, ...
      -0.69102518531179037846547422974771E-15, ...
      +0.11856047138863349552929139525768E-15, ...
      -0.20341689616261559308154210484223E-16, ...
      +0.34900749686463043850374232932351E-17, ...
      -0.59880146934976711003011081393493E-18, ...
      +0.10273801628080588258398005712213E-18, ...
      -0.17627049424561071368359260105386E-19, ...
      +0.30243228018156920457454035490133E-20, ...
      -0.51889168302092313774286088874666E-21, ...
      +0.89027730345845713905005887487999E-22, ...
      -0.15274742899426728392894971904000E-22, ...
      +0.26207314798962083136358318079999E-23, ...
      -0.44964642738220696772598388053333E-24, ...
      +0.77147129596345107028919364266666E-25, ...
      -0.13236354761887702968102638933333E-25, ...
      +0.22709994362408300091277311999999E-26, ...
      -0.38964190215374115954491391999999E-27, ...
      +0.66851981388855302310679893333333E-28, ...
      -0.11469986654920864872529919999999E-28, ...
      +0.19679385886541405920515413333333E-29, ...
      -0.33764488189750979801907200000000E-30, ...
      +0.57930703193214159246677333333333E-31 ]';

    ntpsi = r8_inits ( psics, 42, 0.1 * r8_mach ( 3 ) );
    ntapsi = r8_inits ( apsics, 16, 0.1 * r8_mach ( 3 ) );
    xbig = 1.0 / sqrt ( r8_mach ( 3 ) );
    dxrel = sqrt ( r8_mach ( 4 ) );

  end

  y = abs ( x );

  if ( y < 10.0 )

    n = r8_aint ( x );
    if ( x < 0.0 )
      n = n - 1;
    end
    y = x - n;
    n = n - 1;
    value = r8_csevl ( 2.0 * y - 1.0, psics, ntpsi );

    if ( n == 0 )

      return

    elseif ( n < 0 )

      n = - n;

      if ( x == 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8_PSI - Fatal error!\n' );
        fprintf ( 1, '  X is zero.\n' );
        error ( 'R8_PSI - Fatal error!' )
      end

      if ( x < 0.0 && x + n - 2 == 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8_PSI - Fatal error!\n' );
        fprintf ( 1, '  X is a negative integer.\n' );
        error ( 'R8_PSI - Fatal error!' )
      end

      if ( x < - 0.5 && ...
        abs ( ( x - r8_aint ( x - 0.5 ) ) / x ) < dxrel )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R8_PSI - Warning!\n' );
        fprintf ( 1, '  Answer is less than half precision\n' );
        fprintf ( 1, '  because X is near a negative integer.\n' );
      end

      for i = 1 : n
        value = value - 1.0 / ( x + i - 1 );
      end

    elseif ( 0 < n )

      for i = 1 : n
        value = value + 1.0 / ( y + i );
      end

    end

  else

    if ( y < xbig )
      aux = r8_csevl ( 8.0 / y / y - 1.0, apsics, ntapsi );
    else
      aux = 0.0;
    end

    if ( x < 0.0 )

      value = log ( abs ( x ) ) - 0.5 / x + aux - pi * r8_cot ( pi * x );

    elseif ( 0.0 < x )

      value = log ( x ) - 0.5 / x + aux;

    end

  end

  return
end