BUG: fix romb in scipy.integrate.quadrature

Routine romb for Romberg integration using 2**k +1 samples of a function
was only using 2**(k-1) +1 samples due to wrong ranges.
Renumbered the values R to the common textbook numbering scheme starting
with 0s. Note that romb with 2**k +1 should give the same values as
romberg with divmax=k. Or with 3 values romb is the same as simps.