function square_symq_rule_test03 ( p_lo, p_hi )

%*****************************************************************************80
%
%% square_symq_rule_test03() tests absolute and relative precision.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    16 July 2023
%
%  Author:
%
%    John Burkardt.
%
%  Input:
%
%    integer p_lo, p_hi: the lowest and highest rules to check.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'square_symq_rule_test03():\n' );
  fprintf ( 1, '  Check rules of precision p = %d through %d\n', p_lo, p_hi );
  fprintf ( 1, '  for error in approximating integrals of monomials.\n' );

  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '              maximum                   maximum\n' );
  fprintf ( 1, '   p          absolute                  relative\n' );
  fprintf ( 1, '              error                     error\n' );
  fprintf ( 1, '' );

  for p = p_lo : p_hi

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

    max_abs = 0.0;
    max_rel = 0.0;

    for degree = 0 : p

      expon = [];
      more = 0;
      h = 0;
      t = 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_abs = max ( max_abs, quad_error );
        if ( exact ~= 0.0 )
          max_rel = max ( max_rel, quad_error / abs ( exact ) );
        end

        if ( ~ more )
          break
        end

      end

    end

    fprintf ( 1, '  %2d  %24.16g  %24.16g\n', p, max_abs, max_rel );

  end

  return
end
