# Return definite integral of a bivariate 
# polynomial in BB-form over the unit square:
# int_{unit sqare} = sum b_{ij} / ((du+1)!/(du!*1!)) / ((dv+1)!/(dv!*1!))
#
# d = degree of polynomial
# b = cooeficients of polynomial in table form

bb4S := proc(b,du,dv)

local out,i,j:
    out := sum(sum(b[i,du-i][j,dv-j],j=0..du),i=0..dv);
    out := out / ((du+1)*(dv+1))
end:
#--------------------
test := 0:
if test=1 then
    a[0,1][0,1] := 0:
    a[1,0][0,1] := 0:
    a[0,1][1,0] := 1:
    a[1,0][1,0] := 1:
    bbS := bb4S(a,1,1);
    print(bbS);
fi:
if test=2 then
    a[0,1][0,1] := 0:
    a[1,0][0,1] := 1:
    a[0,1][1,0] := 1:
    a[1,0][1,0] := 2:
    bbS := bb4S(a,1,1);
    print(bbS);
fi: