(*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 of 
degree 2 
           n + 1 for integral of predefined
function f[
            x, y] over region between
predefined functions
cc (x) and dd (x),  
      for a ≤ x ≤ 
        b; 
  II is the estimate*) 
G2dim[n_, a_, 
        b_]:=(
Gaussian[n, a, 
          b]; 
For[i = 0, 
          i ≤ n, i++, 
xx[i] = 
            x[i]; 
xweight[i] =  
            weights[[i + 1]]]; 
For[ii = 0, 
          ii ≤ n, ii++, 
 
          Gaussian[n,  cc[xx[ii]], dd[xx[ii]]]; 
For[
            jj = 0, jj ≤ 
              n, jj++, 
Ytab[ii, jj] = 
              x[jj]; 
Wtab[ii, jj]  =  
              weights[[jj + 1]]]]; 
II = 
          0; 
For[iii = 0, 
          iii ≤ n, iii++, 
Ii  =  
            0; 
For[jjj = 0, 
            jjj ≤ n, jjj++, 
Ii  =  
              Ii  +   Wtab[iii, jjj] f[xx[iii], 
                    Ytab[iii, jjj]]]; 
   II  =  
            II  +  xweight[
                  iii] Ii];
)

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

a = 1

b = 2

cc[x_] := 3;

dd[x_] := 5;

n = 3

cc[xx[1]]

dd[xx[1]]

NIntegrate[f[x, y], {x, a, b}, {y, 3, 5}]

G2dim[n, a, b]

II

f[x_, y_] := 1;

a = -1

b = 1

cc[x_] := -Sqrt[1 - x^2];

dd[x_] := Sqrt[1 - x^2];

n = 7

NIntegrate[
      f[x, y], {x, -1, 1}, {y, -Sqrt[1 - x^2], 
        Sqrt[1 - x^2]}]

G2dim[n, a, b]

II

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

cc[x_] =  -x;

dd[x_] :=  x^2

a = 2

b = 5

n = 3

NIntegrate[f[x, y], {x, 2, 5}, {y, cc[x], dd[x]}]

G2dim[n, a, b]

II

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

NIntegrate[f[x, y], {x, 2, 5}, {y, cc[x], dd[x]}]

G2dim[3, a, b]

II

n = 4

G2dim[n, a, b]

II

G2dim[5, a, b]

II

G2dim[6, a, b]

II

Abs[Out[145] - Out[165]]/Out[145]