# include <cmath>
# include <cstdlib>
# include <cstring>
# include <ctime>
# include <fstream>
# include <iomanip>
# include <iostream>

using namespace std;

# include "square_symq_rule.hpp"

int main ( );
void test01 ( int p );
void test02 ( int p );
void test03 ( int p_lo, int p_hi );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    square_symq_rule_test() tests square_symq_rule().
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    16 July 2023
//
//  Author:
//
//    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
//    This version by John Burkardt.
//
//  Reference:
//
//    Hong Xiao, Zydrunas Gimbutas,
//    A numerical algorithm for the construction of efficient quadrature
//    rules in two and higher dimensions,
//    Computers and Mathematics with Applications,
//    Volume 59, 2010, pages 663-676.
//
{
  int p;
  int p_hi;
  int p_lo;

  timestamp ( );
  cout << "\n";
  cout << "square_symq_rule_test():\n";
  cout << "  C++ version\n";
  cout << "  Test square_symq_rule().\n";

  p = 5;
  test01 ( p );

  p = 5;
  test02 ( p );

  p_lo = 0;
  p_hi = 20;
  test03 ( p_lo, p_hi );
//
//  Terminate.
//
  cout << "\n";
  cout << "square_symq_rule_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void test01 ( int p )

//****************************************************************************80
//
//  Purpose:
//
//    test01() prints a rule of precision P.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    16 July 2023
//
//  Author:
//
//    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
//    This version by John Burkardt.
//
//  Reference:
//
//    Hong Xiao, Zydrunas Gimbutas,
//    A numerical algorithm for the construction of efficient quadrature
//    rules in two and higher dimensions,
//    Computers and Mathematics with Applications,
//    Volume 59, 2010, pages 663-676.
//
//  Input:
//
//    int p: the precision of the rule.
//
{
  int j;
  int n;
  double *w;
  double w_sum;
  double *x;
  double *y;

  cout << "\n";
  cout << "test01():\n";
  cout << "  Symmetric quadrature rule for a square.\n";
  cout << "  Precision = " << p << "\n";
//
//  Retrieve and print a symmetric quadrature rule.
//
  n = rule_order ( p );

  x = new double[n];
  y = new double[n];
  w = new double[n];

  square_symq ( p, n, x, y, w );

  cout << "\n";
  cout << "  Number of nodes N = " << n << "\n";

  cout << "\n";
  cout << "     J  W       X       Y\n";
  cout << "\n";
  for ( j = 0; j < n; j++ )
  {
    cout << setw(4) << j << "  "
         << setw(14) << w[j] << "  "
         << setw(14) << x[j] << "  "
         << setw(14) << y[j] << "\n";
  }

  w_sum = r8vec_sum ( n, w );

  cout << "  Weight sum  " << w_sum << "\n";

  delete [] w;
  delete [] x;
  delete [] y;

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

void test02 ( int p )

//****************************************************************************80
//
//  Purpose:
//
//    test02() tests a rule of precision P.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    16 July 2023
//
//  Author:
//
//    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
//    This version by John Burkardt.
//
//  Reference:
//
//    Hong Xiao, Zydrunas Gimbutas,
//    A numerical algorithm for the construction of efficient quadrature
//    rules in two and higher dimensions,
//    Computers and Mathematics with Applications,
//    Volume 59, 2010, pages 663-676.
//
//  Input:
//
//    int p: the precision of the rule.
//
{
  int degree;
  const int dim_num = 2;
  double exact;
  int expon[2];
  int h;
  int i;
  int ij;
  double max_error;
  bool more;
  int n;
  double quad;
  double quad_error;
  int t;
  double *v;
  double *w;
  double *x;
  double *xy;
  double *y;

  cout << "\n";
  cout << "test02():\n";
  cout << "  Get a quadrature rule for the symmetric square.\n";
  cout << "  Test its accuracy.\n";
  cout << "  Precision = " << p << "\n";
//
//  Retrieve the quadrature rule.
//
  n = rule_order ( p );

  w = new double[n];
  x = new double[n];
  y = new double[n];

  square_symq ( p, n, x, y, w );
//
//  Pack the x, y vectors as columns of an array.
//
  xy = new double[n*2];
  ij = 0;
  for ( i = 0; i < n; i++ )
  {
    xy[ij] = x[i];
    ij = ij + 1;
    xy[ij] = y[i];
    ij = ij + 1;
  }

  cout << "\n";
  cout << "  Number of quadrature points = " << n << "\n";
  cout << "\n";
  cout << "  Degree  Maximum error\n";
  cout << "\n";

  v = new double[n];

  for ( degree = 0; degree <= p + 2; degree++ )
  {
    more = 0;
    h = 0;
    t = 0;
    max_error = 0.0;

    while ( true )
    {
      comp_next ( degree, dim_num, expon, more, h, t );

      v = monomial_value ( dim_num, n, expon, xy );

      quad = 0.0;
      for ( i = 0; i < n; i++ )
      {
        quad = quad + w[i] * v[i];
      }
      quad = quadrilateral_unit_area ( ) * quad;

      exact = quadrilateral_unit_monomial_integral ( expon );

      quad_error = fabs ( quad - exact );
      max_error = fmax ( max_error, quad_error );

      if ( ! more )
      {
        break;
      }
    }
    cout << "  " << setw(2) << degree
         << "  " << setw(24) << max_error << "\n";
  }
//
//  Free memory.
//
  delete [] v;
  delete [] w;
  delete [] x;
  delete [] xy;
  delete [] y;

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

void test03 ( int p_lo, int p_hi )

//****************************************************************************80
//
//  Purpose:
//
//    test03() tests absolute and relative precision.
//
//  Licensing:
//
//    This code is distributed under the GNU GPL license.
//
//  Modified:
//
//    16 July 2023
//
//  Author:
//
//    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
//    This version by John Burkardt.
//
//  Reference:
//
//    Hong Xiao, Zydrunas Gimbutas,
//    A numerical algorithm for the construction of efficient quadrature
//    rules in two and higher dimensions,
//    Computers and Mathematics with Applications,
//    Volume 59, 2010, pages 663-676.
//
//  Input:
//
//    int p_lo, p_hi: the lowest and highest rules to check.
//
{
  int degree;
  const int dim_num = 2;
  double exact;
  int expon[2];
  int h;
  int i;
  int ij;
  double max_abs;
  double max_rel;
  bool more;
  int n;
  int p;
  double quad;
  double quad_error;
  int t;
  double *v;
  double *w;
  double *x;
  double *xy;
  double *y;

  cout << "\n";
  cout << "test03():\n";
  cout << "  Test the precision of quadrature rules for the unit quadrilateral.\n";
  cout << "  Check rules of precision p = " << p_lo << " through " << p_hi << "\n";
  cout << "  for error in approximating integrals of monomials.\n";

  cout << "\n";
  cout << "              maximum                   maximum\n";
  cout << "   p          absolute                  relative\n";
  cout << "              error                     error\n";
  cout << "\n";

  for ( p = p_lo; p <= p_hi; p++ )
  {
    n = rule_order ( p );

    w = new double[n];
    x = new double[n];
    y = new double[n];

    square_symq ( p, n, x, y, w );
//
//  Pack the x, y vectors as columns of an array.
//
    xy = new double[n*2];
    ij = 0;
    for ( i = 0; i < n; i++ )
    {
      xy[ij] = x[i];
      ij = ij + 1;
      xy[ij] = y[i];
      ij = ij + 1;
    }

    v = new double[n];

    max_abs = 0.0;
    max_rel = 0.0;

    for ( degree = 0; degree <= p; degree++ )
    {
      more = 0;
      h = 0;
      t = 0;

      while ( true )
      {
        comp_next ( degree, dim_num, expon, more, h, t );

        v = monomial_value ( dim_num, n, expon, xy );

        quad = 0.0;
        for ( i = 0; i < n; i++ )
        {
          quad = quad + w[i] * v[i];
        }
        quad = quadrilateral_unit_area ( ) * quad;

        exact = quadrilateral_unit_monomial_integral ( expon );

        quad_error = fabs ( quad - exact );
        max_abs = fmax ( max_abs, quad_error );
        if ( exact != 0.0 )
        {
          max_rel = fmax ( max_rel, quad_error / fabs ( exact ) );
        }

        if ( ! more )
        {
          break;
        }
      }

    }
    cout << "  " << setw(2) << p
         << "  " << setw(24) << max_abs
         << "  " << setw(24) << max_rel << "\n";

    delete [] v;
    delete [] w;
    delete [] x;
    delete [] xy;
    delete [] y;
  }

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

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}