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.
TST: integrate: add test for romb accuracy improvement
Merge pull request #3731 from 'pr-3731'

BUG: interpolate: fix wrong romb sample usage

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.