# include <cmath>
# include <cstdlib>

using namespace std;

# include "sinc.hpp"

//****************************************************************************80

void cisi ( double x, double &ci, double &si )

//****************************************************************************80
//
//  Purpose:
//
//    cisi() computes cosine Ci(x) and sine integrals 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:
//
//    27 March 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 of Ci(x) and Si(x).
//
//  Output:
//
//    double &ci, &si: the values of Ci(x) and Si(x).
//
{
  double *bj;
  double el;
  double eps;
  int k;
  int m;
  double p2;
  double x2;
  double xa;
  double xa0;
  double xa1;
  double xcs;
  double xf;
  double xg;
  double xg1;
  double xg2;
  double xr;
  double xs;
  double xss;

  p2 = 1.570796326794897;
  el = 0.5772156649015329;
  eps = 1.0E-15;
  x2 = x * x;

  if ( x == 0.0 )
  {
    ci = -1.0e+300;
    si = 0.0;
  }
  else if ( x <= 16.0 )
  {
    xr = -0.25 * x2;
    ci = el + log ( x ) + xr;
    for ( k = 2; k <= 40; k++ )
    {
      xr = -0.5 * xr * ( k - 1 ) / ( k * k * ( 2 * k - 1 ) ) * x2;
      ci = ci + xr;
      if ( fabs ( xr ) < fabs ( ci ) * eps )
      {
        break;
      }
    }

    xr = x;
    si = x;
    for ( k = 1; k <= 40; k++ )
    {
      xr = -0.5 * xr * ( 2 * k - 1 ) / k / ( 4 * k * k + 4 * k + 1 ) * x2;
      si = si + xr;
      if ( fabs ( xr ) < fabs ( si ) * eps )
      {
        return;
      }
    }
  }
  else if ( x <= 32.0 )
  {
    m = ( int ) ( 47.2 + 0.82 * x );
    bj = new double[m];
    xa1 = 0.0;
    xa0 = 1.0E-100;
    for ( k = m; 1 <= k; k-- )
    {
      xa = 4.0 * k * xa0 / x - 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 ) ) * x;
      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 ) ) * x;
      xg2 = xg2 + bj[k-1] * xr;
    }

    xcs = cos ( x / 2.0 );
    xss = sin ( x / 2.0 );
    ci = el + log ( x ) - x * xss * xg1 + 2.0 * xcs * xg2 - 2.0 * xcs * xcs;
    si = x * xcs * xg1 + 2.0 * xss * xg2 - sin ( x );

    delete [] 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 / x;
    xg = xr;
    for ( k = 1; k <= 8; k++ )
    {
      xr = -2.0 * xr * ( 2 * k + 1 ) * k / x2;
      xg = xg + xr;
    }
    ci = xf * sin ( x ) / x - xg * cos ( x ) / x;
    si = p2 - xf * cos ( x ) / x - xg * sin ( x ) / x;
  }

  return;
}
//****************************************************************************80

double sincn_antideriv ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincn_antideriv(): antiderivative of the normalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincn_antideriv: the value of the antiderivative.
//
{
  double ci;
  static double r8_pi = 3.141592653589793;
  double si;
  double value;

  cisi ( r8_pi * x, ci, si );
  value = si / r8_pi;

  return value;
}
//****************************************************************************80

double sincn_deriv2 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincn_deriv2(): second derivative of the normalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincn_deriv2: the value of the second derivative.
//
{
  static double r8_pi = 3.141592653589793;
  double value;

  if ( fabs ( x ) <= 1.0e-10 )
  {
    value = - pow ( r8_pi, 2 ) / 3.0;
  }
  else
  {
    value = 
      ( ( 2.0 - pow ( r8_pi * x, 2 ) ) * sin ( r8_pi * x ) 
        - 2.0 *       r8_pi * x        * cos ( r8_pi * x ) ) 
      / ( r8_pi * pow ( x, 3 ) );
  }

  return value;
}
//****************************************************************************80

double sincn_deriv ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincn_deriv() evaluates the derivative of the normalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincn_deriv: the value of the derivative.
//
{
  static double r8_pi = 3.141592653589793;
  double value;

  if ( x == 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = ( r8_pi * x * cos ( r8_pi * x )
                        - sin ( r8_pi * x ) ) / ( r8_pi * pow ( x, 2 ) );
  }

  return value;
}
//****************************************************************************80

double sincn_fun ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincn_fun() evaluates the normalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincn_fun: the value of the function.
//
{
  static double r8_pi = 3.141592653589793;
  double value;

  if ( x == 0.0 )
  {
    value = 1.0;
  }
  else
  {
    value = sin ( r8_pi * x ) / ( r8_pi * x );
  }

  return value;
}
//****************************************************************************80

double sincu_antideriv ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincu_antideriv(): antiderivative of the unnormalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincu_antideriv: the value of the antiderivative.
//
{
  double ci;
  double si;
  double value;

  cisi ( x, ci, si );
  value = si;

  return value;
}
//****************************************************************************80

double sincu_deriv2 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincu_deriv2(): second derivative of the unnormalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincu_deriv2: the value of the second derivative.
//
{
  double value;

  if ( x == 0.0 )
  {
    value = - 1.0 / 3.0;
  }
  else
  {
    value = ( ( 2.0 - pow ( x, 2 ) ) * sin ( x )
              - 2.0 *       x        * cos ( x ) ) / pow ( x, 3 );
  }

  return value;
}
//****************************************************************************80

double sincu_deriv ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincu_deriv() evaluates the derivative of the unnormalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincu_deriv: the value of the derivative.
//
//
{
  double value;

  if ( x == 0.0 )
  {
    value = 0.0;
  }
  else
  {
    value = ( x * cos ( x ) - sin ( x ) ) / pow ( x, 2 );
  }

  return value;
}
//****************************************************************************80

double sincu_fun ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    sincu_fun() evaluates the unnormalized sinc() function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real x: the argument of the function.
//
//  Output:
//
//    real sincu_fun: the value of the function.
//
{
  double value;

  if ( x == 0.0 )
  {
    value = 1.0;
  }
  else
  {
    value = sin ( x ) / x;
  }

  return value;
}