<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
<p><div class="lev1 toc-item"><a href="#Benchmark-between-Python-and-Julia" data-toc-modified-id="Benchmark-between-Python-and-Julia-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Benchmark between Python and Julia</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method" data-toc-modified-id="The-Romberg-method-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>The Romberg method</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method,-naive-recursive-version-in-Python" data-toc-modified-id="The-Romberg-method,-naive-recursive-version-in-Python-12"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>The Romberg method, naive recursive version in Python</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method,-dynamic-programming-version-in-Python" data-toc-modified-id="The-Romberg-method,-dynamic-programming-version-in-Python-13"><span class="toc-item-num">1.3&nbsp;&nbsp;</span>The Romberg method, dynamic programming version in Python</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method,-better-dynamic-programming-version-in-Python" data-toc-modified-id="The-Romberg-method,-better-dynamic-programming-version-in-Python-14"><span class="toc-item-num">1.4&nbsp;&nbsp;</span>The Romberg method, better dynamic programming version in Python</a></div><div class="lev2 toc-item"><a href="#First-benchmark" data-toc-modified-id="First-benchmark-15"><span class="toc-item-num">1.5&nbsp;&nbsp;</span>First benchmark</a></div><div class="lev2 toc-item"><a href="#Using-Pypy-for-speedup" data-toc-modified-id="Using-Pypy-for-speedup-16"><span class="toc-item-num">1.6&nbsp;&nbsp;</span>Using Pypy for speedup</a></div><div class="lev2 toc-item"><a href="#Numba-version-for-Python" data-toc-modified-id="Numba-version-for-Python-17"><span class="toc-item-num">1.7&nbsp;&nbsp;</span>Numba version for Python</a></div><div class="lev2 toc-item"><a href="#Naive-Julia-version" data-toc-modified-id="Naive-Julia-version-18"><span class="toc-item-num">1.8&nbsp;&nbsp;</span>Naive Julia version</a></div><div class="lev2 toc-item"><a href="#Benchmark-between-Python,-Pypy-and-Julia" data-toc-modified-id="Benchmark-between-Python,-Pypy-and-Julia-19"><span class="toc-item-num">1.9&nbsp;&nbsp;</span>Benchmark between Python, Pypy and Julia</a></div><div class="lev2 toc-item"><a href="#Second-benchmark" data-toc-modified-id="Second-benchmark-110"><span class="toc-item-num">1.10&nbsp;&nbsp;</span>Second benchmark</a></div><div class="lev2 toc-item"><a href="#Conclusion" data-toc-modified-id="Conclusion-111"><span class="toc-item-num">1.11&nbsp;&nbsp;</span>Conclusion</a></div><div class="lev3 toc-item"><a href="#Remark" data-toc-modified-id="Remark-1111"><span class="toc-item-num">1.11.1&nbsp;&nbsp;</span>Remark</a></div>
# Benchmark between Python and Julia
This small [Jupyter notebook](http
Step1:  Let try it with this function $f(x)$ on $[a,b]=[1993,2015]$
Step2:  The first value is the numerical value of the integral $\int_{a}^{b} f(x) \mathrm{d}x$ and the second value is an estimate of the numerical error.
$0.4\%$ is not much, alright.
The Romberg method, naive recursive version in Python
See https
Step3:  It converges quite quickly to the "true" value as given by scipy.integrate.quad.
The Romberg method, dynamic programming version in Python
See https
Step4:  It converges quite quickly to the "true" value as given by scipy.integrate.quad.
The Romberg method, better dynamic programming version in Python
Instead of using a dictionary, which gets filled up dynamically (and so, slowly), let us use a numpy arrays, as we already know the size of the array we need ($n+1 \times m+1$).
Note that only half of the array is used, so we could try to use sparse matrices maybe, for triangular matrices? From what I know, it is not worth it, both in term of memory information (if the sparsity measure is $\simeq 1/2$, you don't win anything from LIL or other sparse matrices representation...
We could use numpy.tri but this uses a dense array so nope.
Step5:  It converges quite quickly to the "true" value as given by scipy.integrate.quad.
First benchmark
Step6:  We already see that the recursive version is much slower than the dynamic programming one!
But there is not much difference between the one using dictionary (romberg()) and the one using a numpy array of a known size (romberg_better()).
Using Pypy for speedup
Step7:  And now the same code executed by an external Pypy interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)
Step8:  This version uses the improved memoization trick (no dictionary), but uses nested lists and not numpy arrays, I didn't bother to install numpy on my Pypy installation (even thought it should be possible).
Numba version for Python
Step9:  It fails! Almost as always when trying Numba, it fails cryptically, too bad. I don't want to spend time debugging this.
Naive Julia version
Thanks to this page for a nice and short introduction to Julia.
Step10:  It seems to work well, like the Python implementation. We get the same numerical result
Step11:  Let try a less naive version using a fixed-size array instead of a dictionary. (as we did before for the Python version)
Step12:  Benchmark between Python, Pypy and Julia
First with Python
Step13:  And now the same code executed by an external Pypy interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)
Step14:  And finally with Julia
Step15:  On this first test, it doesn't look faster than Pypy...
But what if we use the improved version, with an array instead of dictionary?
Step16:  Oh, this time it finally seems faster. Really faster? Yes, about 3 to 4 time faster than Pypy.
Remark also that this last cells compared by using the magic %%pypy and %%script julia, so they both need a warm-up time (opening the pipe, the sub-process, initializing the JIT compiler etc).
But it is fair to compare Pypy to Julia this way.
Second benchmark
Let try the same numerical algorithm but with a different integrand function.
$$\int_{0}^{1} \exp(-x^2) \mathrm{d}x \approx 0.842700792949715$$
First with Python
Step17:  And now the same code executed by an external Pypy interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)
Step18:  And finally with Julia
Step19:  Still not faster than Pypy... So what is the goal of Julia?
<ASSISTANT_TASK:>
Python Code:
from scipy.integrate import quad
Explanation: Table of Contents
<p><div class="lev1 toc-item"><a href="#Benchmark-between-Python-and-Julia" data-toc-modified-id="Benchmark-between-Python-and-Julia-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Benchmark between Python and Julia</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method" data-toc-modified-id="The-Romberg-method-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>The Romberg method</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method,-naive-recursive-version-in-Python" data-toc-modified-id="The-Romberg-method,-naive-recursive-version-in-Python-12"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>The Romberg method, naive recursive version in Python</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method,-dynamic-programming-version-in-Python" data-toc-modified-id="The-Romberg-method,-dynamic-programming-version-in-Python-13"><span class="toc-item-num">1.3&nbsp;&nbsp;</span>The Romberg method, dynamic programming version in Python</a></div><div class="lev2 toc-item"><a href="#The-Romberg-method,-better-dynamic-programming-version-in-Python" data-toc-modified-id="The-Romberg-method,-better-dynamic-programming-version-in-Python-14"><span class="toc-item-num">1.4&nbsp;&nbsp;</span>The Romberg method, better dynamic programming version in Python</a></div><div class="lev2 toc-item"><a href="#First-benchmark" data-toc-modified-id="First-benchmark-15"><span class="toc-item-num">1.5&nbsp;&nbsp;</span>First benchmark</a></div><div class="lev2 toc-item"><a href="#Using-Pypy-for-speedup" data-toc-modified-id="Using-Pypy-for-speedup-16"><span class="toc-item-num">1.6&nbsp;&nbsp;</span>Using Pypy for speedup</a></div><div class="lev2 toc-item"><a href="#Numba-version-for-Python" data-toc-modified-id="Numba-version-for-Python-17"><span class="toc-item-num">1.7&nbsp;&nbsp;</span>Numba version for Python</a></div><div class="lev2 toc-item"><a href="#Naive-Julia-version" data-toc-modified-id="Naive-Julia-version-18"><span class="toc-item-num">1.8&nbsp;&nbsp;</span>Naive Julia version</a></div><div class="lev2 toc-item"><a href="#Benchmark-between-Python,-Pypy-and-Julia" data-toc-modified-id="Benchmark-between-Python,-Pypy-and-Julia-19"><span class="toc-item-num">1.9&nbsp;&nbsp;</span>Benchmark between Python, Pypy and Julia</a></div><div class="lev2 toc-item"><a href="#Second-benchmark" data-toc-modified-id="Second-benchmark-110"><span class="toc-item-num">1.10&nbsp;&nbsp;</span>Second benchmark</a></div><div class="lev2 toc-item"><a href="#Conclusion" data-toc-modified-id="Conclusion-111"><span class="toc-item-num">1.11&nbsp;&nbsp;</span>Conclusion</a></div><div class="lev3 toc-item"><a href="#Remark" data-toc-modified-id="Remark-1111"><span class="toc-item-num">1.11.1&nbsp;&nbsp;</span>Remark</a></div>
# Benchmark between Python and Julia
This small [Jupyter notebook](http://jupyter.org/) shows a simple benchmark comparing various implementations in Python and one in Julia of a specific numerical algorithm, the [Romberg integration method](https://en.wikipedia.org/wiki/Romberg%27s_method).
For Python:
- a recursive implementation,
- a dynamic programming implementation,
- also using Pypy instead,
- (maybe a Numba version of the dynamic programming version)
  + (maybe a Cython version too)
For Julia:
- a dynamic programming implementation will be enough.
----
## The Romberg method
> For mathematical explanations, see [the Wikipedia page](https://en.wikipedia.org/wiki/Romberg%27s_method)
We will use [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function to compare the result of our manual implementations.
End of explanation
import math
f = lambda x: (12*x+1)/(1+math.cos(x)**2)
a, b = 1993, 2017
quad(f, a, b)
Explanation: Let try it with this function $f(x)$ on $[a,b]=[1993,2015]$:
$$ f(x) := \frac{12x+1}{1+\cos(x)^2} $$
End of explanation
def romberg_rec(f, xmin, xmax, n=8, m=None):
    if m is None:  # not m was considering 0 as None
        m = n
    assert n >= m
    if n == 0 and m == 0:
        return ((xmax - xmin) / 2.0) * (f(xmin) + f(xmax))
    elif m == 0:
        h = (xmax - xmin) / float(2**n)
        N = (2**(n - 1)) + 1
        term = math.fsum(f(xmin + ((2 * k) - 1) * h) for k in range(1, N))
        return (term * h) + (0.5) * romberg_rec(f, xmin, xmax, n - 1, 0)
    else:
        return (1.0 / ((4**m) - 1)) * ((4**m) * romberg_rec(f, xmin, xmax, n, m - 1) - romberg_rec(f, xmin, xmax, n - 1, m - 1))
romberg_rec(f, a, b, n=0)  # really not accurate!
romberg_rec(f, a, b, n=1)  # alreay pretty good!
romberg_rec(f, a, b, n=2)
romberg_rec(f, a, b, n=3)
romberg_rec(f, a, b, n=8)  # Almost the exact value.
romberg_rec(f, a, b, n=10)  # Almost the exact value.
romberg_rec(f, a, b, n=12)  # Almost the exact value.
Explanation: The first value is the numerical value of the integral $\int_{a}^{b} f(x) \mathrm{d}x$ and the second value is an estimate of the numerical error.
$0.4\%$ is not much, alright.
The Romberg method, naive recursive version in Python
See https://mec-cs101-integrals.readthedocs.io/en/latest/_modules/integrals.html#romberg_rec for the code and https://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.romberg_rec for the doc
End of explanation
def romberg(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = {(0, 0): 0.5 * (xmax - xmin) * (f(xmax) + f(xmin))}
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / float(2**i)
        xsamples = [xmin + ((2 * k - 1) * h_i) for k in range(1, 1 + 2**(i - 1))]
        r[(i, 0)] = (0.5 * r[(i - 1, 0)]) + h_i * math.fsum(f(x) for x in xsamples)
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            try:
                r[(i, j)] = (((4**j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / float((4**j) - 1)
            except:
                raise ValueError("romberg() with n = {}, m = {} and i = {}, j = {} was an error.".format(n, m, i, j))
    return r[(n, m)]
romberg(f, a, b, n=0)  # really not accurate!
romberg(f, a, b, n=1)  # alreay pretty good!
romberg(f, a, b, n=2)
romberg(f, a, b, n=3)
romberg(f, a, b, n=8)  # Almost the exact value.
romberg(f, a, b, n=10)  # Almost the exact value.
romberg(f, a, b, n=12)  # Almost the exact value.
Explanation: It converges quite quickly to the "true" value as given by scipy.integrate.quad.
The Romberg method, dynamic programming version in Python
See https://mec-cs101-integrals.readthedocs.io/en/latest/_modules/integrals.html#romberg for the code and https://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.romberg for the doc.
It is not hard to make this function non-recursive, by storing the intermediate results.
End of explanation
import numpy as np
def romberg_better(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = np.zeros((n+1, m+1))
    r[0, 0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i, 0] = (0.5 * r[i - 1, 0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i, j] = (((4.**j) * r[i, j - 1]) - r[i - 1, j - 1]) / ((4.**j) - 1.)
    return r[n, m]
romberg_better(f, a, b, n=0)  # really not accurate!
romberg_better(f, a, b, n=1)  # alreay pretty good!
romberg_better(f, a, b, n=2)
romberg_better(f, a, b, n=3)
romberg_better(f, a, b, n=8)  # Almost the exact value.
romberg_better(f, a, b, n=10)  # Almost the exact value.
romberg_better(f, a, b, n=12)  # Almost the exact value.
Explanation: It converges quite quickly to the "true" value as given by scipy.integrate.quad.
The Romberg method, better dynamic programming version in Python
Instead of using a dictionary, which gets filled up dynamically (and so, slowly), let us use a numpy arrays, as we already know the size of the array we need ($n+1 \times m+1$).
Note that only half of the array is used, so we could try to use sparse matrices maybe, for triangular matrices? From what I know, it is not worth it, both in term of memory information (if the sparsity measure is $\simeq 1/2$, you don't win anything from LIL or other sparse matrices representation...
We could use numpy.tri but this uses a dense array so nope.
End of explanation
%timeit quad(f, a, b)
%timeit romberg_rec(f, a, b, n=10)
%timeit romberg(f, a, b, n=10)
%timeit romberg_better(f, a, b, n=10)
Explanation: It converges quite quickly to the "true" value as given by scipy.integrate.quad.
First benchmark
End of explanation
%%time
import numpy as np
import math
import random
f = lambda x: (12*x+1)/(1+math.cos(x)**2)
# Same code
def romberg(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = np.zeros((n+1, m+1))
    r[0, 0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i, 0] = (0.5 * r[i - 1, 0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i, j] = (((4.**j) * r[i, j - 1]) - r[i - 1, j - 1]) / ((4.**j) - 1.)
    return r[n, m]
for _ in range(100000):
    a = random.randint(-2000, 2000)
    b = a + random.randint(0, 100)
    romberg(f, a, b)
Explanation: We already see that the recursive version is much slower than the dynamic programming one!
But there is not much difference between the one using dictionary (romberg()) and the one using a numpy array of a known size (romberg_better()).
Using Pypy for speedup
End of explanation
%%time
%%pypy
import math
import random
f = lambda x: (12*x+1)/(1+math.cos(x)**2)
# Same code
def romberg_pypy(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = [[0 for _ in range(n+1)] for _ in range(m+1)]
    r[0][0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i][0] = (0.5 * r[i - 1][0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i][j] = (((4.**j) * r[i][j - 1]) - r[i - 1][j - 1]) / ((4.**j) - 1.)
    return r[n][m]
for _ in range(100000):
    a = random.randint(-2000, 2000)
    b = a + random.randint(0, 100)
    romberg_pypy(f, a, b)
Explanation: And now the same code executed by an external Pypy interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)
End of explanation
from numba import jit
@jit
def romberg_numba(f, xmin, xmax, n=8):
    assert xmin <= xmax
    m = n
    # First value:
    r = {(0, 0): 0.5 * (xmax - xmin) * (f(xmax) + f(xmin))}
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / float(2**i)
        xsamples = [xmin + ((2 * k - 1) * h_i) for k in range(1, 1 + 2**(i - 1))]
        r[(i, 0)] = (0.5 * r[(i - 1, 0)]) + h_i * math.fsum(f(x) for x in xsamples)
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            try:
                r[(i, j)] = (((4**j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / float((4**j) - 1)
            except:
                raise ValueError("romberg() with n = {}, m = {} and i = {}, j = {} was an error.".format(n, m, i, j))
    return r[(n, m)]
romberg_numba(f, a, b, n=8)  # Almost the exact value.
Explanation: This version uses the improved memoization trick (no dictionary), but uses nested lists and not numpy arrays, I didn't bother to install numpy on my Pypy installation (even thought it should be possible).
Numba version for Python
End of explanation
%%script julia
function f(x)
    (12*x + 1) / (1 + cos(x)^2)
end
a = 1993
b = 2017
function romberg_julia(f, xmin, xmax; n=8)
    m = n
    # First value:
    r = Dict()
    r[(0, 0)] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in 1 : n
        h_i = (xmax - xmin) / (2^i)
        sum_f_x = 0
        for k in 1 : 2^(i - 1)
            sum_f_x += f(xmin + ((2 * k - 1) * h_i))
        end
        r[(i, 0)] = (r[(i - 1, 0)] / 2.) + (h_i * sum_f_x)
    end
    # All the other values:
    for j in 1 : m
        for i in j : n
            r[(i, j)] = (((4^j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / (4^j - 1.)
        end
    end
    r[(n, m)]
end
println(romberg_julia(f, a, b, n=0))  # really not accurate!
println(romberg_julia(f, a, b, n=1))  # alreay pretty good!
println(romberg_julia(f, a, b, n=2))
println(romberg_julia(f, a, b, n=3))
println(romberg_julia(f, a, b, n=8))  # Almost the exact value.
println(romberg_julia(f, a, b, n=10))  # Almost the exact value.
println(romberg_julia(f, a, b, n=12))  # Almost the exact value.
Explanation: It fails! Almost as always when trying Numba, it fails cryptically, too bad. I don't want to spend time debugging this.
Naive Julia version
Thanks to this page for a nice and short introduction to Julia.
End of explanation
f = lambda x: (12*x+1)/(1+math.cos(x)**2)
a, b = 1993, 2017
quad(f, a, b)
romberg(f, a, b, n=12)
Explanation: It seems to work well, like the Python implementation. We get the same numerical result:
End of explanation
%%script julia
function f(x)
    (12*x + 1) / (1 + cos(x)^2)
end
a = 1993
b = 2017
function romberg_julia_better(f, xmin, xmax; n=8)
    m = n
    # First value:
    r = zeros((n+1, m+1))  # https://docs.julialang.org/en/latest/stdlib/arrays/#Base.zeros
    r[1, 1] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in 1 : n
        h_i = (xmax - xmin) / (2^i)
        sum_f_x = 0
        for k in 1 : 2^(i - 1)
            sum_f_x += f(xmin + ((2 * k - 1) * h_i))
        end
        r[i + 1, 1] = (r[i, 1] / 2.) + (h_i * sum_f_x)
    end
    # All the other values:
    for j in 1 : m
        for i in j : n
            r[i + 1, j + 1] = (((4.^j) * r[i + 1, j]) - r[i, j]) / (4.^j - 1.)
        end
    end
    r[n + 1, m + 1]
end
println(romberg_julia_better(f, a, b, n=0))  # really not accurate!
println(romberg_julia_better(f, a, b, n=1))  # alreay pretty good!
println(romberg_julia_better(f, a, b, n=2))
println(romberg_julia_better(f, a, b, n=3))
println(romberg_julia_better(f, a, b, n=8))  # Almost the exact value.
println(romberg_julia_better(f, a, b, n=10))  # Almost the exact value.
println(romberg_julia_better(f, a, b, n=12))  # Almost the exact value.
Explanation: Let try a less naive version using a fixed-size array instead of a dictionary. (as we did before for the Python version)
End of explanation
%%time
import numpy as np
import math
import random
f = lambda x: (12*x+1)/(1+math.cos(x)**2)
# Same code
def romberg(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = np.zeros((n+1, m+1))
    r[0, 0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i, 0] = (0.5 * r[i - 1, 0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i, j] = (((4.**j) * r[i, j - 1]) - r[i - 1, j - 1]) / ((4.**j) - 1.)
    return r[n, m]
for _ in range(100000):
    a = random.randint(-2000, 2000)
    b = a + random.randint(0, 100)
    romberg(f, a, b)
Explanation: Benchmark between Python, Pypy and Julia
First with Python:
End of explanation
%%time
%%pypy
import math
import random
f = lambda x: (12*x+1)/(1+math.cos(x)**2)
# Same code
def romberg_pypy(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = [[0 for _ in range(n+1)] for _ in range(m+1)]
    r[0][0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i][0] = (0.5 * r[i - 1][0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i][j] = (((4.**j) * r[i][j - 1]) - r[i - 1][j - 1]) / ((4.**j) - 1.)
    return r[n][m]
for _ in range(100000):
    a = random.randint(-2000, 2000)
    b = a + random.randint(0, 100)
    romberg_pypy(f, a, b)
Explanation: And now the same code executed by an external Pypy interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)
End of explanation
%%time
%%script julia
function f(x)
    (12*x + 1) / (1 + cos(x)^2)
end
function romberg_julia(f, xmin, xmax; n=8)
    m = n
    # First value:
    r = Dict()
    r[(0, 0)] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in 1 : n
        h_i = (xmax - xmin) / (2^i)
        sum_f_x = 0
        for k in 1 : 2^(i - 1)
            sum_f_x += f(xmin + ((2 * k - 1) * h_i))
        end
        r[(i, 0)] = (r[(i - 1, 0)] / 2.) + (h_i * sum_f_x)
    end
    # All the other values:
    for j in 1 : m
        for i in j : n
            r[(i, j)] = (((4^j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / (4^j - 1.)
        end
    end
    r[(n, m)]
end
for _ in 1:100000
    a = rand(-2000:2000)
    b = a + rand(0:100)
    romberg_julia(f, a, b)
end
Explanation: And finally with Julia:
End of explanation
%%time
%%script julia
function f(x)
    (12*x + 1) / (1 + cos(x)^2)
end
function romberg_julia_better(f, xmin, xmax; n=8)
    m = n
    # First value:
    r = zeros((n+1, m+1))  # https://docs.julialang.org/en/latest/stdlib/arrays/#Base.zeros
    r[1, 1] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in 1 : n
        h_i = (xmax - xmin) / (2^i)
        sum_f_x = 0
        for k in 1 : 2^(i - 1)
            sum_f_x += f(xmin + ((2 * k - 1) * h_i))
        end
        r[i + 1, 1] = (r[i, 1] / 2.) + (h_i * sum_f_x)
    end
    # All the other values:
    for j in 1 : m
        for i in j : n
            r[i + 1, j + 1] = (((4.^j) * r[i + 1, j]) - r[i, j]) / (4.^j - 1.)
        end
    end
    r[n + 1, m + 1]
end
for _ in 1:100000
    a = rand(-2000:2000)
    b = a + rand(0:100)
    romberg_julia_better(f, a, b)
end
Explanation: On this first test, it doesn't look faster than Pypy...
But what if we use the improved version, with an array instead of dictionary?
End of explanation
%%time
import numpy as np
import math
import random
f = lambda x: (2.0 / math.sqrt(math.pi)) * math.exp(-x**2)
# Same code
def romberg(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = np.zeros((n+1, m+1))
    r[0, 0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i, 0] = (0.5 * r[i - 1, 0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i, j] = (((4.**j) * r[i, j - 1]) - r[i - 1, j - 1]) / ((4.**j) - 1.)
    return r[n, m]
for _ in range(100000):
    a = 0
    b = 1
    romberg(f, a, b)
print(romberg(f, a, b))
Explanation: Oh, this time it finally seems faster. Really faster? Yes, about 3 to 4 time faster than Pypy.
Remark also that this last cells compared by using the magic %%pypy and %%script julia, so they both need a warm-up time (opening the pipe, the sub-process, initializing the JIT compiler etc).
But it is fair to compare Pypy to Julia this way.
Second benchmark
Let try the same numerical algorithm but with a different integrand function.
$$\int_{0}^{1} \exp(-x^2) \mathrm{d}x \approx 0.842700792949715$$
First with Python:
End of explanation
%%time
%%pypy
import math
import random
f = lambda x: (2.0 / math.sqrt(math.pi)) * math.exp(-x**2)
# Same code
def romberg_pypy(f, xmin, xmax, n=8, m=None):
    assert xmin <= xmax
    if m is None:
        m = n
    assert n >= m >= 0
    # First value:
    r = [[0 for _ in range(n+1)] for _ in range(m+1)]
    r[0][0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in range(1, n + 1):
        h_i = (xmax - xmin) / 2.**i
        r[i][0] = (0.5 * r[i - 1][0]) + h_i * math.fsum(
            f(xmin + ((2 * k - 1) * h_i))
            for k in range(1, 1 + 2**(i - 1))
        )
    # All the other values:
    for j in range(1, m + 1):
        for i in range(j, n + 1):
            r[i][j] = (((4.**j) * r[i][j - 1]) - r[i - 1][j - 1]) / ((4.**j) - 1.)
    return r[n][m]
for _ in range(100000):
    a = 0
    b = 1
    romberg_pypy(f, a, b)
print(romberg_pypy(f, a, b))
Explanation: And now the same code executed by an external Pypy interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)
End of explanation
%%time
%%script julia
function f(x)
    (2.0 / sqrt(pi)) * exp(-x^2)
end
function romberg_julia(f, xmin, xmax; n=8)
    m = n
    # First value:
    r = Dict()
    r[(0, 0)] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in 1 : n
        h_i = (xmax - xmin) / (2^i)
        sum_f_x = 0
        for k in 1 : 2^(i - 1)
            sum_f_x += f(xmin + ((2 * k - 1) * h_i))
        end
        r[(i, 0)] = (r[(i - 1, 0)] / 2.) + (h_i * sum_f_x)
    end
    # All the other values:
    for j in 1 : m
        for i in j : n
            r[(i, j)] = (((4^j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / (4^j - 1.)
        end
    end
    r[(n, m)]
end
for _ in 1:100000
    a = 0
    b = 1
    romberg_julia(f, a, b)
end
println(romberg_julia(f, 0, 1))
Explanation: And finally with Julia:
End of explanation
%%time
%%script julia
function f(x)
    (2.0 / sqrt(pi)) * exp(-x^2)
end
function romberg_julia_better(f, xmin, xmax; n=8)
    m = n
    # First value:
    r = zeros((n+1, m+1))  # https://docs.julialang.org/en/latest/stdlib/arrays/#Base.zeros
    r[1, 1] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.
    # One side of the triangle:
    for i in 1 : n
        h_i = (xmax - xmin) / (2^i)
        sum_f_x = 0
        for k in 1 : 2^(i - 1)
            sum_f_x += f(xmin + ((2 * k - 1) * h_i))
        end
        r[i + 1, 1] = (r[i, 1] / 2.) + (h_i * sum_f_x)
    end
    # All the other values:
    for j in 1 : m
        for i in j : n
            r[i + 1, j + 1] = (((4.^j) * r[i + 1, j]) - r[i, j]) / (4.^j - 1.)
        end
    end
    r[n + 1, m + 1]
end
for _ in 1:100000
    a = 0
    b = 1
    romberg_julia_better(f, a, b)
end
println(romberg_julia_better(f, 0, 1))
Explanation: Still not faster than Pypy... So what is the goal of Julia?
End of explanation
<END_TASK>