Our main result is a significantly improved algorithm for this problem and the d-dimensional analog. Our algorithm outputs an R with the same approximation guarantees but it runs in time B*poly(d,log(N),log(M),1/\epsilon).
We need two crucial ideas to achieve this bound: bulk sampling and
estimation for multipoint polynomial evaluation using an
unevenly-spaced Fourier tranform, and construction and use of
arithmetic-progression independent random variables.
Paper Available at:
ftp://dimacs.rutgers.edu/pub/dimacs/TechnicalReports/TechReports/2004/2004-49.ps.gz