MidpointOpenNewtonCotesRule:=n->factor(int(interp([seq((i-1/2)*h,
              i=1..n)], [seq(f((i-1/2)*h),i=1..n)],z),z=0..n*h)):