function value = r8_spence ( x )

%*****************************************************************************80
%
%% R8_SPENCE evaluates a form of Spence's function for an R8 argument.
%
%  Discussion:
%
%    This function evaluates a form of Spence's function defined by
%
%      f(x) = Integral ( 0 <= y <= x ) - log ( abs ( 1 - y ) ) / y dy
%
%  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:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions, page 1004,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    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.
%
%    K Mitchell,
%    Tables of the function Integral ( 0 < y < x ) - log | 1 - y | dy / y
%    with an account of some properties of this and related functions,
%    Philosophical Magazine,
%    Volume 40, pages 351-368, 1949.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, Spence's function evaluated at X.
%
  persistent nspenc
  persistent pi26
  persistent spencs
  persistent xbig

  pi26 = +1.644934066848226436472415166646025189219;

  if ( isempty ( nspenc ) )

    spencs = [ ...
      +0.1527365598892405872946684910028, ...
      +0.8169658058051014403501838185271E-01, ...
      +0.5814157140778730872977350641182E-02, ...
      +0.5371619814541527542247889005319E-03, ...
      +0.5724704675185826233210603054782E-04, ...
      +0.6674546121649336343607835438589E-05, ...
      +0.8276467339715676981584391689011E-06, ...
      +0.1073315673030678951270005873354E-06, ...
      +0.1440077294303239402334590331513E-07, ...
      +0.1984442029965906367898877139608E-08, ...
      +0.2794005822163638720201994821615E-09, ...
      +0.4003991310883311823072580445908E-10, ...
      +0.5823462892044638471368135835757E-11, ...
      +0.8576708692638689278097914771224E-12, ...
      +0.1276862586280193045989483033433E-12, ...
      +0.1918826209042517081162380416062E-13, ...
      +0.2907319206977138177795799719673E-14, ...
      +0.4437112685276780462557473641745E-15, ...
      +0.6815727787414599527867359135607E-16, ...
      +0.1053017386015574429547019416644E-16, ...
      +0.1635389806752377100051821734570E-17, ...
      +0.2551852874940463932310901642581E-18, ...
      +0.3999020621999360112770470379519E-19, ...
      +0.6291501645216811876514149171199E-20, ...
      +0.9933827435675677643803887752533E-21, ...
      +0.1573679570749964816721763805866E-21, ...
      +0.2500595316849476129369270954666E-22, ...
      +0.3984740918383811139210663253333E-23, ...
      +0.6366473210082843892691326293333E-24, ...
      +0.1019674287239678367077061973333E-24, ...
      +0.1636881058913518841111074133333E-25, ...
      +0.2633310439417650117345279999999E-26, ...
      +0.4244811560123976817224362666666E-27, ...
      +0.6855411983680052916824746666666E-28, ...
      +0.1109122433438056434018986666666E-28, ...
      +0.1797431304999891457365333333333E-29, ...
      +0.2917505845976095173290666666666E-30, ...
      +0.4742646808928671061333333333333E-31 ]';

    nspenc = r8_inits ( spencs, 38, 0.1 * r8_mach ( 3 ) );
    xbig = 1.0 / r8_mach ( 3 );
  end

  if ( x <= - xbig )

    aln = log ( 1.0 - x );
    value = - pi26  - 0.5 * aln * ( 2.0 * log ( - x ) - aln );

  elseif ( x <= - 1.0 )

    aln = log ( 1.0 - x );

    value = - pi26 - 0.5 * aln * ( 2.0 ...
      * log ( - x ) - aln ) + ( 1.0 + r8_csevl ( ...
      4.0 / ( 1.0 - x ) - 1.0, spencs, nspenc ) ) / ( 1.0 - x );

  elseif ( x <= 0.0 )

    value = - 0.5 * log ( 1.0 - x ) ...
      * log ( 1.0 - x ) - x * ( 1.0 + r8_csevl ( ...
      4.0 * x / ( x - 1.0 ) - 1.0, spencs, nspenc ) ) ...
      / ( x - 1.0 );

  elseif ( x <= 0.5 )

    value = x * ( 1.0 + r8_csevl ( 4.0 * x - 1.0, spencs, nspenc ) );

  elseif ( x < 1.0 )

    value = pi26 - log ( x ) * log ( 1.0 - x ) ...
      - ( 1.0 - x ) * ( 1.0 + r8_csevl ( 4.0 ...
      * ( 1.0 - x ) - 1.0, spencs, nspenc ) );

  elseif ( x == 1.0 )

    value = pi26;

  elseif ( x <= 2.0 )

    value = pi26 - 0.5 * log ( x ) ...
      * log ( ( x - 1.0 ) * ( x - 1.0 ) / x ) ...
      + ( x - 1.0 ) * ( 1.0 + r8_csevl ( 4.0 ...
      * ( x - 1.0 ) / x - 1.0, spencs, nspenc ) ) / x;

  elseif ( x < xbig )

    value = 2.0 * pi26 - 0.5 * log ( x ) * log ( x ) ...
      - ( 1.0 + r8_csevl ( 4.0 / x - 1.0, spencs, ...
      nspenc ) ) / x;

  else

    value = 2.0 * pi26 - 0.5 * log ( x ) * log ( x );

  end

  return
end