rectLeft(f, a, b, n)={
  sum(i=0,n-1,f(a+(b-a)*i/n), 0.)*(b-a)/n
};
rectMid(f, a, b, n)={
  sum(i=1,n,f(a+(b-a)*(i-.5)/n), 0.)*(b-a)/n
};
rectRight(f, a, b, n)={
  sum(i=1,n,f(a+(b-a)*i/n), 0.)*(b-a)/n
};
trapezoidal(f, a, b, n)={
  sum(i=1,n-1,f(a+(b-a)*i/n), f(a)/2+f(b)/2.)*(b-a)/n
};
Simpson(f, a, b, n)={
  my(h=(b - a)/n, s);
  s = 2*sum(i=1,n-1,
    2*f(a + h * (i+1/2)) + f(a + h * i)
  , 0.) + 4*f(a + h/2) + f(a) + f(b);
  s * h / 6
};
test(f, a, b, n)={
  my(v=[rectLeft, rectMid, rectRight, trapezoidal, Simpson]);
  print("Testing function "f" on ",[a,b]," with "n" intervals:");
  for(i=1,#v, print("\t"v[i](f, a, b, n)))
};
# \\ Turn on timer
test(x->x^3, 0, 1, 100)
test(x->1/x, 1, 100, 1000)
test(x->x, 0, 5000, 5000000)
test(x->x, 0, 6000, 6000000)