# include <cfloat>
# include <cmath>
# include <cstdlib>

using namespace std;

# include "sncndn.hpp"

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

double jacobi_cn ( double u, double m )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_cn() evaluates the Jacobi elliptic function CN(U,M).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2025
//
//  Author:
//
//    Original ALGOL version by Roland Bulirsch.
//    This version by John Burkardt
//
//  Reference:
//
//    Roland Bulirsch,
//    Numerical calculation of elliptic integrals and elliptic functions,
//    Numerische Mathematik,
//    Volume 7, Number 1, 1965, pages 78-90.
//
//  Input:
//
//    double U, M, the arguments.
//
//  Output:
//
//    double JACOBI_CN, the function value.
//
{
  double cn;
  double dn;
  double sn;

  sncndn ( u, m, sn, cn, dn );

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

double jacobi_dn ( double u, double m )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_dn() evaluates the Jacobi elliptic function DN(U,M).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2025
//
//  Author:
//
//    Original ALGOL version by Roland Bulirsch.
//    This version by John Burkardt
//
//  Reference:
//
//    Roland Bulirsch,
//    Numerical calculation of elliptic integrals and elliptic functions,
//    Numerische Mathematik,
//    Volume 7, Number 1, 1965, pages 78-90.
//
//  Input:
//
//    double U, M, the arguments.
//
//  Output:
//
//    double JACOBI_DN, the function value.
//
{
  double cn;
  double dn;
  double sn;

  sncndn ( u, m, sn, cn, dn );

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

double jacobi_sn ( double u, double m )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_sn() evaluates the Jacobi elliptic function SN(U,M).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2025
//
//  Author:
//
//    Original ALGOL version by Roland Bulirsch.
//    This version by John Burkardt
//
//  Reference:
//
//    Roland Bulirsch,
//    Numerical calculation of elliptic integrals and elliptic functions,
//    Numerische Mathematik,
//    Volume 7, Number 1, 1965, pages 78-90.
//
//  Input:
//
//    double U, M, the arguments.
//
//  Output:
//
//    double JACOBI_SN, the function value.
//
{
  double cn;
  double dn;
  double sn;

  sncndn ( u, m, sn, cn, dn );

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

void sncndn ( double u, double m, double &sn, double &cn, double &dn )

//****************************************************************************80
//
//  Purpose:
//
//    sncndn() evaluates Jacobi elliptic functions SN, CN, and DN.
//
//  Discussion:
//
//    Evaluation
//      a = - cn / sn;
//    corrected to
//      a = cn / sn;
//    after comparing to reference, 26 June 2025.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2025
//
//  Author:
//
//    Original ALGOL version by Roland Bulirsch.
//    This version by John Burkardt
//
//  Reference:
//
//    Roland Bulirsch,
//    Numerical calculation of elliptic integrals and elliptic functions,
//    Numerische Mathematik,
//    Volume 7, Number 1, 1965, pages 78-90.
//
//  Input:
//
//    double U, M, the arguments.
//
//  Output:
//
//    double &SN, &CN, &DN, the value of the Jacobi
//    elliptic functions sn(u,m), cn(u,m), and dn(u,m).
//
{
# define MAXIT 25

  double a;
  double b;
  double c;
  double ca;
  double d;
  double m_array[MAXIT];
  double n_array[MAXIT];
  int i;
  int l;
  double m_comp;
  double u_copy;

  m_comp = 1.0 - m;
  u_copy = u;

  if ( m_comp == 0.0 )
  {
    cn = 1.0 / cosh ( u_copy );
    dn = cn;
    sn = tanh ( u_copy );
    return;
  }

  if ( 1.0 < m )
  {
    d = 1.0 - m_comp;
    m_comp = - m_comp / d;
    d = sqrt ( d );
    u_copy = d * u_copy;
  }

  ca = sqrt ( DBL_EPSILON );

  a = 1.0;
  dn = 1.0;
  l = MAXIT;

  for ( i = 1; i <= MAXIT; i++ )
  {
    m_array[i-1] = a;
    m_comp = sqrt ( m_comp );
    n_array[i-1] = m_comp;
    c = 0.5 * ( a + m_comp );

    if ( fabs ( a - m_comp ) <= ca * a )
    {
      l = i;
      break;
    }

    m_comp = a * m_comp;
    a = c;
  }

  u_copy = c * u_copy;
  sn = sin ( u_copy );
  cn = cos ( u_copy );

  if ( sn != 0.0 )
  {
    a = cn / sn;
    c = a * c;

    for ( i = l; 1 <= i; i-- )
    {
      b = m_array[i-1];
      a = c * a;
      c = dn * c;
      dn = ( n_array[i-1] + a ) / ( b + a );
      a = c / b;
    }

    a = 1.0 / sqrt ( c * c + 1.0 );

    if ( sn < 0.0 )
    {
      sn = - a;
    }
    else
    {
      sn = a;
    }
    cn = c * sn;
  }

  if ( 1.0 < m )
  {
    a = dn;
    dn = cn;
    cn = a;
    sn = sn / d;
  }

  return;
# undef MAXIT
}
//****************************************************************************80

void jacobi_cn_values ( int &n_data, double &u, double &a, double &k, 
  double &m, double &fx )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_cn_values() returns values of the Jacobi elliptic function CN(U,M).
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      JacobiCN[ u, m ]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 November 2020
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Input:
//
//    int &N_DATA.  The user sets N_DATA to 0 before the first call.  
//
//  Output:
//
//    int &N_DATA.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    double &U the argument of the function.
//
//    double &A, &K, &M, the parameters of the function.
//
//    double &FX, the value of the function.
//
{
# define N_MAX 20

  static double m_vec[N_MAX] = {
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00 };

  static double fx_vec[N_MAX] = {
      0.9950041652780258E+00,
      0.9800665778412416E+00,
      0.8775825618903727E+00,
      0.5403023058681397E+00,
     -0.4161468365471424E+00,
      0.9950124626090582E+00,
      0.9801976276784098E+00,
      0.8822663948904403E+00,
      0.5959765676721407E+00,
     -0.1031836155277618E+00,
      0.9950207489532265E+00,
      0.9803279976447253E+00,
      0.8868188839700739E+00,
      0.6480542736638854E+00,
      0.2658022288340797E+00,
      0.3661899347368653E-01,
      0.9803279976447253E+00,
      0.8868188839700739E+00,
      0.6480542736638854E+00,
      0.2658022288340797E+00 };

  static double u_vec[N_MAX] = {
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      4.0E+00,
     -0.2E+00,
     -0.5E+00,
     -1.0E+00,
     -2.0E+00 };

  if ( n_data < 0 )
  {
    n_data = 0;
  }

  n_data = n_data + 1;

  if ( N_MAX < n_data )
  {
    n_data = 0;
    a = 0.0;
    k = 0.0;
    m = 0.0;
    u = 0.0;
    fx = 0.0;
  }
  else
  {
    m = m_vec[n_data-1];
    k = sqrt ( m );
    a = asin ( k );
    u = u_vec[n_data-1];
    fx = fx_vec[n_data-1];
  }

  return;
# undef N_MAX
}
//****************************************************************************80

void jacobi_dn_values ( int &n_data, double &u, double &a, double &k, 
  double &m, double &fx )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_dn_values() returns values of the Jacobi elliptic function DN(U,M).
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      JacobiDN[ u, m ]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 November 2020
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Input:
//
//    int &N_DATA.  The user sets N_DATA to 0 before the first call.  
//
//  Output:
//
//    int &N_DATA.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    double &U the argument of the function.
//
//    double &A, &K, &M, the parameters of the function.
//
//    double &FX, the value of the function.
//
{
# define N_MAX 20

  static double m_vec[N_MAX] = {
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00 };

  static double fx_vec[N_MAX] = {
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.1000000000000000E+01,
     0.9975093485144243E+00,
     0.9901483195224800E+00,
     0.9429724257773857E+00,
     0.8231610016315963E+00,
     0.7108610477840873E+00,
     0.9950207489532265E+00,
     0.9803279976447253E+00,
     0.8868188839700739E+00,
     0.6480542736638854E+00,
     0.2658022288340797E+00,
     0.3661899347368653E-01,
     0.9803279976447253E+00,
     0.8868188839700739E+00,
     0.6480542736638854E+00,
     0.2658022288340797E+00  };

  static double u_vec[N_MAX] = {
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      4.0E+00,
     -0.2E+00,
     -0.5E+00,
     -1.0E+00,
     -2.0E+00 };

  if ( n_data < 0 )
  {
    n_data = 0;
  }

  n_data = n_data + 1;

  if ( N_MAX < n_data )
  {
    n_data = 0;
    a = 0.0;
    k = 0.0;
    m = 0.0;
    u = 0.0;
    fx = 0.0;
  }
  else
  {
    m = m_vec[n_data-1];
    k = sqrt ( m );
    a = asin ( k );
    u = u_vec[n_data-1];
    fx = fx_vec[n_data-1];
  }

  return;
# undef N_MAX
}
//****************************************************************************80

void jacobi_sn_values ( int &n_data, double &u, double &a, double &k, 
  double &m, double &fx )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_sn_values() returns values of the Jacobi elliptic function SN(U,M).
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      JacobiSN[ u, m ]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 November 2020
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Input:
//
//    int &N_DATA.  The user sets N_DATA to 0 before the first call.  
//
//  Output:
//
//    int &N_DATA.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    double &U the argument of the function.
//
//    double &A, &K, &M, the parameters of the function.
//
//    double &FX, the value of the function.
//
{
# define N_MAX 20

  static double m_vec[N_MAX] = {
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.0E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     0.5E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00,
     1.0E+00 };

  static double fx_vec[N_MAX] = {
      0.9983341664682815E-01,
      0.1986693307950612E+00,
      0.4794255386042030E+00,
      0.8414709848078965E+00,
      0.9092974268256817E+00,
      0.9975068547462484E-01,
      0.1980217429819704E+00,
      0.4707504736556573E+00,
      0.8030018248956439E+00,
      0.9946623253580177E+00,
      0.9966799462495582E-01,
      0.1973753202249040E+00,
      0.4621171572600098E+00,
      0.7615941559557649E+00,
      0.9640275800758169E+00,
      0.9993292997390670E+00,
     -0.1973753202249040E+00,
     -0.4621171572600098E+00,
     -0.7615941559557649E+00,
     -0.9640275800758169E+00  };

  static double u_vec[N_MAX] = {
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      0.1E+00,
      0.2E+00,
      0.5E+00,
      1.0E+00,
      2.0E+00,
      4.0E+00,
     -0.2E+00,
     -0.5E+00,
     -1.0E+00,
     -2.0E+00 };

  if ( n_data < 0 )
  {
    n_data = 0;
  }

  n_data = n_data + 1;

  if ( N_MAX < n_data )
  {
    n_data = 0;
    a = 0.0;
    k = 0.0;
    m = 0.0;
    u = 0.0;
    fx = 0.0;
  }
  else
  {
    m = m_vec[n_data-1];
    k = sqrt ( m );
    a = asin ( k );
    u = u_vec[n_data-1];
    fx = fx_vec[n_data-1];
  }

  return;
# undef N_MAX
}