# include <math.h>
# include <stdlib.h>

# include "sine_integral.h"

/******************************************************************************/

double si ( double x )

/******************************************************************************/
/*
  Purpose:

    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:

    10 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:

    double x: the argument.

  Output:

    double *si: the value of Si(x).
*/
{
  double *bj;
  double eps;
  int k;
  int m;
  double p2;
  double value;
  double x2;
  double xa;
  double xa0;
  double xa1;
  double xabs;
  double xcs;
  double xf;
  double xg;
  double xg1;
  double xg2;
  double xr;
  double xs;
  double xsign;
  double xss;

  p2 = 1.570796326794897;
  eps = 1.0E-15;
  x2 = x * x;
  xabs = fabs ( x );

  if ( x < 0.0 )
  {
    xsign = -1.0;
  }
  else
  {
    xsign = 1.0;
  }

  if ( xabs == 0.0 )
  {
    value = 0.0;
  }
  else if ( xabs <= 16.0 )
  {
    xr = xabs;
    value = xabs;
    for ( k = 1; k <= 40; k++ )
    {
      xr = - 0.5 * xr * ( 2 * k - 1 ) / k / ( 4 * k * k + 4 * k + 1 ) * x2;
      value = value + xr;
      if ( fabs ( xr ) < fabs ( value ) * eps )
      {
        value = xsign * value;
        return value;
      }
    }
  }
  else if ( xabs <= 32.0 )
  {
    m = ( int ) ( 47.2 + 0.82 * xabs );
    bj = ( double * ) malloc ( m * sizeof ( double ) );
    xa1 = 0.0;
    xa0 = 1.0E-100;
    for ( k = m; 1 <= k; k-- )
    {
      xa = 4.0 * k * xa0 / xabs - xa1;
      bj[k-1] = xa;
      xa1 = xa0;
      xa0 = xa;
    }
    xs = bj[0];
    for ( k = 3; k <= m; k = k + 2 )
    {
      xs = xs + 2.0 * bj[k-1];
    }
    bj[0] = bj[0] / xs;
    for ( k = 2; k <= m; k++ )
    {
      bj[k-1] = bj[k-1] / xs;
    }
    xr = 1.0;
    xg1 = bj[0];
    for ( k = 2; k <= m; k++ )
    {
      xr = 0.25 * xr * pow ( 2.0 * k - 3.0, 2 )
        / ( ( k - 1.0 ) * pow ( 2.0 * k - 1.0, 2 ) ) * xabs;
      xg1 = xg1 + bj[k-1] * xr;
    }

    xr = 1.0;
    xg2 = bj[0];
    for ( k = 2; k <= m; k++ )
    {
      xr = 0.25 * xr * pow ( 2.0 * k - 5.0, 2 )
        / ( ( k - 1.0 ) * pow ( 2.0 * k - 3.0, 2 ) ) * xabs;
      xg2 = xg2 + bj[k-1] * xr;
    }

    xcs = cos ( xabs / 2.0 );
    xss = sin ( xabs / 2.0 );
    value = xsign * ( xabs * xcs * xg1 + 2.0 * xss * xg2 - sin ( xabs ) );

    free ( bj );
  }
  else
  {
    xr = 1.0;
    xf = 1.0;
    for ( k = 1; k <= 9; k++ )
    {
      xr = - 2.0 * xr * k * ( 2 * k - 1 ) / x2;
      xf = xf + xr;
    }
    xr = 1.0 / xabs;
    xg = xr;
    for ( k = 1; k <= 8; k++ )
    {
      xr = - 2.0 * xr * ( 2 * k + 1 ) * k / x2;
      xg = xg + xr;
    }
    value = xsign * ( p2 - xf * cos ( xabs ) / xabs - xg * sin ( xabs ) / xabs );
  }

  return value;
}