# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "square_symq_rule.h"

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

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

int main ( )

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

    square_symq_rule_test() tests square_symq_rule().

  Licensing:

    This code is distributed under the GNU GPL license.

  Modified:

    15 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 ( );
  printf ( "\n" );
  printf ( "square_symq_rule_test():\n" );
  printf ( "  C version\n" );
  printf ( "  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.
*/
  printf ( "\n" );
  printf ( "square_symq_rule_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void test01 ( int p )

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

    test01() prints a rule of precision P.

  Licensing:

    This code is distributed under the GNU GPL license.

  Modified:

    15 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;

  printf ( "\n" );
  printf ( "test01():\n" );
  printf ( "  Symmetric quadrature rule for a square.\n" );
  printf ( "  precision = %d\n", p );
/*
  Retrieve and print a symmetric quadrature rule.
*/
  n = rule_order ( p );

  x = ( double * ) malloc ( n * sizeof ( double ) );
  y = ( double * ) malloc ( n * sizeof ( double ) );
  w = ( double * ) malloc ( n * sizeof ( double ) );

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

  printf ( "\n" );
  printf ( "  Number of nodes N = %d\n", n );

  printf ( "\n" );
  printf ( "     J  W       X       Y\n" );
  printf ( "\n" );
  for ( j = 0; j < n; j++ )
  {
    printf ( "  %4d  %14.6g  %14.6g  %14.6g\n",
      j, w[j], x[j], y[j] );
  }

  w_sum = r8vec_sum ( n, w );

  printf ( "  Weight sum  %g\n", w_sum );

  free ( x );
  free ( y );
  free ( w );

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

void test02 ( int p )

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

    test02() tests a rule of precision P.

  Licensing:

    This code is distributed under the GNU GPL license.

  Modified:

    15 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.  0 <= DEGREE <= 50.
*/
{
  int degree;
  const int dim_num = 2;
  double exact;
  int expon[2];
  int h;
  int i;
  int ij;
  double max_error;
  int more;
  int n;
  double quad;
  double quad_error;
  int t;
  double *v;
  double *w;
  double *x;
  double *xy;
  double *y;

  printf ( "\n" );
  printf ( "test02():\n" );
  printf ( "  Test a rule of precision %d.\n", p );
/*
  Retrieve a symmetric quadrature rule.
*/
  n = rule_order ( p );

  x = ( double * ) malloc ( n * sizeof ( double ) );
  y = ( double * ) malloc ( n * sizeof ( double ) );
  w = ( double * ) malloc ( n * sizeof ( double ) );

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

  printf ( "\n" );
  printf ( "  Stated precision of rule    = %d\n", p );
  printf ( "  Number of quadrature points = %d\n", n );
  printf ( "\n" );
  printf ( "  Degree  Maximum error\n" );
  printf ( "\n" );

  v = ( double * ) malloc ( n * sizeof ( double ) );

  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;
      }
    }
    printf ( "  %2d  %24.16g\n", degree, max_error );
  }
/*
  Free memory.
*/
  free ( v );
  free ( w );
  free ( x );
  free ( xy );
  free ( y );

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

void test03 ( int p_lo, int p_hi )

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

    test03() tests absolute and relative precision.

  Licensing:

    This code is distributed under the GNU GPL license.

  Modified:

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

    integer 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;
  int more;
  int n;
  int p;
  double quad;
  double quad_error;
  int t;
  double *v;
  double *w;
  double *x;
  double *xy;
  double *y;

  printf ( "\n" );
  printf ( "test02():\n" );
  printf ( "  Test the precision of a quadrature rule for the unit square.\n" );
  printf ( "  Check rules of precision p = %d through' %d\n", p_lo, p_hi );
  printf ( "  for error in approximating integrals of monomials.\n" );

  printf ( "\n" );
  printf ( "              maximum                   maximum\n" );
  printf ( "   p          absolute                  relative\n" );
  printf ( "              error                     error\n" );
  printf ( "\n" );

  for ( p = p_lo; p <= p_hi; p++ )
  {
/*
  Retrieve a symmetric quadrature rule.
*/
    n = rule_order ( p );

    x = ( double * ) malloc ( n * sizeof ( double ) );
    y = ( double * ) malloc ( n * sizeof ( double ) );
    w = ( double * ) malloc ( n * sizeof ( double ) );

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

    v = ( double * ) malloc ( n * sizeof ( double ) );

    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;
        }
      }
    }
    printf ( "  %2d  %24.16g  %24.16g\n", p, max_abs, max_rel );

    free ( v );
    free ( w );
    free ( x );
    free ( xy );
    free ( y );
  }

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

void timestamp ( )

/******************************************************************************/
/*
  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:

    24 September 2003

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

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

  now = time ( NULL );
  tm = localtime ( &now );

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

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}