(*Compute the densities and nodes 
for Gaussian quadrature 
of degree (precision) 2  n + 
        1 over [a, b]; 
  generates n + 
        1 nodes and densities . 
Based on R . 
            Curto and L . 
            Fialkow
[
              Houston J .  Math .  17 (1991),  
              603 - 635], 
      MR 93  a : 
          47016*)         
Gaussian[
        n_, a_, b_] := (
For[i = 0, 
          i ≤ 
            2  n + 1, i++, 
bb[
              i]  =  (b^(i + 1) - a^(i + 1))/(i + 
                  1)]; 
For[i = 0, i ≤ n, 
          i = i + 1, 
For[j = 0, j ≤ n, 
            j = j + 1, 
Mtab[i, j] = 
              bb[i + j]]]; 
H = 
          Table[Mtab[r, s], {r, 0, n}, {s, 0, n}]; 
J  =  
          N[Inverse[H]]; 
For[i = 0,  i ≤ n, 
          i = i + 1, 
v[i]  =  
            bb[n + 1 + i]]; 
w  = 
          Table[v[i], {i, 0, n}]; 
z  =  
          J . w; 
p  =  0; 
For[i = 0, 
          i ≤ n, 
          i = i + 1, 
p  =  
            p  +  z[[i + 1]] t^i]; 
q  =  
          t^(n + 1) - p; 
roots  =  
          Solve[q == 0, t]; 
For[i = 1, 
          i ≤ n + 1, 
          i = i + 1, 
x[
              i - 1]  =  roots[[i]][[1]][[2]]]; 
For[ i = 0, i ≤ n, 
          i = i + 1, 
For[j = 0, j ≤ n, 
            j = j + 1, 
TTab[i, j] =  
              x[j]^i]]; 
V =  
          Table[TTab[rr, ss], {rr, 0, n}, {ss, 0, n}]; 
For[
          i = 0,  i ≤ n, 
          i = i + 1, 
vvv[i]  =  
            bb[i]]; 
w  = 
          Table[vvv[i], {i, 0, n}]; 
weights = 
          N[Inverse[V] . w];)

(*Iterated Gaussian quadrature estimate of 
degree 2  n + 
        1 for integral of predefined
function f[x, 
            y] over rectangle [a, b] x [c, 
            d]; 
  exact if f is a polynomial whose degrees in x and y
 are each  ≤  
        2  n + 1;  
      the value is stored in II*) 
G2dim[n_, 
        a_, b_, c_, 
        d_]:=(
Gaussian[n, a, 
          b]; 
For[i = 0, 
          i ≤ n, i++, 
xx[i] = 
            x[i]; 
xweight[i] =  
            weights[[i + 1]]]; 
Gaussian[n, c, 
          d]; 
For[i = 0, 
          i ≤ n, i++, 
yy[i] = 
            x[i]; 
yweight[i]  =  
            weights[[i + 1]]]; 
II = 
          0; 
For[j = 0, 
          j ≤ n, j++, 
For[i = 0, 
            i ≤ n, i++, 
II  =  
              II  +  xweight[j] yweight[i] f[xx[j], 
                    yy[i]]]];
)

f[x_, y_] :=  x^7 y^7;

a = 1

b = 2

c = 3

d = 5

n = 3

NIntegrate[f[x, y], {x, a, b}, {y, c, d}]

G2dim[n, a, b, c, d]

II

f[x_, y_] :=  x^9 y^3;

NIntegrate[f[x, y], {x, a, b}, {y, c, d}]

G2dim[n, a, b, c, d]

II

f[x_, y_] :=  x^10 y^4

NIntegrate[f[x, y], {x, a, b}, {y, c, d}]

G2dim[n, a, b, c, d]

II

f[x_, y_] :=  Exp[Sin[x]  + y];

NIntegrate[f[x, y], {x, a, b}, {y, c, d}]

G2dim[n, a, b, c, d]

II

f[x_, y_] :=  Sin[Exp[x^2  +  y^2]];

NIntegrate[f[x, y], {x, a, b}, {y, c, d}]

G2dim[n, a, b, c, d]

II