The problem¶

It is a very well known fact that integration can be a really, really hard process; the usual way of computing the exact value of a definite inegral requires finding an antiderivative for the integrand, evaluating it at the boundary points, and then subtracting the results. The problem is that this first step can be very difficult, if not impossible¹, since the primitives of many simple functions such as $\e^{-x^2}$ or $x^x$ cannot be written in terms of elementary functions (sines, cosines, exponentials, inverse trigonometric functions, ...).

So what then? Sometimes Taylor-expanding the function and then interchanging the infinite sum with the integral via the dominated convergence theorem works, others fancy theorems such as Glasser's and Ramanujan's master theorems may do the trick, or even moving up to the complex plane and using the residue theorem there. When all the previous fail, though, there is no choice but to resort to numerical integration, that is, giving up the exact result and instead try to find a very precise estimate with the aid of a computer or calculator.

The aim of this post is to explain the integration method developed by Turing Award laureate William Kahan for HP calculators in the 1980s (HP-34C, HP-15C, HP-42S), but this will be the culmination of an introduction to numerical analysis (with a special emphasis on numerical integration).

A word² about errors¶

Computers are undoubtedly very poweful machines, capable of performing computations that would take days (if not months or years) to a regular human in a fraction of a second. However, because of their own nature as physical, material objects, they are still constrained by space and time, so they can only manipulate numbers with a finite amount of digits.

This is a problem, and a big one since the vast majority of real numbers are irrational, including ubiquitous constants such as $\pi$ or the aforementioned $\e$, and even periodic rational numbers like $1/3$ or $2/7$ cannot be represented with a finite precision.

Floating-point numbers¶

The current standard that dictates how the so-called floating-point numbers are stored and handled is the IEEE 754 (a number that you might eventually end up memorizing), which works like scientific notation but in binary. For instance, let us examine the 32 bits of the IEEE 754 representation of $2.875$: $$ 2.875 = 1.4375 \cdot 2^1 = \texttt{0b} \underbrace{\color{blue} \texttt{0}}_{\text{Sign bit}}\: \underbrace{\color{orange} \texttt{10000000}}_{\text{Exponent}}\: \underbrace{\color{green} \texttt{01110000000000000000000}}_{\text{Mantissa}}; $$ the mantissa encodes the decimal part of the factor, since $0.4375=0\cdot2^{-1} + 1\cdot2^{-2} + 1\cdot2^{-3} + 1\cdot2^{-4}$, whereas the exponent bits store, well, the exponent of the power of $2$ multiplying; finally, the sign bit is set to 0 when the number is positive, and 1 when it is negative.

There are also 64-bit floats, and if your computer has an x86 CPU with a floating-point unit (it most likely does) it even supports 80-bit floating-point numbers; if you know some C, the three formats described above correspond to the float, double and long double types when available. Most interpreted languages such as Python and JavaScript (the latter is the one running in the previous console) make things easy by keeping all decimal numbers 64-bit floats by default, so why not go up again and try executing IEEE754f64(2.875) this time to see how the number is actually saved in memory?

Truncation and overflow/underflow errors¶

Anyway, this is not the topic of this section, but rather what happens when you operate with two floats and the result cannot be stored within the given size: the number gets truncated and the mantissa rounds to the nearest even number when interpreted as an integer (i.e., it keeps its last bit equal to 0). This means that our computations lose precision, and whilst rounding $2.19856478212575$ to $2.1985647821258$ is not a big deal if done only once, the error can accumulate, and after many iterations end up giving a completely wrong result.

A tragic example of this phenomenon is the Patriot Missile incident in Dharan: Three decades ago, back during the times of the Second Gulf War, the US military had installed a battery of defense missiles in Dharan (South Arabia), but it failed to intercept an Iraqi Scud missile due to computer arithmetic errors, causing the death of 28 soldiers and over a hundred other people. Specifically, the system's internal clock measured time in tenths of a second, so the error came when multiplying by $1/10$ to perform the conversion; all of this in a fixed-point³ 24-bit register.

Another type of behavior floating-point numbers may exhibit when operating with them are the so-called overflow and underflow errors. When adding two very large numbers, their sum may exceed the maximum value of the corresponding float type, and what happens in this case is that the output rolls back to a nonsense result. Likewise, when subtracting two very close floats, if the infimum result cannot be represented precisely, it may round down to zero.

For instance, recall that $$ \frac{\diff f}{\diff x} = \lim_{h\to 0} \frac{f(x + h) - f(x)}h \approx \frac{f(x + \varepsilon) - f(x)}\varepsilon $$ for some very small value of $\varepsilon$. However, doing this in practice conveys the risk of having $f(x+\varepsilon)-f(x)$ underflow, and hence incorrectly yielding $f'(x) = 0$.

Asymptotic notation¶

If errors are unavoidable in numerical computations, can we at least estimate and control them? The answer is thankfully affirmative, and mathematicians have developed a special set of notations for that. Let us first delve into the formal definition of the most famous of them all, the big-O:

Intuitively, having $f=\mathscr O(g)$ means that $f$ is bounded in magnitude from above by $g$ (possibly scaled by a factor), and in particular that when $x$ goes to infinity, $f$ grows as fast or slower than $g$. You might have seen this in computer science when analyzing the time complexity of an algorithm, and indeed it is the exact same concept, where $f$ in this case is the time (or memory, or any other resource) the algorithm needs to run.

As it was previously stated, there are other symbols in the literature that play similar roles (mostly due to Donald Knuth), for instance $f=\Omega(g)$ to denote a lower bound, or $f=\Theta(g)$ to represent an exact bound.⁴

But again, we came here seeking to somewhat understand the error of our algorithms, and to see how the definitions above make sense when doing mathematics, the example of the Taylor polynomial is very enlightening. Recall first Taylor's theorem:

As much as I dislike not proving the results I use, from the fact above it follows, for example, that $$ \begin{alignat*}{5} \sin x &= x &&-\frac{x^3}6 &&+\frac{x^5}{120} &&-\frac{x^7}{5040} &&+\mathscr O(x^9),\\ \cos (\pi x) &= 1 &&-\frac{(\pi x)^2}2 &&+\frac{(\pi x)^4}{24} &&-\frac{(\pi x)^6}{720} &&+\mathscr O((\pi x)^8),\\ \e^{x-4} &= 1 &&+(x-4) &&+\frac{(x-4)^2}2 &&+\frac{(x-4)^3}6 &&+\mathscr O((x-4)^4); \end{alignat*} $$ we have approximated three previous functions whilst still retaining some information about how big can be the difference with respect to the actual value.

If real analysis is the art of bounding, numerical analysis is the art of approximating.

Newton-Cotes formulae¶

Now comes the time to talk about integration itself, but before bringing up the clever methods developed for this purpose, it is worth going back and remembering the concept of Riemann sums.

Riemann sums¶

Let us recall their formal definition:

The previous analytic statement may nevertheless hide its geometric intuition: We want to approximate the area under a curve, so we pick $n$ distinct points on $[a,b]$, compute the area of the rectangle with the value of the function as height and the difference between the points as width, and then take their sum. A plot will explain it better than 100 paragraphs of my prose.

Riemann sums are constructed as part of the apparatus needed to define the Riemann integral⁵, so although they are asymptotically optimal, in practice their convergence is very slow, since we are approximating $f$ on an interval by only one of its values; we can do much better.

Every function is (approximately) a polynomial¶

Nowadays, it is very hard to imagine life without a computer, especially when doing quantitative sciences; it would be insufferable to the majority of the readers (also including the author) to numerically estimate an integral by hand, or multiply two very long decimal numbers, or anything along those lines. Nonetheless, if you go through the notes of old geniuses such as Newton, you will find immense patience to perform these kinds of computations, and those endeavors in the long term ended up leading to new insights on how to do the calculations faster; the focus of this section are the Newton-Cotes formulas, one of those brilliant discoveries.

As said before, integrating is hard, but integrating polynomials is something that can be done quite easily; if any complicated function $\mathfrak f$ could be rewritten in this form, the problem would be solved. Unfortunately, this is not possible because, well, it would imply that any function is a combination of polynomials (and hence a polynomial itself), which is certainly not the case; however, we can approximate $\mathfrak f$ by a function of this kind $\tilde{\mathfrak f}$ to get an estimate of the value of the integral, and this is the fundamental idea behind these formulas.

What are our choices? The Taylor polynomial has appeared before, but it has two significant downsides:

1. It is expensive to compute, because it requires knowing (or estimating, which takes even more resources) $n$ derivatives of $\mathfrak f$ at the center in order to find the Taylor polynomial of degree $n$.

2. It is not really that good; you can look at the plot below and convince yourself that the error is quite noticeable once you move away from the center, and to explain this behavior you may recall that among the properties of the Taylor polynomial were included equality on the first, second, and so forth up until the $n$^th derivative, which is a constraint that caps its precision.

How about an interpolating polynomial? This looks much more promising, and it was Newton who developed a method to interpolate $n$ points rapidly (in $\mathscr O(n^2)$ time complexity, using the notation from the first part) and with very high numerical stability; the way it works is as follows:

Newton interpolation¶

Let $x_1, x_2, \cdots, x_n$ be the (distinct) abscissae of the points to interpolate, and $\mathfrak f$ our complicated function. Instead of using the canonical basis for degree-$n$ polynomials ($\{1, x, x^2, \ldots, x^n\}$), we will choose the more clever basis $\{1,\ (x-x_1),\ (x-x_1)(x-x_2),\ldots,\ (x-x_1)\cdots (x-x_{n-1})\}$; call them $\phi_1,\ldots,\phi_n$ respectively. For $p = \lambda_1 \phi_1 + \cdots + \lambda_n \phi_n$ to interpolate the $\mathfrak f(x)$'s, it must satisfy $$ \begin{cases} p(x_1) &= \mathfrak f(x_1),\\ p(x_2) &= \mathfrak f(x_2),\\ &\vdots\\ p(x_n) &= \mathfrak f(x_n), \end{cases} \implies \begin{cases} {\color{red} \lambda_1} &= \mathfrak f(x_1),\\ {\color{red} \lambda_1} + {\color{blue} \lambda_2} (x_2 - x_1) &= \mathfrak f(x_2),\\ &\vdots\\ {\color{red} \lambda_1} + {\color{blue} \lambda_2} (x_n - x_1) + \cdots + {\color{orange} \lambda_n} \prod_{i=1}^{n-1} (x_n - x_i) &= \mathfrak f(x_n), \end{cases} $$ which is a lower-triangular system of equations for the $\lambda$'s, and it can be solved in the claimed quadratic time complexity via forward substitution⁶. There is another route (which was the main discovery of Newton) and, although asymptotically equal, in practice it is faster because it does not require computing each $(x_i-x_1)\cdots(x_i-x_{i-1})$; it is also more stable (less prone to floating-point errors) and it does not require allocating $n\times n$ floats to store the matrix.

The method is called divided differences, and I like to think about it as a discrete version of the Taylor polynomial where, instead of having derivatives, we have approximations that make it equal to the function at the given points. The procedure involves constructing a table⁷ with entries below the diagonal, whose first two columns are the abscissae and the values of the function at those points, and the entry at row $i$ column $j$, called $\nabla^{j-1}_{i-j+1}\mathfrak f$ is computed recursively as $$ \nabla^j_i\mathfrak f = \frac{\nabla^{j-1}_{i+1}\mathfrak f - \nabla^{j-1}_i\mathfrak f} {x_{i+j-1} - x_i}, $$ which is sort of a $(j-1)$^th derivative of $\mathfrak f$ at $x_{i+j-1}$. My notation is arguably kind of cumbersome, but the algorithm is as simple as subtracting the element at its left, and the one at its up and left, and dividing the result by the difference of $x_i$'s that are found when moving always left from the current cell, minus that one found when moving left and up until finding an $\mathfrak f(x_i)$, then the cell to its left. Finally, the elements in the diagonal are the seeked $\lambda$'s for the interpolating polynomial.

Probably the most clever fact about Newton interpolation is that it allows for a very good estimation of the error. Notice that for a point $x_{n+1}$ that is not one of the given nodes, we can write $\mathfrak f(x_{n+1})$ as $p(x_{n+1}) + \varepsilon$, but also as the Newton interpolant with $x_{n+1}$ as node, which by the previously described recurrence formula can be written in terms of $p$ with a correction factor, yielding $$ \mathfrak f(x_{n+1}) = p(x_{n+1}) + \underbrace{\nabla_1^n\mathfrak f \prod_{i=1}^n (x_{n+1} - x_i)}_\varepsilon. $$ The reason why this equation is any useful is that the intuition about the divided differences being related to the derivatives of the function is not only a perception but an actual result. Having $a := \min x_i, b := \max x_i$, it was proven that $$ \text{exists $\xi \in [a, b]$ such that}\qquad \nabla_1^n\mathfrak f = \frac{\mathfrak f^{(n)}(\xi)}{n!},\tag{1} $$ just like in the continuous case with the Taylor polynomial. Hence, at least we can bound the magnitude of the error, as⁸ $$ |\varepsilon_n(x)| \leq \left(\sup_{\zeta\in[a,b]} \frac{|\mathfrak f^{(n)}(\zeta)|}{n!}\right) \left(\sup_{\xi\in[a,b]} \prod_{i=1}^n |\xi - x_i|\right) \leq \lVert f^{(n)}\rVert_\infty (b - a)^n. $$

Trapezoid and Simpson's rule¶

We have wandered along a very long tangent, but the tools we have acquired along the way will pay off here. As it was stated a few paragraphs above, it is a good idea to replace $\mathfrak f$ by its interopolating polynomial and take the integral of the latter instead; the Newton-Cotes formulas are the resulting approximations of the integral when this is done with equally-spaced nodes (including the endpoints). The trapezoid rule is the simplest of them, where we use only a first degree polynomial, so we should perhaps start by studying this scenario.

Using Newton's procedure yields $p(x) = \mathfrak f(a) + \dfrac{\mathfrak f(b)-\mathfrak f(a)}{b - a}x$, and taking the integral gives $$ \int_a^b \mathfrak f(x)\diff x \approx \int_a^b \left[\mathfrak f(a) + \frac{\mathfrak f(b)-\mathfrak f(a)}{b-a}x\right] = \frac{b-a}2(\mathfrak f(a) + \mathfrak f(b)); $$ or geometrically, the area bounded by the line interpolant and the boundary points is a trapezoid, whose area is well known to be $h(a+b)/2$, and produces the same result:

This appears to be a very poor approximation, but this view is due to the fact that the interval along which we are integrating is very large, so $\mathfrak f$ may change heavily between the endpoints; as it will be shown later, the smaller the length of the interval, the better the trapezoid rule represents the original integral. In real life, however, this is almost never the case, but it suffices to divide the original integration region into small pieces to which the trapezoid rule can be applied accurately. This divide and conquer approach is called the composite trapezoid rule.

Next comes Simpson's rule¹⁰, named after English mathematician Thomas Simpson, which involves a quadratic interpolating polynomial this time using also the midpoint of $a$ and $b$ as a node. Proceeding as before, we obtain $p(x)=\mathfrak f(a)+2\dfrac{\mathfrak f(m)-\mathfrak f(a)}{b-a}(x-a) +2\dfrac{\mathfrak f(b)-2\mathfrak f(m)+\mathfrak f(a)}{(b-a)^2} (x-a)(x-m)$, with $m = (a+b)/2$; therefore $$ \int_a^b \mathfrak f(x)\diff x \approx \int_a^b \left[\mathfrak f(a) + 2\frac{\mathfrak f(m)-\mathfrak f(a)}{b-a}(x-a) + 2\frac{\mathfrak f(b)-2\mathfrak f(m)+\mathfrak f(a)}{(b-a)^2} (x-a)(x-m)\right]\diff x = \frac{b-a}{6}\left[\mathfrak f(a) + 4\mathfrak f\left(\frac{a+b}2\right) + \mathfrak f(b)\right]. $$ The plot in this case is very similar to the previous one:

And mixing the previous two gives raise to the composite Simpson's rule.

As a prudent scientist should always do, we need to understand the error of our simplified model, but we are in the right position to ask such a question since in the previous section we studied the error committed by the interpolating polynomial. Thus, integrating the expression from equation $(1)$ and invoking the integral mean value theorem gives back some useful information: $$ {\rm E} = \int_a^b \varepsilon_n(x) \diff x = \int_a^b \left(\frac{\mathfrak f^{(n)}(\phi(x))}{n!} \prod_{k=0}^{n-1} (x-(a+hk))\right) \diff x = \frac{\mathfrak f^{(n)}(\xi)}{n!} \int_a^b \left(\prod_{k=0}^{n-1} (x-(a+hk))\right) \diff x. $$ For instance, when $n=2$ (trapezoidal rule) this reduces down to $-h^3\mathfrak f''(\xi)/12$, and when $n=3$ it can be checked to be $-h^5\mathfrak f^{(\text{iv})}(\xi)/90$; notice that the index of the derivative of $f$ always seems to be an even integer, but we will discuss this interesting property in a moment.

The error of the composite formulas follows immediately from their definition and the previous result; for example, in the case of Simpson's rule: $$ \bar{\rm E} = \sum_{i = 0}^{n - 1} \int\limits_{a + ih}^{a + (i+1)h} \varepsilon_3(x)\diff x = \sum_{i = 0}^{n - 1} \left(-\frac{\mathfrak f^{(\text{iv})}(\xi)}{90}h^5\right) = -\frac{\mathfrak f^{(\text{iv})}(\xi')}{90} h^5 n = -\frac{\mathfrak f^{(\text{iv})}(\xi')}{90} h^4 (b - a) = \mathscr O(h^4), $$ and naturally we had to pay a fee for splitting the integral, that being that the order of convergence (the exponent of the $h$) dropped from 5 to 4; contrary to regular algorithm complexity, the higher the exponent the better, since $h$ is getting smaller and smaller, approaching 0.

The unreasonable effectiveness of Simpson's rule¶

We have seen that the error for the $n$^th Newton-Cotes rule depends at least on the $n$^th derivative of the function, so polynomials of degree up to $n$ are not only approximated, but integrated precisely; the largest degree of a polynomial that is exactly integrated by the method is called its degree of precision. It makes sense that the trapezoid rule is correct for lines, because it uses an interpolant that is equal to the function, and likewise for Simpson with parabolas; however, notice that the error term of the latter depended on $\mathfrak f^{(\text{iv})}$, which vanishes even for cubics, so in theory the degree of precision of Simpson's rule is 4, not 3. Why is that so?

When I first learned this result I believed it was quite nice, and that is the reason I have included it in this addendum:

The crux of the proof was that quadratic interpolants share three points with the original function, and that the points were evenly spaced; in simple terms, because of the last fact the error is distributed as an odd function around the midpoint of the interval, which is annihilated when taking the integral.

Notice, though, that the previous properties are not exclusive of second-degree interpolants, but can be generalized to any even ones, so even-numbered Newton-Cotes rules have an extra degree of precision. Is it not beautiful? Try some examples in the console above to help convince yourself!

Richardson extrapolation and Romberg's method¶

There is one last concept to cover before diving into the integration method of the old HP calculators, namely Richardson extrapolation as an acceleration method, and Romberg's method as its application to integrals. So far we have studied the order of convergence of our methods (that being how fast the error decreases when $h$ shrinks to 0), and although it is not bad for Simpson's rule and the like, it would be good to improve it in some way.

Richardson extrapolation achieves precisely that but in a more general fashion, for any numerical method rather than focusing on integrals. Assume we have a quantity $Q$, which can be computed by a method $\hat Q(h)$ that has order of convergence $n$ (that is, $\hat Q(h) = Q + \mathscr O(h^n)$); we can improve the $n$ to an $n+1$ by doing $$ \tilde Q(h) = \frac{\lambda^n\hat Q(h/\lambda) - \hat Q(h)}{\lambda^n - 1} \qquad\text{for any $\lambda \in \RR$}, \tag2 $$ which works because the $h^n$ cancel out, as it can be seen below: $$ \tilde Q(h) = \frac {\lambda^n(Q + \cancel{th^n \lambda^{-n}} + \mathscr O(h^{n+1})) - (Q + \cancel{th^n} + \mathscr O(h^{n+1}))} {\lambda^n - 1} = \frac{(\lambda^n-1)Q + \mathscr O(h^{n+1})}{\lambda^n - 1} = Q + \mathscr O(h^{n+1}). $$

But we are pretty much only interested in computing integrals, so let us apply it for this purpose. The idea is to build a lower-triangular table, where the first column will represent the estimates using a composite rule with $h_n = (b-a)/2^n$ and the rest will be computed from Richardson extrapolation by following $(2)$ with $\lambda=4$. Calling the integral at the $i$^th row and $j$^th column $R_{ij}$, we see from the description that $$ \begin{align*} R_{i1} &= \texttt{composite\_integral}\left(f, a, b, \frac{b-a}{2^i}\right),\\ R_{ij} &= \frac{4^{j-1} R_{i,j-1} - R_{i-1,j-1}}{4^{j-1} - 1}. \end{align*} $$

As always, regarding the error, it can be checked that $R_{ij}$ has a fabulous $\mathscr O(h_n^{2m})$, so this is a very noticeable speedup; evidently, the final approximation is stored at the last diagonal entry.

Kahan's brilliant idea¶

We have learned a lot today, especially on how to find integrals efficientely, so why did HP calculators use a different method if Romberg-Simpson works perfectly fine? So far, we have only considered functions that are well-behaved in the interval of integration, but having a singularity at one of the endpoints does not mean that the integral does not converge. For instance, $$ \int_0^1 \frac{\diff x}{\sqrt x} = \lim_{x \to 0^+}\int_x^1 \frac{\diff x}{\sqrt x} = \lim_{x \to 0^+} (2\sqrt x)\bigg|_x^1 = 2, $$ yet typing composite_trapezoid(x => 1/Math.sqrt(x), 0, 1) in the console returns null, meaning that a division by zero has occurred.

This is not the behavior you would expect from a calculator, and the brilliant idea William Velvel Kahan came up with in order to circumvent this issue was to apply the following cubic substitution: $$ \int_a^b \mathfrak f(x) \diff x = \left\{\begin{aligned} x &= \frac{a+b}2 + \frac{b-a}4 y (3 - y^2),\\ \diff x &= \frac34(b-a)(1-y^2) \diff y \end{aligned}\right\} = \int_{-1}^1 \mathfrak f\underbrace{\left(\frac{a+b}2 + \frac{b-a}4 y (3 - y^2)\right)}_{x(y)} \cdot \underbrace{\frac34(b-a)(1-y^2)}_{\diff x/\diff y} \diff y, $$ which has two marvelous properties:

1. The integration region is always $[-1, 1]$.
2. More importantly, $\diff x/\diff y$ vanishes at the endpoints, so those need not be computed, avoiding any possible singularity.

Afterwards, our beloved trapezoid rule with Romberg's method does the trick, always being careful to return 0 as soon as we verify that the Jacobian cancels out.

References¶

—. Michel Jean, "Romberg–Kahan Integrator," AxiomeLab, Tech. Rep., 2026. url: https://axiomelab.gitlab.io/hp34c/.

—. Michel Jean, "The Secret Behind the ∫ Key on Your HP Calculator," AxiomeLab, Tech. Rep., 2026. url: https://axiomelab.gitlab.io/hp-integral/article_hp_en.html.

William Kahan. "Handheld calculator evaluates integrals," in Hewlett-Packard Journal, vol. 31, 1980.

A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics. Springer Berlin Heidelberg, 2007. doi: 10.1007/b98885.

¹ As long as the function is integrable it is never truly impossible, since for a function $\mathfrak f$ of this kind, one could simply define $\mathfrak F$ to be its antiderivative; however, this is useless unless it leads to an easy way of computing $\mathfrak F(x)$.

² In computer science, a word is the smallest unit of data handled in the instruction set of a processor.

³ In contrast to floating-point numbers, where the separator between the integer and fractional part can be "moved" (choosing a different exponent), fixed-point numbers have a fixed number of bits allocated for each (integer and fractional) part.

⁴ Formally, we say $f = \Omega(g)$ if $g = \mathscr O(f)$, and $f = \Theta(g)$ if $f = \mathscr O(g)$ and $g = \mathscr O(f)$.

⁵ The Riemann integral is essentially the limit of the Riemann sums as we let $n$ go to infinity. In a real analysis textbook, however, the upper and lower are first defined taking $f(x_i^*) = \sup f([x_{i-1}, x_i])$ and $f(x_i^*) = \inf f([x_{i-1}, x_i])$ respectively (that is, they are the upper and lower bounds of any Riemann sum), and if the limit as $n$ approaches infinity of the two is the same quantity $L$, we say that $f$ is Riemann-integrable on $[a, b]$, and write $\int_a^b f(x) \diff x = L$.

⁶ From the first equation one obtains the value of $\lambda_1$, then it is substituted into the second one to find $\lambda_2$, and so forth until $\lambda_n$. This procedure (generalized) is named forward substitution, and the reader might benefit from thinking about why it requires $n(n-1)$ floating-point operations. There is also a backward substitution, which is identical but for upper-triangular systems rather than lower-triangular ones.

⁷ In practice there is no need to build the full table, since the only important values are the diagonal ones; it suffices to allocate $2n$ floats to store the first and second columns, and then update the latter at every iteration. Since the third column only has $n-1$ elements, the remaining space can be used for the first diagonal entry, and so forth for the rest of the table.

⁸ For those not familiar with the concept of $p$-norms, the particular case of the norm infinity ($\lVert\cdot\rVert_\infty$) can be thought of as the largest value of a sequence (formally, of a vector from the space the norm acts upon) in absolute values.

⁹ console.log(solution1) will only reveal the numerical value of the integral rather than a closed-form expression that can be understood. From my perspective, the easiest way to solve the integral is by substituting $x$ into a hyperbolic cosine to simplify the square root in the denominator; hence $$ \int_2^4 \frac{\diff x}{\sqrt{x^2 - 1}} = \left\{\begin{aligned} x &= \cosh t\\ \diff x &= \sinh t \diff t \end{aligned}\right\} = \int\limits_{\arccosh 2}^{\arccosh 4} \frac{\sinh t}{\sqrt{\cosh^2 t - 1}} \diff t = \int\limits_{\arccosh 2}^{\arccosh 4} \cancel{\frac{\sinh t}{\sinh t}} \diff t = {\color{green} \boxed{\arccosh 4 - \arccosh 2.}} $$

¹⁰ In reality, there are two rules by Simpson; they are the Newton-Cotes formulas for second and third degree interpolants, and they receive the names $1/3$- and $3/8$-Simpson's rule when it is necessary to distinguish between them (which is almost never the case, and this article will not be an exception).

¹¹ Once again, solution2 only contains the numerical answer, but a symbolic one can be achieved through integration by parts. Indeed, notice that $$ \int \log x \diff x = \left\{\begin{aligned} u &= \log x\\ \diff v &= \diff x \end{aligned}\right\} = x\log x - \int \cancel{\frac xx} \diff x = x(\log x - 1), $$ and through a simple change of variables $$ \int_{-2}^1 \log(x + 2) \diff x = \lim_{t \to 0^+} \int_t^3 \log x \diff x = \lim_{t \to 0^+} \left[x(\log x - 1)\right]_t^3 = 3(\log 3 - 1) - \underbrace{\lim_{t \to 0^+} t(\log t - 1)}_{\rm L}, $$ but we can solve $\rm L$ by rewriting it conveniently and using L'Hôpital's rule: $$ {\rm L} = \lim_{t \to 0^+} t(\log t - 1) = \lim_{t \to 0^+} \frac{\log t - 1}{1/t} \stackrel{\text{L'H}}= \lim_{t \to 0^+} \frac{1/t}{-1/t^2} = 0; $$ putting everything together, the result turns out to be $$ \int_{-2}^1 \log(x + 2) \diff x = 3(\log 3 - 1) - {\rm L} = {\color{green}\boxed{3(\log 3 - 1).}} $$ Also, I write $\log x = \ln x \neq \log_{10} x$; I have not yet found a convincing argument to have $\log x = \log_{10} x$.