function value = error_fc ( x )

%*****************************************************************************80
%
%% ERROR_FC computes the complementary error function.
%
%  Discussion:
%
%    This function was renamed "ERROR_FC" from "ERFC", to avoid a conflict
%    with the name of a corresponding routine often, but not always,
%    supplied as part of the math support library.
%
%    The definition of the complementary error function is:
%
%      ERFC(X) = 1 - ERF(X)
%
%    where ERF(X) is the error function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    FORTRAN77 original version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    David Kahaner, Cleve Moler, Steven Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1989,
%    ISBN: 0-13-627258-4,
%    LC: TA345.K34.
%
%  Parameters:
%
%    Input, real X, the argument of the error function.
%
%    Output, real VALUE, the value of the complementary 
%    error function at X.
%
  persistent nterc2
  persistent nterf
  persistent nterfc
  persistent sqeps
  persistent xmax
  persistent xsml

  erfcs = [ ...
   -0.049046121234691808, ...
   -0.14226120510371364,  ...
    0.010035582187599796, ...
   -0.000576876469976748, ...
    0.000027419931252196, ...
   -0.000001104317550734, ...
    0.000000038488755420, ...
   -0.000000001180858253, ...
    0.000000000032334215, ...
   -0.000000000000799101, ...
    0.000000000000017990, ...
   -0.000000000000000371, ...
    0.000000000000000007 ];

  erfccs = [ ...
    0.0715179310202925, ...
   -0.026532434337606719, ...
    0.001711153977920853, ...
   -0.000163751663458512, ...
    0.000019871293500549, ...
   -0.000002843712412769, ...
    0.000000460616130901, ...
   -0.000000082277530261, ...
    0.000000015921418724, ...
   -0.000000003295071356, ...
    0.000000000722343973, ...
   -0.000000000166485584, ...
    0.000000000040103931, ...
   -0.000000000010048164, ...
    0.000000000002608272, ...
   -0.000000000000699105, ...
    0.000000000000192946, ...
   -0.000000000000054704, ...
    0.000000000000015901, ...
   -0.000000000000004729, ...
    0.000000000000001432, ...
   -0.000000000000000439, ...
    0.000000000000000138, ...
   -0.000000000000000048 ];

  erc2cs = [ ...
   -0.069601346602309501, ...
   -0.041101339362620893, ...
    0.003914495866689626, ...
   -0.000490639565054897, ...
    0.000071574790013770, ...
   -0.000011530716341312, ...
    0.000001994670590201, ...
   -0.000000364266647159, ...
    0.000000069443726100, ...
   -0.000000013712209021, ...
    0.000000002788389661, ...
   -0.000000000581416472, ...
    0.000000000123892049, ...
   -0.000000000026906391, ...
    0.000000000005942614, ...
   -0.000000000001332386, ...
    0.000000000000302804, ...
   -0.000000000000069666, ...
    0.000000000000016208, ...
   -0.000000000000003809, ...
    0.000000000000000904, ...
   -0.000000000000000216, ...
    0.000000000000000052 ];

  sqrtpi = 1.7724538509055160;

  if ( size ( nterf ) == 0 )

    eta = 0.1 * r8_epsilon ( );
    nterf = inits ( erfcs, 13, eta );
    nterfc = inits ( erfccs, 24, eta );
    nterc2 = inits ( erc2cs, 23, eta );

    xsml = -sqrt ( - log ( sqrtpi * r8_epsilon ( ) ) );
    xmax = sqrt ( - log ( sqrtpi * r8_tiny ( ) ) );
    xmax = xmax - 0.5 * log ( xmax ) / xmax - 0.01;
    sqeps = sqrt ( 2.0 * r8_epsilon ( ) );

  end

  if ( x <= xsml )
    value = 2.0;
    return
  end
%
%  X so big that ERFC will underflow.
%
  if ( xmax < x )
    value = 0.0;
    return
  end

  y = abs ( x );
%
%  erfc(x) = 1.0 - erf(x) for -1 <= x <= 1.
%
  if ( y <= 1.0 )

    if ( y < sqeps )
      value = 1.0 - 2.0 * x / sqrtpi;
    elseif ( sqeps <= y )
      value = 1.0 - x * ( 1.0 + csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
    end

    return

  end
%
%  For 1 < |x| <= xmax, erfc(x) = 1.0 - erf(x)
%
  y = y * y;

  if ( y <= 4.0 )
    value = exp ( -y ) / abs ( x ) * ( 0.5 ...
      + csevl ( ( 8.0 / y - 5.0 ) / 3.0, erc2cs, nterc2 ) );
  else
    value = exp ( -y ) / abs ( x ) * ( 0.5 ...
      + csevl ( 8.0 / y - 1.0, erfccs, nterfc ) );
  end

  if ( x < 0.0 )
    value = 2.0 - value;
  end

  return
end