function square_symq_rule_test04 ( p )

%*****************************************************************************80
%
%% square_symq_rule_test04() tests the precision of a quadrilateral quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    08 June 2023
%
%  Author:
%
%    John Burkardt.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'square_symq_rule_test04():\n' );
  fprintf ( 1, '  Test the precision of a quadrature rule for the square.\n' );

  dim_num = 2;
%
%  Retrieve the rule.
%
  n = rule_order ( p );
  [ x, y, w ] = square_symq ( p, n );
%
%  Pack the x, y vectors as columns of an array.
%
  xy = [ x', y' ];

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

  for degree = 0 : p + 2

    expon = [];
    more = 0;
    h = 0;
    t = 0;
    max_error = 0.0;

    while ( true )

      [ expon, more, h, t ] = comp_next ( degree, dim_num, expon, more, h, t );

      v = monomial_value ( expon, xy );

      quad = quadrilateral_unit_area ( ) * dot ( w, v );

      exact = quadrilateral_unit_monomial_integral ( expon );

      quad_error = abs ( quad - exact );
      max_error = max ( max_error, quad_error );

      if ( ~ more )
        break
      end

    end

    fprintf ( 1, '  %2d  %24.16f\n', degree, max_error );

  end

  return
end
