function [ f, g ] = r4_sifg ( x )

%*****************************************************************************80
%
%% R4_SIFG is a utility routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 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 F, G.
%
  persistent f1cs
  persistent f2cs
  persistent g1cs
  persistent g2cs
  persistent nf1
  persistent nf2
  persistent ng1
  persistent ng2
  persistent xbig
  persistent xbnd
  persistent xmaxf
  persistent xmaxg

  if ( isempty ( nf1 ) )

    f1cs = [ ...
      -0.1191081969051363610E+00, ...
      -0.0247823144996236248E+00, ...
       0.0011910281453357821E+00, ...
      -0.0000927027714388562E+00, ...
       0.0000093373141568271E+00, ...
      -0.0000011058287820557E+00, ...
       0.0000001464772071460E+00, ...
      -0.0000000210694496288E+00, ...
       0.0000000032293492367E+00, ...
      -0.0000000005206529618E+00, ...
       0.0000000000874878885E+00, ...
      -0.0000000000152176187E+00, ...
       0.0000000000027257192E+00, ...
      -0.0000000000005007053E+00, ...
       0.0000000000000940241E+00, ...
      -0.0000000000000180014E+00, ...
       0.0000000000000035063E+00, ...
      -0.0000000000000006935E+00, ...
       0.0000000000000001391E+00, ...
      -0.0000000000000000282E+00 ]';
     f2cs = [ ...
      -0.0348409253897013234E+00, ...
      -0.0166842205677959686E+00, ...
       0.0006752901241237738E+00, ...
      -0.0000535066622544701E+00, ...
       0.0000062693421779007E+00, ...
      -0.0000009526638801991E+00, ...
       0.0000001745629224251E+00, ...
      -0.0000000368795403065E+00, ...
       0.0000000087202677705E+00, ...
      -0.0000000022601970392E+00, ...
       0.0000000006324624977E+00, ...
      -0.0000000001888911889E+00, ...
       0.0000000000596774674E+00, ...
      -0.0000000000198044313E+00, ...
       0.0000000000068641396E+00, ...
      -0.0000000000024731020E+00, ...
       0.0000000000009226360E+00, ...
      -0.0000000000003552364E+00, ...
       0.0000000000001407606E+00, ...
      -0.0000000000000572623E+00, ...
       0.0000000000000238654E+00, ...
      -0.0000000000000101714E+00, ...
       0.0000000000000044259E+00, ...
      -0.0000000000000019634E+00, ...
       0.0000000000000008868E+00, ...
      -0.0000000000000004074E+00, ...
       0.0000000000000001901E+00, ...
      -0.0000000000000000900E+00, ...
       0.0000000000000000432E+00 ]';
    g1cs = [ ...
      -0.3040578798253495954E+00, ...
      -0.0566890984597120588E+00, ...
       0.0039046158173275644E+00, ...
      -0.0003746075959202261E+00, ...
       0.0000435431556559844E+00, ...
      -0.0000057417294453025E+00, ...
       0.0000008282552104503E+00, ...
      -0.0000001278245892595E+00, ...
       0.0000000207978352949E+00, ...
      -0.0000000035313205922E+00, ...
       0.0000000006210824236E+00, ...
      -0.0000000001125215474E+00, ...
       0.0000000000209088918E+00, ...
      -0.0000000000039715832E+00, ...
       0.0000000000007690431E+00, ...
      -0.0000000000001514697E+00, ...
       0.0000000000000302892E+00, ...
      -0.0000000000000061400E+00, ...
       0.0000000000000012601E+00, ...
      -0.0000000000000002615E+00, ...
       0.0000000000000000548E+00 ]';
    g2cs = [ ...
      -0.0967329367532432218E+00, ...
      -0.0452077907957459871E+00, ...
       0.0028190005352706523E+00, ...
      -0.0002899167740759160E+00, ...
       0.0000407444664601121E+00, ...
      -0.0000071056382192354E+00, ...
       0.0000014534723163019E+00, ...
      -0.0000003364116512503E+00, ...
       0.0000000859774367886E+00, ...
      -0.0000000238437656302E+00, ...
       0.0000000070831906340E+00, ...
      -0.0000000022318068154E+00, ...
       0.0000000007401087359E+00, ...
      -0.0000000002567171162E+00, ...
       0.0000000000926707021E+00, ...
      -0.0000000000346693311E+00, ...
       0.0000000000133950573E+00, ...
      -0.0000000000053290754E+00, ...
       0.0000000000021775312E+00, ...
      -0.0000000000009118621E+00, ...
       0.0000000000003905864E+00, ...
      -0.0000000000001708459E+00, ...
       0.0000000000000762015E+00, ...
      -0.0000000000000346151E+00, ...
       0.0000000000000159996E+00, ...
      -0.0000000000000075213E+00, ...
       0.0000000000000035970E+00, ...
      -0.0000000000000017530E+00, ...
       0.0000000000000008738E+00, ...
      -0.0000000000000004487E+00, ...
       0.0000000000000002397E+00, ...
      -0.0000000000000001347E+00, ...
       0.0000000000000000801E+00, ...
      -0.0000000000000000501E+00 ]';

    tol = 0.1 * r4_mach ( 3 );
    nf1 = r4_inits ( f1cs, 20, tol );
    nf2 = r4_inits ( f2cs, 29, tol );
    ng1 = r4_inits ( g1cs, 21, tol );
    ng2 = r4_inits ( g2cs, 34, tol );
    xbig = sqrt ( 1.0 / r4_mach ( 3 ) );
    xmaxf = exp ( min ( - log ( r4_mach ( 1 ) ), ...
      log ( r4_mach ( 2 ) ) ) - 0.01 );
    xmaxg = 1.0 / sqrt ( r4_mach ( 1 ) );
    xbnd = sqrt ( 50.0 );

  end

  if ( x < 4.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_SIFG - Fatal error!\n' );
    fprintf ( 1, '  Approximation invalid for X < 4.\n' );
    error ( 'R4_SIFG - Fatal error!' )
  end

  if ( x <= xbnd )

    f = ( 1.0 + r4_csevl ( ( 1.0 / x / x - 0.04125 ) ...
      / 0.02125, f1cs, nf1 ) ) / x;

    g = ( 1.0 + r4_csevl ( ( 1.0 / x / x - 0.04125 ) ...
      / 0.02125, g1cs, ng1 ) ) / x / x;

  elseif ( x <= xbig )

    f = ( 1.0 + r4_csevl ( 100.0 / x / x - 1.0, f2cs, nf2 ) ) / x;
    g = ( 1.0 + r4_csevl ( 100.0 / x / x - 1.0, g2cs, ng2) ) / x / x;

  else

    if ( x < xmaxf )
      f = 1.0 / x;
    else
      f = 0.0;
    end

    if ( x < xmaxg )
      g = 1.0 / x / x;
    else
      g = 0.0;
    end

  end

  return
end