function si ( x )

!*****************************************************************************80
!
!! si() computes the sine integral Si(x).
!
!  Licensing:
!
!    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However,
!    they give permission to incorporate this routine into a user program
!    provided that the copyright is acknowledged.
!
!  Modified:
!
!    09 August 2025
!
!  Author:
!
!    Original Fortran77 version by Shanjie Zhang, Jianming Jin.
!    This version by John Burkardt.
!
!  Reference:
!
!    Shanjie Zhang, Jianming Jin,
!    Computation of Special Functions,
!    Wiley, 1996,
!    ISBN: 0-471-11963-6,
!    LC: QA351.C45.
!
!  Input:
!
!    real ( kind = rk8 ) x: the argument.
!
!  Output:
!
!    real ( kind = rk8 ) si: the value of Si(x).
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) bj(101)
  real ( kind = rk8 ) el
  real ( kind = rk8 ) eps
  integer k
  integer m
  real ( kind = rk8 ) p2
  real ( kind = rk8 ) si
  real ( kind = rk8 ) x
  real ( kind = rk8 ) x2
  real ( kind = rk8 ) xa
  real ( kind = rk8 ) xa0
  real ( kind = rk8 ) xa1
  real ( kind = rk8 ) xabs
  real ( kind = rk8 ) xcs
  real ( kind = rk8 ) xf
  real ( kind = rk8 ) xg
  real ( kind = rk8 ) xg1
  real ( kind = rk8 ) xg2
  real ( kind = rk8 ) xr
  real ( kind = rk8 ) xs
  real ( kind = rk8 ) xsign
  real ( kind = rk8 ) xss

  p2 = 1.570796326794897D+00
  el = 0.5772156649015329D+00
  eps = 1.0D-15
  x2 = x * x
  xabs = abs ( x )
  if ( x < 0.0D+00 ) then
    xsign = -1.0D+00
  else
    xsign = +1.0D+00
  end if

  if ( xabs == 0.0D+00 ) then

    si = 0.0D+00

  else if ( xabs <= 16.0D+00 ) then

    xr = xabs
    si = xabs
    do k = 1, 40
      xr = - 0.5D+00 * xr * ( 2 * k - 1 ) / k / ( 4 * k * k + 4 * k + 1 ) * x2
      si = si + xr
      if ( abs ( xr ) < abs ( si ) * eps ) then
        si = xsign * si
        return
      end if
    end do

  else if ( xabs <= 32.0D+00 ) then

    m = int ( 47.2D+00 + 0.82D+00 * xabs )
    xa1 = 0.0D+00
    xa0 = 1.0D-100
    do k = m, 1, -1
      xa = 4.0D+00 * k * xa0 / xabs - xa1
      bj(k) = xa
      xa1 = xa0
      xa0 = xa
    end do
    xs = bj(1)
    do k = 3, m, 2
      xs = xs + 2.0D+00 * bj(k)
    end do
    bj(1) = bj(1) / xs
    do k = 2, m
      bj(k) = bj(k) / xs
    end do
    xr = 1.0D+00
    xg1 = bj(1)
    do k = 2, m
      xr = 0.25D+00 * xr * ( 2.0D+00 * k - 3.0D+00 )**2 &
        / ( ( k - 1.0D+00 ) * ( 2.0D+00 * k - 1.0D+00 )**2 ) * xabs
      xg1 = xg1 + bj(k) * xr
    end do

    xr = 1.0D+00
    xg2 = bj(1)
    do k = 2, m
      xr = 0.25D+00 * xr * ( 2.0D+00 * k - 5.0D+00 )**2 &
        / ( ( k - 1.0D+00 ) * ( 2.0D+00 * k - 3.0D+00 )**2 ) * xabs
      xg2 = xg2 + bj(k) * xr
    end do

    xcs = cos ( xabs / 2.0D+00 )
    xss = sin ( xabs / 2.0D+00 )
    si = xsign * ( xabs * xcs * xg1 + 2.0 * xss * xg2 - sin ( xabs ) )

  else

    xr = 1.0D+00
    xf = 1.0D+00
    do k = 1, 9
      xr = - 2.0D+00 * xr * k * ( 2 * k - 1 ) / x2
      xf = xf + xr
    end do
    xr = 1.0D+00 / xabs
    xg = xr
    do k = 1, 8
      xr = - 2.0D+00 * xr * ( 2 * k + 1 ) * k / x2
      xg = xg + xr
    end do
    si = xsign * ( p2 - xf * cos ( xabs ) / xabs - xg * sin ( xabs ) / xabs )

  end if

  return
end