An aggregation is some kind of summary of a set of data. This can be the sum, length, minimum, etc. It is quite common to want to calculate such a summary repeatedly, e.g. “the maximum noise level in dB for the past 30 seconds” for a nuisance detector. In such a case we say there is a sliding window over our data, and we want to aggregate over our window. If our aggregation is a binary operator with an inverse, like integer sums, there is a very easy solution using a double-ended queue: from collections import deque class SlidingWindowSum: def __init__(self): self.sum = 0 self.elems = deque() def push(self, x): self.sum += x self.elems.append(x) def pop(self): self.sum -= self.elems.popleft() def eval(self): return self.sum But what if our operator has no inverse? This is actually the case for most interesting summaries such as minimum, quantile, approximate unique count (for example using HyperLogLog), etc. In fact, even something as simple as a floating-point sum suffers from the fact that floating-point addition is not invertible. For example, if you ever have a NaN in your input data with the above naive algorithm your sum will forever remain NaN, even long after the bad value has left your window. Six years ago I came up with an algorithm for maintaining just the minimum/maximum in a sliding window and posted it to cs.stackexchange. I now consider this algorithm pointless, because it turns out there is a simple and efficient algorithm that solves this problem for a very wide class of aggregations. I’m writing this blog post to spread the word, because I feel it should be more widely known. Folklore I came across this algorithm while reading a far more advanced paper, Low-Latency Sliding-Window Aggregation in Worst-Case Constant Time by Tangwongsan et al. Why is this paper titled low-latency? Because it does the same as what I’m about to describe, but in O(1) time for each step. However, in it they also described a “two-stack” algorithm, which does it in amortized O(1), and is far, far simpler. Amortized O(1) means that across many operations the total amount of work per element is constant, but an individual operation can take much longer. This is almost always fine, unless you absolutely need a low upper bound on latency. Funnily enough that paper attributes this algorithm to “adamax” from a 2011 Stack Overflow post. They in turn credit a 2001 lecture note by D. Sleator for the inspiration. However, this lecture note does not describe a sliding window aggregate, it describes the classical two-stack algorithm for implementing a FIFO queue and does amortized analysis on it. Ultimately I would not be surprised to find that this algorithm was already described in an obscure paper from the 1970s, seeing how simple and brilliant it is. Two stacks Like the authors of the paper, I will generalize the two-stack algorithm to arbitrary associative aggregation functions. By abstracting the aggregation as a set of functions, empty(), unit(x), combine(x, y) and finalize(x), you can describe many possible aggregations, for example a mean: empty = lambda: (0, 0) unit = lambda x: (x, 1) combine = lambda x, y: (x[0] + y[0], x[1] + y[1]) finalize = lambda x: x[0] / x[1] if x[1] else None I’d like to note here that these functions have the following signatures: fn empty() -> Agg; fn unit(x: Value) -> Agg; fn combine(x: Agg, y: Agg) -> Agg; fn finalize(x: Agg) -> Out; I’m making a distinction here between Value, Agg and Out because while they seem superficially similar for something like an integer sum, for an approximate unique count on strings you would have (Value, Agg, Out) = (String, HyperLogLogSketch, u64), three wildly different types. Without further ado, the algorithm: class TwoStackAgg: def __init__(self): self.values = [] self.values_agg = empty() self.cum_aggs = [] def push(self, x): self.values.append(x) self.values_agg = combine(self.values_agg, unit(x)) def pop(self): if not self.cum_aggs: cum_agg = empty() while self.values: cum_agg = combine(unit(self.values.pop()), cum_agg) self.cum_aggs.append(cum_agg) self.values_agg = empty() self.cum_aggs.pop() def eval(self): return finalize( combine(self.cum_aggs[-1], self.values_agg) if self.cum_aggs else self.values_agg ) That’s it, the entire algorithm. There’s two stacks (values and cum_aggs) and one more aggregate, values_agg. At any point in time values_agg holds the aggregate of values, and cum_aggs contains the cumulative aggregates of all values in our window that aren’t in values, in reverse order. From this we can get the aggregate over our entire window in constant time by by combining the last value of cum_aggs with values_agg. The neat part is that (assuming w is our window size) every wth operation we drain all of values and maintain a running aggregate while pushing the partial cumulative aggregates onto cum_aggs. This is what makes it amortized O(1), doing O(w) internal operations every wth pop bounds the total amount of work per element to O(1), even though a singular operation might not be constant time. I think this is best visualized. Suppose we sum [1, 2, ..., 10] with a fixed-size sliding window of four elements, then the state on each eval() call would look like this (values_agg not shown as it is simply the aggregate of the values): cum_aggs values out [] [] = 0 [] [1] = 1 [] [1, 2] = 1 + 2 [] [1, 2, 3] = 1 + 2 + 3 [] [1, 2, 3, 4] = 1 + 2 + 3 + 4 [4, 3 + 4, 2 + 3 + 4] [5] = 2 + 3 + 4 + 5 [4, 3 + 4] [5, 6] = 3 + 4 + 5 + 6 [4] [5, 6, 7] = 4 + 5 + 6 + 7 [] [5, 6, 7, 8] = 5 + 6 + 7 + 8 [8, 7 + 8, 6 + 7 + 8] [9] = 6 + 7 + 8 + 9 [8, 7 + 8] [9, 10] = 7 + 8 + 9 + 10 [8] [9, 10] = 8 + 9 + 10 [] [9, 10] = 9 + 10 [10] [] = 10 [] [] = 0 In total the memory usage is O(w), where w is your maximum window size. Note that for simplicity of analysis and the example I assumed a fixed-size window w, but there is nothing about the two-stack algorithm that requires this. You can call push(x) and pop() as many times as you’d like between each eval(), growing and shrinking the window size as needed. Floating-point non-associativity Note that we required above that our aggregate combine is associative, meaning: combine(combine(x, y), z) = combine(x, combine(y, z)) Technically speaking, floating-point addition doesn’t respect this. Nevertheless, the above algorithm is still very useful because the results closely match the expected outcome, even more so if you use a compensated summation algorithm like Kahan summation. Another neat thing about the two-stack algorithm is that it doesn’t require commutativity, if you follow the above implementation precisely. The order of operands is maintained, which can matter for things like string concatenation. However, there is a second very useful property of the above algorithm. Each aggregate is strictly a combination of the elements in the window, and none outside the window. This means if your window contains a NaN or infinity (or some other outlier), that value only poisons the windows that contain it rather than the rest of your computation. But even without NaN or infinity it is useful, due to not propagating errors endlessly. E.g. if your sliding window starts with [1e20, 1], this is what would happen with a naive rolling sum: >>> 1e20 + 1 - 1e20 - 1 -1.0 Compensated summation will reduce these effects, but not making your result depend on values outside of the window will eliminate long-term error accumulation entirely.

The following incredibly small sorting algorithm has an $O(n^{4/3})$ worst-case runtime: def fibonacci_sort(v): a, b = 1, 1 while a * b < len(v): a, b = b, a + b while a > 0: a, b = b - a, a g = a * b for i in range(g, len(v)): while i >= g and v[i - g] > v[i]: v[i], v[i - g] = v[i - g], v[i] i -= g As the name implies, it uses the Fibonacci sequence ($1, 1, 2, 3, 5, \dots$) to sort the elements. In this article I will explain how it works, show an interesting divisibility property of the Fibonacci numbers and use that to prove its complexity. I will also explain how I ended up with one of Donald Knuth’s coveted reward checks. Shellsort The above sorting algorithm is an instance of the more generic Shellsort. Shellsort, named after Donald L. Shell, does not directly sort all the elements in one step. Rather, it first only sorts subsequences of the array in a process known as $k$-sorting. $k$-sorting A subsequence of an array can be any of its elements, but they must remain in the original order. Unlike a substring, the subsequence can have gaps—the elements need not be contiguous. When $k$-sorting an array you split the elements of the array up in groups, and sort those groups independently from each other. Each group is formed by taking every $k$th value ($k$ is also referred to as the gap), differing only by their starting point. For example, $3$-sorting an array with eight elements looks like this: After this step we say the array is $k$-ordered. This means that for all $i$ we have ${A[i] \leq A[i + k]}$. Nothing is (directly) known about the relative order of elements which aren’t separated by a multiple of $k$. Since everything is a multiple of $1$, we see that $1$-ordered is just… sorted. Confusingly, there appears to be another definition for $k$-sorted which states that all $i$ and all $j \geq k$ we have ${A[i] \leq A[i + j]}$, not just at the multiples of $k$. That is not the case here, in the context of Shellsort all literature uses the earlier definition. There is a fascinating property of $k$-sorting that is key to how Shellsort operates: Theorem 1. If an array is $h$-ordered and then it is $k$-sorted, it remains $h$-ordered as well. It is surprisingly tricky to prove this simple statement. A proof sketch paraphrased from a formal proof due to Knuth (TAOCP, Vol. 3, Chapter 5.2.1, Theorem K) goes as follows: Lemma 1. If the last $r$ elements from array $Y$ are bigger or equal to respectively the first $r$ elements from array $X$, then this remains true after sorting both arrays. Proof sketch. There are at least $r$ elements in $X$ which are smaller than or equal to elements in $Y$, thus the maximum element in $Y$ dominates at least $r$ elements, and thus the last element of $Y$ dominates the $r$th element of $X$ after sorting both. Apply a similar argument for $r - 1$ and the second largest element, etc. A visual representation of this lemma really helps the understanding I believe: Now we can use this lemma to prove our original theorem: Proof sketch of Theorem 1. We have some array $A$ which is $h$-ordered, and thus we have $A[i] \leq A[i + h]$ for valid $i$. Then we $k$-sort it and now have $A[i] \leq A[i + k]$ instead. For any particular choice of $i$, define $X$ as $A[i + sk]$ for all valid integer $s$, and $Y$ as $A[i + tk + h]$ for all valid integer $t$. Then you can apply Lemma 1 to show that $A[i] \leq A[i + h]$ still remains true after $k$-sorting. Once again I think a visual representation really helps here, especially to see the parallels with Lemma 1. Suppose we have a $3$-ordered array and then $5$-sort it. The proof sketch of Theorem 1 for $i = 2$ (and in fact for any $i \equiv 2 \pmod 5$) that we still have $A[i] \leq A[i + 3]$ can then be visualized as such: I chose to only show a proof sketch rather than a full formal proof in this blog post, as I quote: “This is much harder to write down than to understand.” — Donald Knuth The result of this theorem is that as you apply more and more steps of $k$-sorting, all the work compounds into a larger overall ordering. Even with just two steps we can see very powerful results: Lemma 2. If an array is both $h$-ordered and $k$-ordered where $k$ and $h$ are relatively prime (they don’t share any divisors other than 1), then two elements in the array which are at least $s \geq (h-1)(k-1)$ steps apart are in a correct relative order. Proof. It is possible to write $s = \alpha h + \beta k$ with integer $\alpha, \beta \geq 0$ (for a proof of that see here). This means we can do $\alpha$ steps of $h$ (each time using $A[i] \leq A[i + h]$ since we’re $h$-ordered) and similarly $\beta$ steps of $k$ to see that $A[i] \leq A[i + \alpha h + \beta k] \leq A[i + s]$. For example, here is a visualization of the relative order implied by the combination of $4$-ordering and $9$-ordering: We see that indeed starting from offset $(4 - 1)(9 - 1) = 24$ each element is ordered relative to the element at 0. More generally, we find that if we $k$-sort with some set $\{k_i, \dots, k_j\}$ and $\gcd(k_i, \dots, k_j) = 1$, then there exists some upper bound such that all numbers greater than it are representable as sums of non-negative multiples of $k_i, \dots, k_j$, meaning all elements beyond that point are ordered correctly relative to the starting point. Finding this upper bound is known as the Frobenius problem or coin problem, and it is a hard number theoretic problem. Some formulae like the above one for two coprime integers are known but the general problem for arbitrary sets of $k$ is NP-hard. It is this problem that forms the surprising link between Shellsort and number theory which will ultimately lead us towards the Fibonacci numbers. Insertion sort Shellsort uses insertion sort to $k$-sort the subsequences. This sort repeatedly swaps the last unsorted element with the element before it until it falls into place before continuing with the next unsorted element. This is generally speaking a rather inefficient algorithm for large arrays, as it has a worst-case of $O(n^2)$. However, this worst-case complexity can be refined. If each element is at most $m$ steps away from its final sorted position, insertion sort takes $O(nm)$ time. And this is the key to Shellsort’s subquadratic worst-case complexity. You choose a clever gap sequence such that gaps start off very large leading to small subsequences and thus small $O(n^2)$ terms. Then, when Shellsort starts $k$-sorting with smaller gaps, you can use the fact that no element is very far from its final position to prove that it is still efficient. For example, based on our earlier observations, if the array is $4$-ordered and $9$-ordered and we do a $1$-sort (that is, a regular sort with no gaps) we know insertion sort can not run slower than $O(24n) = O(n)$ because each element is no further than $m = 24$ steps away from its final position. Gap sequence Choosing the gap sequence is thus key to Shellsort’s performance. The Wikipedia page lists many known gap sequences with different complexities. For example Hibbard’s 1963 sequence $k_i = 2^{i} - 1$ with $O(n^{3/2})$ complexity, or Pratt’s from 1971 which chooses all numbers of the form $2^p3^q$ for a complexity of $O(n\,(\log n)^2)$ (still the best known sequence to date, at least asymptotically). However, all ‘new’ sequences listed after 1986 are empirically established, their complexities are unknown. They perform very well in practice, but they could possibly have much slower than expected performance for some inputs. A search on Google Scholar reveals no interesting new sequences either. With that in mind I’m happy to announce that Shellsort with the gap sequence $$k_i = F_i \cdot F_{i+1} = (1, 2, 6, 15, 40, 104, 273, \dots)$$ where $F_n$ is the $n$th Fibonacci number has a worst-case runtime of $O(n^{4/3})$. An interesting Fibonacci property Before we can prove this worst-case runtime of gap sequence we’re going to need to take a look at the Fibonacci numbers in more detail: $$F_0 = 0, \quad F_1 = 1, \quad F_n = F_{n-1} + F_{n-2}$$ This very well-known series of numbers $0, 1, 1, 2, 3, 5, 8, 13, \dots$ has all kind of interesting properties, but today we’ll need a very curious property where for $a, b > 0$, $$\gcd(F_a, F_b) = F_{\gcd(a, b)},$$ where $\gcd$ is the greatest common divisor function. This property in plain words states that the greatest common divisor of the $a$th and $b$th Fibonacci number can be found by computing the greatest common divisor of $a$ and $b$, and looking up that index in the Fibonacci series. Addition formula To prove that we first need to prove the addition formula of the Fibonacci numbers, where $n, m > 0$, $$F_{n+m} = F_{n+1}F_m + F_nF_{m-1}.$$ A fairly simple way to do this is using the matrix exponential form of the Fibonacci sequence, $$\begin{pmatrix}F_{n+1}&F_n\\F_n&F_{n-1}\end{pmatrix} = {\begin{pmatrix}1&1\\1&0\end{pmatrix}}^n.$$ This form can be easily proven as correct using induction. With $n = 1$ you can directly see the equality holds, and to prove it holds for $n + 1$ assuming it holds for $n$ we have \begin{align*} {\begin{pmatrix}1&1\\1&0\end{pmatrix}}^{n + 1} &= {\begin{pmatrix}1&1\\1&0\end{pmatrix}}{\begin{pmatrix}1&1\\1&0\end{pmatrix}}^n = {\begin{pmatrix}1&1\\1&0\end{pmatrix}}\begin{pmatrix}F_{n+1}&F_n\\F_n&F_{n-1}\end{pmatrix}\\ &= \begin{pmatrix}F_{n+1}+F_{n}&F_n+F_{n-1}\\F_{n+1}&F_{n}\end{pmatrix} = \begin{pmatrix}F_{n+2}&F_{n+1}\\F_{n+1}&F_{n}\end{pmatrix}. \end{align*} Then, using the fact that in matrix exponentiation $A^{n+m} = A^n \times A^m$ we find: \begin{align*} \begin{pmatrix}F_{n+m+1}&F_{n+m}\\F_{n+m}&F_{n+m-1}\end{pmatrix} &= \begin{pmatrix}F_{n+1}&F_n\\F_n&F_{n-1}\end{pmatrix} \times \begin{pmatrix}F_{m+1}&F_m\\F_m&F_{m-1}\end{pmatrix}\\ &= \begin{pmatrix}F_{n+1}F_{m+1}+F_nF_m&F_{n+1}F_m+F_nF_{m-1}\\F_nF_{m+1}+F_{n-1}F_m&F_nF_m + F_{n-1}F_{m-1}\end{pmatrix}. \end{align*} Our desired addition formula can then be read from the top right element of both sides of the matrix equation. If we generalize the Fibonacci formula a bit to allow negative indices we find that $F_{-1} = 1$ (which is the only choice preserving $F_{n} = F_{n-1} + F_{n-2}$), and that the above addition formula also holds for $n > 0, m = 0$ which we’ll need in a bit. Divisibility and coprimality We also need to prove that $F_{kn} \equiv 0 \bmod F_{n}$ for all $k, n > 0$. If we do induction on $k$ it obviously holds true for $k = 1$, and we see that if it holds for $k$ we have \begin{align*} F_{(k + 1)n} \equiv F_{kn + n} &\equiv F_{kn+1}F_{n} + F_{kn}F_{n - 1}\\ &\equiv F_{kn+1} \cdot 0 + 0 \cdot F_{n - 1} \equiv 0 &&\mod F_{n}. \end{align*} In other words, $F_n$ divides $F_{kn}$, also written $F_{n} \mathrel{|} F_{kn}$. Finally we need to prove that $\gcd(F_n, F_{n+1}) = 1$, meaning consecutive Fibonacci numbers don’t share any prime factors (they are coprime). This is once again proven through induction, the property holds for $n = 1$ and since $\gcd(a, b) = \gcd(a, a + b)$, we have $$\gcd(F_n, F_{n+1}) = \gcd(F_n + F_{n+1} , F_{n+1}) = \gcd(F_{n+1}, F_{n+2}).$$ Euclidean algorithm The Euclidean algorithm is an algorithm to compute the greatest common divisor based on the identity $$\gcd(a, b) = \gcd(a, b - a),$$ where $1 < a < b$. By repeatedly applying this identity (swapping $a, b$ when needed to keep $a < b$) until $a = 1$ you’ll find the greatest common divisor. You can speed up the algorithm by replacing the repeated subtraction with the modulo operator: $$\gcd(a, b) = \gcd(a, b \bmod a).$$ Now, suppose we are computing $\gcd(F_a, F_b)$ with $0 < a < b$. We can write $b = qa + r$ with $q \geq 1$ and $0 \leq r < a$, this is just splitting up $b$ into the quotient $q$ and remainder $r$ when dividing by $a$. Then, using the addition formula we find $$\gcd(F_a, F_b) = \gcd(F_a, F_{qa + r}) = \gcd(F_a, F_{qa + 1}F_r + F_{qa}F_{r-1})$$ We proved earlier that $F_{qa} \equiv 0 \mod F_a$, so applying the above modular identity we can eliminate the $F_{qa}F_{r-1}$ term, $$\gcd(F_a, F_b) = \gcd(F_a, F_{qa + 1}F_r).$$ We also know that $F_{qa+1}$ and $F_{qa}$ are coprime. Since all factors of $F_a$ are factors of $F_{qa}$ we can conclude that $F_a$ and $F_{qa+1}$ share no (prime) factors either. This property lets us eliminate that factor from the greatest common divisor calculation entirely: $$\gcd(F_a, F_b) = \gcd(F_a, F_r).$$ Since $r = b \bmod a$ we note the following symmetry when $1 < a < b$: $$\gcd(a, b) = \gcd(a, b \bmod a),$$ $$\gcd(F_a, F_b) = \gcd(F_a, F_{b \bmod a}).$$ In other words, by repeatedly applying the above identity we can perform the Euclidean algorithm on the indices of the Fibonacci numbers, giving $$\gcd(F_a, F_b) = F_{\gcd(a, b)}.$$ Fibonacci sort’s complexity Now, for the actual proof I have to stand on top of the shoulders of giants. Recall that our gap sequence was as follows: $$k_i = F_{i} \cdot F_{i + 1}.$$ In the 1986 paper A New Upper Bound for Shellsort by Robert Sedgewick, a formula due to Johnson is shown (Theorem 4) which lets us bound the Frobenius term $g(a, b, c)$ for three numbers $a, b, c$, assuming that they are independent. The Frobenius term $g(a, b, c)$ is what we discussed earlier in our analysis of Shellsort. It means that for all $n \geq g(a, b, c)$ there exist some combination of non-negative integers $\alpha, \beta, \gamma$ such that $$\alpha a + \beta b + \gamma c = n.$$ Conversely, independence means none of the numbers can be written as a linear combination of the others with non-negative integer coefficients. Assuming $i \geq 2$, for $\{k_i, k_{i+1}, k_{i+2}\}$ it is easy to prove the independence using our earlier GCD formula. From this we establish that $$\gcd(F_{i}, F_{i+1}) = F_{\gcd(i, i+1)} = F_{1} = 1,$$ $$\gcd(F_{i}, F_{i+2}) = F_{\gcd(i, i+2)} = F_{2} = 1,$$ and thus $$\gcd(k_{i}, k_{i+1}) = \gcd(F_{i}\cdot F_{i+1}, F_{i+1}\cdot F_{i+2}) = F_{i+1}\cdot \gcd(F_{i}, F_{i+2}) = F_{i+1}.$$ This directly proves that $k_{i+1}$ is not a multiple of $k_{i}$. This leaves one other possibility, that $k_{i+2}$ is dependent on $\{k_i, k_{i+1}\}$: $$\alpha k_{i} + \beta k_{i+1} = k_{i+2}$$ $$\alpha F_i F_{i+1} + \beta F_{i+1} F_{i+2} = F_{i+2}F_{i+3}$$ $$F_{i+1}(\alpha F_i + \beta F_{i+2}) = F_{i+2}F_{i+3}$$ This can only be true if $F_{i+2}F_{i+3}$ is a multiple of $F_{i+1}$, but $F_{i+1}$ is coprime with both $F_{i+2}$ and $F_{i+3}$ as we saw earlier, and thus this rules out any dependence. Frobenius bound Since we have proven independence we can use Johnson’s formula from the above paper. It states that $$g(a, b, c) = d \cdot g\left(\frac{a}{d}, \frac{b}{d}, c\right) + (d-1)c$$ where $d = \gcd(a, b)$. This lets us use the formula for coprime $a, b$ we saw in Lemma 2, because when $a, b < c$ and $a, b$ are coprime we have $$g(a, b, c) \leq g(a, b) = (a - 1)(b - 1) \leq ab.$$ In our case this means that $$g(k_{i+1}, k_{i+2}, k_{i+3}) \leq F_{i+2} \cdot g(F_{i+1}, F_{i+3}) + F_{i+2}F_{i+3}F_{i+4},$$ $$g(k_{i+1}, k_{i+2}, k_{i+3}) \leq F_{i+1}F_{i+2}F_{i+3} + F_{i+2}F_{i+3}F_{i+4}.$$ If we now use the asymptotic approximation $F_n = \Theta(\phi^n)$ where ${\phi = (\sqrt{5} + 1)/2}$ we have $$k_i = \Theta(\phi^{2n}), \quad g(k_{i+1}, k_{i+2}, k_{i+3}) \leq \Theta(\phi^{3n}).$$ In other words, $g(k_{i+1}, k_{i+2}, k_{i+3}) = O(k_i^{3/2})$. See my previous blog post on magical Fibonacci formulae to see how the asymptotic approximation is derived. Back to Shellsort We proved just above that after having processed gaps $k_{i+3}, k_{i+2}, k_{i+1}$ the maximum distance for an element from its true location is no more than $O(k_i^{3/2})$. Then if we perform a $k_i$-sort next, we have $k_i$ subsequences of length $n / k_i$ in which the maximum distance of an element from its sorted location within that subsequence is no more than $O(k_i^{3/2} / k_i) = O(k_i^{1/2})$. Recall from earlier that insertion sort on a subsequence of length $n$ where each element is at most $m$ places from its sorted location has complexity $O(nm)$. This means the $k_i$-sort step for all subsequences combined takes at most $O(k_i \cdot (n / k_i) \cdot k_i^{1/2}) = O(n \cdot k_i^{1/2})$ time. This bound is good when $k_i$ is small, but we also have a different upper bound which is good when $k_i$ is big, using the fact that insertion sort takes at most $O(n^2)$ time. This gives us the bound $O(k_i \cdot (n / k_i)^2) = O(n^2 / k_i)$ as well. We note that both bounds are equal when $k_i = \Theta(n^{2/3})$, giving us $O(n^{4/3})$ overall by consistently choosing the smaller of the two bounds. To make this a bit more formal, let $t$ be the maximum index such that $k_t \leq n$, and let $s$ represent the index on which we switch our bound. Then we have $$cT \leq \sum_{i=1}^{s - 1} \left(n \cdot k_i^{1/2}\right) + \sum_{i=s}^t \left(n^2 / k_i\right),$$ where $T$ is our total runtime and $c > 0$ is some ‘constant’ irrelevant for the asymptotic analysis (varying from equation to equation). Then, substituting the asymptotic approximation for $k_i \approx \phi^{2i}$ we get $$cT \leq n \sum_{i=1}^{s - 1} \phi^{i} + n^2\sum_{i=s}^t \phi^{-2i},$$ $$cT \leq n \cdot \frac{\phi^s - \phi}{\phi - 1} + n^2 \cdot \frac{\phi^{2 - 2s} - \phi^{-2t}}{\phi^2 - 1}.$$ Removing constant factors from terms and substituting $s = 2t/3$ gives us $$cT \leq n \cdot (\phi^{2t/3} - 1) + n^2 \cdot (\phi^{- 4t/3} - \phi^{-2t}).$$ Finally we note that within a constant factor $\phi^{2t} \approx k_t \approx n$ to simplify to $$cT \leq n \cdot (n^{1/3} - 1) + n^2 \cdot (n^{- 2/3} - 1/n),$$ $$cT \leq n^{4/3} + n^{4/3} - 2n,$$ giving our overall bound of $O(n^{4/3})$. Generalizing The above proof is analogous to the one found in A New Upper Bound for Shellsort by Sedgewick, I don’t claim credit for it. Ultimately there is nothing fundamental about the choice of the Fibonacci numbers here, any sequence which starts with $1$ and: grows like $\Theta(b^{2n})$ for some base $b$, has a greatest common divisor on the order of $\Theta(b^n)$ for consecutive terms, and where any three consecutive terms are independent, should have a $O(n^{4/3})$ runtime when used as a gap sequence in Shellsort. For example, Sedgewick himself used (among others) the sequence $$k_1 = 1, \quad k_i = (2^i - 3)(2^{i+1} - 3),$$ but it is pretty neat that the product of consecutive Fibonacci numbers can be used here as well. In fact, we can generalize our Fibonacci sequence to a simple way of constructing sequences that match the above set of criteria. Simply start with a sequence $A_n$ where any three consecutive elements of $A$ are pairwise coprime and grow like $\Theta(b^n)$, then define the gap sequence as $$k_i = A_i \cdot A_{i+1}.$$ Real-world performance So… is Fibonacci sort any good in practice? As far as shellsorts go, it is middle-of-the-pack. I also found a similarly structured sequence (with a similar $O(n^{4/3})$ proof) which performs quite a bit better: $$A_1 = 1, \quad A_2 = 2,\quad A_n = 2A_{n-2} + 1$$ $$k_i = A_i \cdot A_{i+1}$$ The advantage of this gap sequence, like the Fibonacci one, is that it is a simple reversible recurrent formula meaning the implementation does not need a lookup table, and can compute the gaps on-the-fly. This is useful, because in my opinion Shellsort only truly has one niche where it shines: tiny code size. I think it’s a solid choice for places where every byte of code matters, while still giving a decent algorithm. I ran a simple benchmark for the number of comparisons performed by the following gap sequences: AlgorithmFormulaComplexity hibbard63$k_i = 2^{i} - 1$$O(n^{3/2})$ pratt71$k = \{2^p3^q\mid p, q \in \mathbb{N}_0\}$$O(n\,(\log n)^2)$ sedge86a$k_1 = 1, k_i = 4^{i} + 3\cdot 2^{i-1}+ 1$$O(n^{4/3})$ sedge86b$k_1 = 1, k_i = (2^i - 3)(2^{i+1} - 3)$$O(n^{4/3})$ sedge86i$k = \{4^j - 3\cdot 2^j + 1\} \cup {\{9\cdot 4^j - 9\cdot 2^j + 1\}}$Unknown lee21$k_i = \lceil\frac{\gamma^i - 1}{\gamma - 1}\rceil, \gamma = 2.243609061420001\dots$Unknown fib25$k_i = F_i \cdot F_{i+1}$$O(n^{4/3})$ orlp25$k_i = A_i \cdot A_{i+1}$$O(n^{4/3})$ I also added heapsort as a $O(n \log n)$ comparison datapoint as it too is a small in-place unstable sort. The results are as follows, with a log scale for size $n$ on the X-axis and the number of comparisons divided by $n \log_2 n$ on the Y-axis: As you can see fib25 doesn’t perform great nor terrible. My other sequence orlp25 performs almost as good as the best experimentally derived sequences, but still provides a worst-case guarantee. Knuth Reward Check The astute observer may have seen sedge86i marked as having an unknown complexity in bold. However, the Wikipedia page for Shellsort at the time of writing lists it as $O(n^{4/3})$. While writing this article I did quite some reading of Sedgewick’s 1986 paper A New Upper Bound for Shellsort, in fact I originally found the paper through a reference on the Wikipedia page for that sequence. While reading the paper I noticed that there are two theorems he proves both of which lead to $O(n^{4/3})$ behavior, Theorem 5 and Theorem 6. Theorem 5 requires as a key property that $k_i, k_{i+1}$, and $k_{i+2}$ are pairwise coprime for all $i$. Theorem 6 requires as a key property that the greatest common divisor of consecutive terms $k_i, k_{i+1}$ always is on the order of $\sqrt{k_i}$ (that’s the theorem I used above). But sedge86i which consists of the merger of two sequences satisfies neither property, nor do the individual sequences before merging! \begin{align*} \{4^j - 3\cdot 2^j + 1\} &= 5, 41, 209, 929, 3905, 16001, 64769, 260609, \dots\\ \{9\cdot 4^j - 9\cdot 2^j + 1\} &= 1, 19, 109, 505, 2161, 8929, 36289, 146305, 587521, \dots \end{align*} As counterexamples in the pre-merged sequences we find $\gcd(209, 3905) = \gcd(36289, 587521) = 11$, and $\gcd(16764929, 37730305) = 29$ in the merged sequence. In fact, if you read the paper, Sedgewick only describes the merged sequence in his conclusion section, as such: The particular sequence used here is a merge of the sequences […]. These occasionally have triples that are not relatively prime, but the combination does better on random inputs than the sequence of Theorem 5 because it has more smaller increments. Sedgewick never claims this particular merged sequence leads to a worst-case of $O(n^{4/3})$. So why is it listed as such on the Wikipedia page? After more digging I realized that the sequence is also found in Knuth’s TAOCP, in Volume 3, Sorting and Searching, Chapter 5.2.1. At the time Knuth wrote: The final examples in Table 6 come from another sequence devised by Sedgewick, based on slightly different heuristics. When these increments $(h_0, h_1, h_2, \dots) = 1, 5, 19, 41, 109, 209, \dots$ are used, Sedgewick proved that the worst-case running time is $O(N^{4/3})$. Except Sedgewick did no such thing (nor did he claim to). I guess whoever edited the Wikipedia article read Knuth’s book and assumed that statement was correct. Contacting Knuth I emailed Knuth with my findings, and a good few months later I received a physical letter containing: my email printed out with a few pencil scribbles, a reward check for finding an error. The errata for Volume 3 now read (emphasis mine): These increments […] combine two sequences that resemble increments for which Sedgewick proved the worst-case time bound $O(N^{4/3})$. I could not resist including a brief synopsis of this blog post with my email, explaining the Fibonacci gap sequence. Knuth wrote back (still scribbled in pencil on my email), “You should experiment with this, and if it performs well or average please publish that fact ASAP!”

It seems that in 2025 a lot of people fall into one of two camps when it comes to AI: skeptic or fanatic. The skeptic thinks AI sucks, that it’s overhyped, it only ever parrots nonsense and it will all blow over soon. The fanatic thinks general human-level intelligence is just around the corner, and that AI will solve almost all our problems. I hope my title is sufficiently ambiguous to attract both camps. The fanatic will be outraged, being ready to jump into the fray to point out why AI isn’t or won’t stay bad. The skeptic will feel validated, and will be eager to read more reasons as to why AI sucks. I’m neither a skeptic nor a fanatic. I see AI more neutrally, as a tool, and from that viewpoint I make the following two observations: AI is bad. It is often incorrect, expensive, racist, trained on data without knowledge or consent, environmentally unfriendly, disruptive to society, etc. AI is useful. Despite the above shortcomings there are tasks for which AI is cheap and effective. I’m no seer, perhaps AI will improve, become more accurate, less biased, cheaper, trained on open access data, cost less electricity, etc. Or perhaps we have plateaued in performance, and there is no political or economic goodwill to address any of the other issues, nor will there be. However, even if AI does not improve in any of the above metrics, it will still be useful, and I hope to show you in this article why. Hence my point: bad AI is here to stay. If you agree with me on this, I hope you’ll also agree with me that we have to stop pretending AI is useless and start taking it and its problems seriously. A formula for query cost Suppose I am a human with some kind of question that can be answered. I know AI could potentially help me with this question, but I wonder if it’s worth it or if I should not use it at all. To help with this we can quantify the risk associated with any potential method of answering the question: $$\mathrm{Risk}_\mathrm{AI} = \mathrm{Cost(query)} + (1 - P(\mathrm{success})) \cdot \mathrm{Cost(bad)}$$ That is, the risk of using any particular method is the cost associated with the method plus the cost of the consequences of a bad answer multiplied by the probability of failure. Here ‘Cost’ is a highly multidimensional object, which can consist of but is not limited to: time, money, environmental impact, ethical concerns, etc. In a lot of cases however we don’t have to blindly trust the answer, and we can verify it. In these cases the consequence of a bad answer is that you’re left in the exact same scenario before trying, except knowing that the AI is of no use. In some scenarios when the AI is non-deterministic it might be worth it to try again as well, but let’s assume for now that you’d have to switch method. In this case the risk is: $$\mathrm{Risk}_\mathrm{AI} = \mathrm{Cost(query)} + \mathrm{Cost(verify)} + (1 - P(\mathrm{success})) \cdot {\mathrm{Risk}}_\mathrm{Other}$$ The cost of a query is usually fairly fixed and known, and although verification cost can vary drastically from task to task, I’d argue that in most cases the cost of verification is also fairly predictable and known. This makes the risk formula applicable in a lot of scenarios, if you have a good idea of the chance of success. The latter, however, can usually only be established empirically, so for one-shot queries without having done any similar queries in the past it can be hard to evaluate whether trying AI is a good idea before doing so. There is one more expansion to the formula I’d like to make before we can look at some examples, and that is to the definition of a successful answer: $$P(\mathrm{success}) = P(\mathrm{correct} \cap \mathrm{relevant})$$ I define a successful answer as one that is both correct and relevant. For example “1 + 1 = 2” might be a correct answer, but irrelevant if we asked about anything else. Relevance is always subjective, but often the correctness of an answer is as well - I’m not assuming here that all questions are about objective facts. Cheap and effective AI queries Because AIs are fallible, usually the biggest cost is in fact the time needed for a human to verify the answer as correct and relevant (or the cost of consequences if left unverified). However, I’ve noticed a real asymmetry between these two properties when it comes to AIs: AIs often give incorrect answers. Worse, they will do so confidently, forcing you to waste time checking their answer instead of them simply stating that they don’t know for sure. AIs almost never give irrelevant answers. If I ask about cheese, the probability a modern AI starts talking about cars is very low. With this in mind I identify five general categories of query for which even bad AI is useful, either by massively reducing or eliminating this verification cost or by leaning on the strong relevance of AI answers: Inspiration, where $\operatorname{Cost}(\mathrm{bad}) \approx 0$, Creative, where $P(\mathrm{correct}) \approx 1$, Planning, where $P(\mathrm{correct}) = P(\mathrm{relevant}) = 1$, Retrieval, where $P(\mathrm{correct}) \approx P(\mathrm{relevant})$, and Objective, where $P(\mathrm{relevant}) = 1$ and correctness verification cost is low. Let’s go over them one by one and look at some examples. Inspiration ($\operatorname{Cost}(\mathrm{bad}) \approx 0$) In this category are the queries where the consequences of a wrong answer are (near) zero. Informally speaking, “it can’t hurt to try”. In my experience these kinds of queries tend to be the ones where you are looking for something but don’t know exactly what; you’ll know it when you see it. For example: “I have leeks, eggs and minced meat in the fridge, as well as a stocked pantry with non-perishable staples. Can you suggest me some dishes I can make with this for a dinner?” “What kind of fun activities can I do with a budget of $100 in New York?” “Suggest some names for a Python function that finds the smallest non-negative number in a list.” “The user wrote this partial paragraph on their phone, suggest three words that are most likely to follow for a quick typing experience.” “Give me 20 synonyms of or similar words to ‘good’.” I think the last query highlights where AI shines or falls for this kind of query. The more localized and personalized your question is, the better the AI will do compared to an alternative. For simple synonyms you can usually just look up the word on a dedicated synonym site, as millions of other people have also wondered the same thing. But the exact contents of your fridge or your exact Python function you’re writing are rather unique to you. Creative ($P(\mathrm{correct}) \approx 1$) In this category are the queries where there are no (almost) no wrong answers. The only thing that really matters is the relevance of the answer, and as I mentioned before, I think AIs are pretty good at being relevant. Examples of queries like these are: “Draw me an image of a polar bear using a computer.” “Write and perform for me a rock ballad about gnomes on tiny bicycles.” “Rephrase the following sentence to be more formal.” “Write a poem to accompany my Sinterklaas gift.” This category does have a controversial aspect to it: it is ‘soulless’, inhuman. Usually if there are no wrong answers we expect the creator to use this opportunity to express their inner thoughts, ideas, experiences and emotions to evoke them in others. If an AI generates art it is not viewed as genuine, even if it evokes the same emotions to those ignorant of the art’s source, because the human to human connection is lost. Current AI models have no inner thoughts, ideas, experiences or emotions, at least not in a way I recognize them. I think it’s fine to use AI art in places where it would otherwise be meaningless (e.g. your corporate presentation slides), fine for humans to use AI-assisted art tools to express themselves, but ultimately defeating the point of art if used as a direct substitute. In the Netherlands we celebrate Sinterklaas which is, roughly speaking, Santa Claus (except we also have Santa Claus, so our children double-dip during the gifting season). Traditionally, gifts from Sinterklaas come accompanied by poems describing the gift and the receiver in a humorous way. It is quite common nowadays, albeit viewed as lazy, to generate such a poem using AI. What’s interesting is that this practice long predates modern LLMs–the poems have such a fixed structure that poem generators have existed a long time. The earliest reference I can find is the 1984 MS-DOS program “Sniklaas”. So people being lazy in supposedly heartfelt art is nothing new. Planning ($P(\mathrm{correct}) = P(\mathrm{relevant}) = 1$) This is a more restrictive form of creativity, where irrelevant answers are absolutely impossible. This often requires some modification of the AI output generation method, where you restrict the output to the valid subdomain (for example yes / no, or binary numbers, etc). However, this is often trivial if you actually have access to the raw model by e.g. masking out invalid outputs, or you are working with a model which outputs the answer directly rather than in natural language or a stream of tokens. One might think in such a restrictive scenario there would be no useful queries, but this isn’t true. The quality of the answer with respect to some (complex) metric might still vary, and AIs might be far better than traditional methods at navigating such domains. For example: “Here is the schema of my database, a SQL query, a small sample of the data and 100 possible query plans. Which query plan seems most likely to execute the fastest? Take into account likely assumptions based on column names and these small data samples.” “What follows is a piece of code. Reformat the code, placing whitespace to maximize readability, while maintaining the exact same syntax tree as per this EBNF grammar.” “Re-order this set of if-else conditions in my code based on your intuition to minimize the expected number of conditions that need to be checked.” “Simplify this math expression using the following set of rewrite rules.” Retrieval ($P(\mathrm{correct}) \approx P(\mathrm{relevant})$) I define retrieval queries as those where the correctness of the answer depends (almost) entirely on its relevance. I’m including classification tasks in this category as well, as one can view it as retrieval of the class from the set of classes (or for binary classification, retrieval of positive samples from a larger set). Then, as long as the cost of verification is low (e.g. a quick glance at a result by a human to see if it interests them), or the consequences of not verifying an irrelevant answer are minimal, AIs can be excellent at this. For example: “Here are 1000 reviews of a restaurant, which ones are overall positive? Which ones mention unsanitary conditions?” “Find me pictures of my dog in my photo collection.” “What are good data structures for maintaining a list of events with dates and quickly counting the number of events in a specified period of time?” “I like Minecraft, can you suggest me some similar games?” “Summarise this 200 page government proposal.” “Which classical orchestral piece starts like ‘da da da daaaaa’”? Objective ($P(\mathrm{relevant}) = 1$, low verification cost) If a problem has an objective answer which can be verified, the relevance of the answer doesn’t really matter or arguably even make sense as a concept. Thus in these cases I’ll define $P(\mathrm{relevant}) = 1$ and leave the cost of verification entirely to correctness. AIs are often incorrect, but not always, so if the primary cost is verification and verification can be done very cheaply or entirely automatically without error, AIs can still be useful despite their fallibility. “What is the mathematical property where a series of numbers can only go up called?” “Identify the car model in this photo.” “I have a list of all Unicode glyphs which are commonly confused with other letters. Can you write an efficient function returning a boolean value that returns true for values in the list but false for all other code points?” “I formalized this mathematical conjecture in Lean. Can you help me write a proof for it?” In a way this category is reminiscent of the $P = NP$ problem. If you have an efficient verification algorithm, is finding solutions still hard in general? The answer seems to be yes, yet the proof eludes us. However, this is only true in general. For specific problems it might very well be possible to use AI to generate provably correct solutions with high probability, even though the the search space is far too large or too complicated for a traditional algorithm. Conclusion Out of the five identified categories, I consider inspiration and retrieval queries to be the strongest use-case for AI where often there is no alternative at all, besides an expensive and slow human that would rather be doing something else. Relevance is highly subjective, complex and fuzzy, which AI handles much better than traditional algorithms. Planning and objective queries are more niche, but absolutely will see use-cases for AI that are hard to replace. Creative queries are both something I think AI is really good at, while simultaneously being the most dangerous and useless category. Art, creativity and human-to-human connections are in my opinion some of the most fundamental aspects of human society, and I think it is incredibly dangerous to mess with them. So dangerous in fact I consider many such queries useless. I wanted the above examples all to be useful queries, so I did not list the following four examples in the “Creative” section despite them belonging there: “Write me ten million personalized spam emails including these links based on the following template.” “Emulate being the perfect girlfriend for me–never disagree with me or challenge my world views like real women do.” “Here is a feed of Reddit threads discussing the upcoming election. Post a comment in each thread, making up a personal anecdote how you are affected by immigrants in a negative way.” “A customer sent in this complaint. Try to help them with any questions they have but if your help is insufficient explain that you are sorry but can not help them any further. Do not reveal you are an AI.” Why do I consider these queries useless, despite them being potentially very profitable or effective? Because their cost function includes such a large detriment to society that only those who ignore its cost to society would ever use them. However, since the cost is “to society” and not to any particular individual, the only way to address this problem is with legislation, as otherwise bad actors are free to harm society for (temporary) personal gain. I wrote this article because I noticed that there are a lot of otherwise intelligent people out there who still believe (or hope) that all AI is useless garbage and that it and its problems will go away by itself. They will not. If you know someone that still believes so, please share this article with them. AI is bad, yes, but bad AI is still useful. Therefore, bad AI is here to stay, and we must deal with it.

Hash functions are incredibly neat mathematical objects. They can map arbitrary data to a small fixed-size output domain such that the mapping is deterministic, yet appears to be random. This “deterministic randomness” is incredibly useful for a variety of purposes, such as hash tables, checksums, monte carlo algorithms, communication-less distributed algorithms, etc, the list goes on. In this article we will take a look at the dark side of hash functions: when things go wrong. Luckily this essentially never happens due to unlucky inputs in the wild (for good hash functions, at least). However, people exist, and some of them may be malicious. Thus we must look towards computer security for answers. I will quickly explain some of the basics of hash function security and then show how easy it is to break this security for some commonly used non-cryptographic hash functions. As a teaser, this article explains how you can generate strings such as these, thousands per second: cityhash64("orlp-cityhash64-D-:K5yx*zkgaaaaa") == 1337 murmurhash2("orlp-murmurhash64-bkiaaa&JInaNcZ") == 1337 murmurhash3("orlp-murmurhash3_x86_32-haaaPa*+") == 1337 farmhash64("orlp-farmhash64-/v^CqdPvziuheaaa") == 1337 I also show how you can create some really funky pairs of strings that can be concatenated arbitrarily such that when concatenating $k$ strings together any of the $2^k$ combinations all have the same hash output, regardless of the seed used for the hash function: a = "xx0rlpx!xxsXъВ" b = "xxsXъВxx0rlpx!" murmurhash2(a + a, seed) == murmurhash2(a + b, seed) murmurhash2(a + a, seed) == murmurhash2(b + a, seed) murmurhash2(a + a, seed) == murmurhash2(b + b, seed) a = "!&orlpՓ" b = "yǏglp$X" murmurhash3(a + a, seed) == murmurhash3(a + b, seed) murmurhash3(a + a, seed) == murmurhash3(b + a, seed) murmurhash3(a + a, seed) == murmurhash3(b + b, seed) Hash function security basics Hash functions play a critical role in computer security. Hash functions are used not only to verify messages over secure channels, they are also used to identify trusted updates as well as known viruses. Virtually every signature scheme ever used starts with a hash function. If a hash function does not behave randomly, we can break the above security constructs. Cryptographic hash functions thus take the randomness aspect very seriously. The ideal hash function would choose an output completely at random for each input, remembering that choice for future calls. This is called a random oracle. The problem is that a random oracle requires a true random number generator, and more problematically, a globally accessible infinite memory bank. So we approximate it using deterministic hash functions instead. These compute their output by essentially shuffling their input really, really well, in such a way that it is not feasible to reverse. To help quantify whether a specific function does a good job of approximating a random oracle, cryptographers came up with a variety of properties that a random oracle would have. The three most important and well-known properties a secure cryptographic hash function should satisfy are: Pre-image resistance. For some constant $c$ it should be hard to find some input $m$ such that $h(m) = c$. Second pre-image resistance. For some input $m_1$ it should be hard to find another input $m_2$ such that $h(m_1) = h(m_2)$. Collision resistance. It should be hard to find inputs $m_1, m_2$ such that $h(m_1) = h(m_2)$. Note that pre-image resistance implies second pre-image resistance which in turn implies collision resistance. Conversely, a pre-image attack breaks all three properties. We generally consider one of these properties broken if there exists a method that produces a collision or pre-image faster than simply trying random inputs (also known as a brute force attack). However, there are definitely gradations in breakage, as some methods are only several orders of magnitude faster than brute force. That may sound like a lot, but a method taking $2^{110}$ steps instead of $2^{128}$ are still both equally out of reach for today’s computers. MD5 used to be a common hash function, and SHA-1 is still in common use today. While both were considered cryptographically secure at one point, generating MD5 collisions now takes less than a second on a modern PC. In 2017 a collaboration of researchers from CWI and Google and announced the first SHA-1 collision. However, as far as I’m aware, neither MD5 nor SHA-1 have practical (second) pre-image attacks, only theoretical ones. Non-cryptographic hash functions Cryptographically secure hash functions tend to have a small problem: they’re slow. Modern hash functions such as BLAKE3 resolve this somewhat by heavily vectorizing the hash using SIMD instructions, as well as parallelizing over multiple threads, but even then they require large input sizes before reaching those speeds. One particular use-case for hash functions is deriving a secret key from a password: a key derivation function. Unlike regular hash functions, being slow is actually a safety feature here to protect against brute forcing passwords. Modern ones such as Argon2 also intentionally use a lot of memory for protection against specialized hardware such as ASICs or FPGAs. A lot of problems don’t necessarily require secure hash functions, and people would much prefer a faster hash speed. Especially when we are computing many small hashes, such as in a hash table. Let’s take a look what common hash table implementations actually use as their hash for strings: C++: there are multiple standard library implementations, but 64-bit clang 13.0.0 on Apple M1 ships CityHash64. Currently libstdc++ ships MurmurHash64A, a variant of Murmur2 for 64-bit platforms. Java: OpenJDK uses an incredibly simple hash algorithm, which essentially just computes h = 31 * h + c for each character c. PHP: the Zend engine uses essentially the same algorithm as Java, just using unsigned integers and 33 as its multiplier. Nim: it used to use MurmurHash3_x86_32. While writing this article they appeared to have switched to use farmhash by default. Zig: it uses wyhash by default, with 0 as seed. Javascript: in V8 they use a custom weak string hash, with a randomly initialized seed. There were some that used stronger hashes by default as well: Go uses an AES-based hash if hardware acceleration is available on x86-64. Even though its construction is custom and likely not full-strength cryptographically secure, breaking it is too much effort and quite possibly beyond my capabilities. If not available, it uses an algorithm inspired by wyhash. Python and Rust use SipHash by default, which is a cryptographically secure pseudorandom function. This is effectively a hash function where you’re allowed to use a secret key during hashing, unlike a hash like SHA-2 where everyone knows all information involved. This latter concept is actually really important, at least for protecting against HashDoS in hash tables. Even if a hash function is perfectly secure over its complete output, hash tables further reduce the output to only a couple bits to find the data it is looking for. For a static hash function without any randomness it’s possible to produce large lists of hashes that collide post-reduction, just by brute force. But for non-cryptographic hashes as we’ll see here we often don’t need brute force and can generate collisions at high speed for the full output, if not randomized by a random seed. Interlude: inverse operations Before we get to breaking some of the above hash functions, I must explain a basic technique I will use a lot: the inverting of operations. We are first exposed to this in primary school, where we might get faced by a question such as “$2 + x = 10$”. There we learn subtraction is the inverse of addition, such that we may find $x$ by computing $10 - 2 = 8$. Most operations on the integer registers in computers are also invertible, despite the integers being reduced modulo $2^{w}$ in the case of overflow. Let us study some: Addition can be inverted using subtraction. That is, x += y can be inverted using x -= y. Seems obvious enough. Multiplication by a constant $c$ is not inverted by division. This would not work in the case of overflow. Instead, we calculate the modular multiplicative inverse of $c$. This is an integer $c^{-1}$ such that $c \cdot c^{-1} \equiv 1 \pmod {m}$. Then we invert multiplication by $c$ simply by multiplying by $c^{-1}$. This constant exists if and only if $c$ is coprime with our modulus $m$, which for us means that $c$ must be odd as $m = 2^n$. For example, multiplication by $2$ is not invertible, which is easy to see as such, as it is equivalent to a bit shift to the left by one position, losing the most significant bit forever. Without delving into the details, here is a snippet of Python code that computes the modular multiplicative inverse of an integer using the extended Euclidean algorithm by calculating $x, y$ such that $$cx + my = \gcd(c, m).$$ Then, because $c$ is coprime we find $\gcd(c, m) = 1$, which means that $$cx + 0 \equiv 1 \pmod m,$$ and thus $x = c^{-1}$. def egcd(a, b): if a == 0: return (b, 0, 1) g, y, x = egcd(b % a, a) return (g, x - (b // a) * y, y) def modinv(c, m): g, x, y = egcd(c, m) assert g == 1, "c, m must be coprime" return x % m Using this we can invert modular multiplication: >>> modinv(17, 2**32) 4042322161 >>> 42 * 17 * 4042322161 % 2**32 42 Magic! XOR can be inverted using… XOR. It is its own inverse. So x ^= y can be inverted using x ^= y. Bit shifts can not be inverted, but two common operations in hash functions that use bit shifts can be. The first is bit rotation by a constant. This is best explained visually, for example a bit rotation to the left by 3 places on a 8-bit word, where each bit is shown as a letter: abcdefghi defghiabc The formula for a right-rotation of k places is (x >> k) | (x << (w - k)), where w is the width of the integer type. Its inverse is a left-rotation, which simply swaps the direction of both shifts. Alternatively, the inverse of a right-rotation of k places is another right-rotation of w-k places. Another common operation in hash functions is the “xorshift”. It is an operation of one of the following forms, with $k > 0$: x ^= x << k // Left xorshift. x ^= x >> k // Right xorshift. How to invert it is entirely analogous between the two, so I will focus on the left xorshift. An important observation is that the least significant $k$ bits are left entirely untouched by the xorshift. Thus by repeating the operation, we recover the least significant $2k$ bits, as the XOR will invert itself for the next $k$ bits. Let’s take a look at the resulting value to see how we should proceed: v0 = (x << k) ^ x // Apply first step of inverse v1 = v0 ^ (v0 << k). v1 = (x << 2*k) ^ (x << k) ^ (x << k) ^ x // Simplify using self-inverse (x << k) ^ (x << k) = 0. v1 = (x << 2*k) ^ x From this we can conclude the following identity: $$\operatorname{xorshift}(\operatorname{xorshift}(x, k), k) = \operatorname{xorshift}(x, 2k)$$ Now we only need one more observation to complete our algorithm: a xorshift of $k \geq w$ where $w$ is the width of our integer is a no-op. Thus we repeatedly apply our doubling identity until we reach large enough $q$ such that $\operatorname{xorshift}(x, 2^q \cdot k) = x$. For example, to invert a left xorshift by 13 for 64-bit integers we apply the following sequence: x ^= x << 13 // Left xorshift by 13. x ^= x << 13 // Inverse step 1. x ^= x << 26 // Inverse step 2. x ^= x << 52 // Inverse step 3. // x ^= x << 104 // Next step would be a no-op. Armed with this knowledge, we can now attack. Breaking CityHash64 Let us take a look at (part of) the source code of CityHash64 from libcxx that’s used for hashing strings on 64-bit platforms: C++ standard library code goes through a process known as 'uglification', which prepends underscores to all identifiers. This is because those identifiers are reserved by the standard to only be used in the standard library, and thus won't be replaced by macros from standards-compliant code. For your sanity's sake I removed them here. static const uint64_t mul = 0x9ddfea08eb382d69ULL; static const uint64_t k0 = 0xc3a5c85c97cb3127ULL; static const uint64_t k1 = 0xb492b66fbe98f273ULL; static const uint64_t k2 = 0x9ae16a3b2f90404fULL; static const uint64_t k3 = 0xc949d7c7509e6557ULL; template<class T> T loadword(const void* p) { T r; std::memcpy(&r, p, sizeof(r)); return r; } uint64_t rotate(uint64_t val, int shift) { if (shift == 0) return val; return (val >> shift) | (val << (64 - shift)); } uint64_t hash_len_16(uint64_t u, uint64_t v) { uint64_t x = u ^ v; x *= mul; x ^= x >> 47; uint64_t y = v ^ x; y *= mul; y ^= y >> 47; y *= mul; return y; } uint64_t hash_len_17_to_32(const char *s, uint64_t len) { const uint64_t a = loadword<uint64_t>(s) * k1; const uint64_t b = loadword<uint64_t>(s + 8); const uint64_t c = loadword<uint64_t>(s + len - 8) * k2; const uint64_t d = loadword<uint64_t>(s + len - 16) * k0; return hash_len_16( rotate(a - b, 43) + rotate(c, 30) + d, a + rotate(b ^ k3, 20) - c + len ); } To break this, let’s assume we’ll always give length 32 inputs. Then the implementation will always call hash_len_17_to_32, and we have full control over variables a, b, c and d by changing our input. Note that d is only used once, in the final expression. This makes it a prime target for attacking the hash. We will choose a, b and c arbitrarily, and then solve for d to compute a desired hash outcome. Using the above modinv function we first compute the necessary modular multiplicative inverses of mul and k0: >>> 0x9ddfea08eb382d69 * 0xdc56e6f5090b32d9 % 2**64 1 >>> 0xc3a5c85c97cb3127 * 0x81bc9c5aa9c72e97 % 2**64 1 We also note that in this case the xorshift is easy to invert, as x ^= x >> 47 is simply its own inverse. Having all the components ready, we can invert the function step by step. We first load a, b and c like in the hash function, and compute uint64_t v = a + rotate(b ^ k3, 20) - c + len; which is the second parameter to hash_len_16. Then, starting from our desired return value of hash_len_16(u, v) we work backwards step by step, inverting each operation to find the function argument u that would result in our target hash. Then once we have found such the unique u we compute our required input d. Putting it all together: static const uint64_t mul_inv = 0xdc56e6f5090b32d9ULL; static const uint64_t k0_inv = 0x81bc9c5aa9c72e97ULL; void cityhash64_preimage32(uint64_t hash, char *s) { const uint64_t len = 32; const uint64_t a = loadword<uint64_t>(s) * k1; const uint64_t b = loadword<uint64_t>(s + 8); const uint64_t c = loadword<uint64_t>(s + len - 8) * k2; uint64_t v = a + rotate(b ^ k3, 20) - c + len; // Invert hash_len_16(u, v). Original operation inverted // at each step is shown on the right, note that it is in // the inverse order of hash_len_16. uint64_t y = hash; // return y; y *= mul_inv; // y *= mul; y ^= y >> 47; // y ^= y >> 47; y *= mul_inv; // y *= mul; uint64_t x = y ^ v; // uint64_t y = v ^ x; x ^= x >> 47; // x ^= x >> 47; x *= mul_inv; // x *= mul; uint64_t u = x ^ v; // uint64_t x = u ^ v; // Find loadword<uint64_t>(s + len - 16). uint64_t d = u - rotate(a - b, 43) - rotate(c, 30); d *= k0_inv; std::memcpy(s + len - 16, &d, sizeof(d)); } The chance that a random uint64_t forms 8 printable ASCII bytes is $\left(94/256\right)^8 \approx 0.033%$. Not great, but cityhash64_preimage32 is so fast that having to repeat it on average ~3000 times to get a purely ASCII result isn’t so bad. For example, the following 10 strings all hash to 1337 using CityHash64, generated using this code: I’ve noticed there’s variants of CityHash64 with subtle differences in the wild. I chose to attack the variant shipped with libc++, so it should work for std::hash there, for example. I also assume a little-endian machine throughout this article, your mileage may vary on a big-endian machine depending on the hash implementation. orlp-cityhash64-D-:K5yx*zkgaaaaa orlp-cityhash64-TXb7;1j&btkaaaaa orlp-cityhash64-+/LM$0 ;msnaaaaa orlp-cityhash64-u'f&>I'~mtnaaaaa orlp-cityhash64-pEEv.LyGcnpaaaaa orlp-cityhash64-v~~bm@,Vahtaaaaa orlp-cityhash64-RxHr_&~{miuaaaaa orlp-cityhash64-is_$34#>uavaaaaa orlp-cityhash64-$*~l\{S!zoyaaaaa orlp-cityhash64-W@^5|3^:gtcbaaaa Breaking MurmurHash2 We can’t let libstdc++ get away after targetting libc++, can we? The default string hash calls an implementation of MurmurHash2 with seed 0xc70f6907. The hash—simplified to only handle strings whose lengths are multiples of 8—is as follows: uint64_t murmurhash64a(const char* s, size_t len, uint64_t seed) { const uint64_t mul = 0xc6a4a7935bd1e995ULL; uint64_t hash = seed ^ (len * mul); for (const char* p = s; p != s + len; p += 8) { uint64_t data = loadword<uint64_t>(p); data *= mul; data ^= data >> 47; data *= mul; hash ^= data; hash *= mul; } hash ^= hash >> 47; hash *= mul; hash ^= hash >> 47; return hash; } We can take a similar approach here as before. We note that the modular multiplicative inverse of 0xc6a4a7935bd1e995 mod $2^{64}$ is 0x5f7a0ea7e59b19bd. As an example, we can choose the first 24 bytes arbitrarily, and solve for the last 8 bytes: void murmurhash64a_preimage32(uint64_t hash, char* s, uint64_t seed) { const uint64_t mul = 0xc6a4a7935bd1e995ULL; const uint64_t mulinv = 0x5f7a0ea7e59b19bdULL; // Compute the hash state for the first 24 bytes as normal. uint64_t state = seed ^ (32 * mul); for (const char* p = s; p != s + 24; p += 8) { uint64_t data = loadword<uint64_t>(p); data *= mul; data ^= data >> 47; data *= mul; state ^= data; state *= mul; } // Invert target hash transformation. // return hash; hash ^= hash >> 47; // hash ^= hash >> 47; hash *= mulinv; // hash *= mul; hash ^= hash >> 47; // hash ^= hash >> 47; // Invert last iteration for last 8 bytes. hash *= mulinv; // hash *= mul; uint64_t data = state ^ hash; // hash = hash ^ data; data *= mulinv; // data *= mul; data ^= data >> 47; // data ^= data >> 47; data *= mulinv; // data *= mul; std::memcpy(s + 24, &data, 8); // data = loadword<uint64_t>(s); } The following 10 strings all hash to 1337 using MurmurHash64A with the default seed 0xc70f6907, generated using this code: orlp-murmurhash64-bhbaaat;SXtgVa orlp-murmurhash64-bkiaaa&JInaNcZ orlp-murmurhash64-ewmaaa(%J+jw>j orlp-murmurhash64-vxpaaag"93\Yj5 orlp-murmurhash64-ehuaaafa`Wp`/| orlp-murmurhash64-yizaaa1x.zQF6r orlp-murmurhash64-lpzaaaZphp&c F orlp-murmurhash64-wsjbaa771rz{z< orlp-murmurhash64-rnkbaazy4X]p>B orlp-murmurhash64-aqnbaaZ~OzP_Tp Universal collision attack on MurmurHash64A In fact, MurmurHash64A is so weak that Jean-Philippe Aumasson, Daniel J. Bernstein and Martin Boßlet published an attack that creates sets of strings which collide regardless of the random seed used. To be fair to CityHash64… just kidding they found universal collisions against it as well, regardless of seed used. CityHash64 is actually much easier to break in this way, as simply doing the above pre-image attack targetting 0 as hash makes the output purely dependent on the seed, and thus a universal collision. To see how it works, let’s take a look at the core loop of MurmurHash64A: uint64_t data = loadword<uint64_t>(p); data *= mul; // Trivially invertible. data ^= data >> 47; // Trivially invertible. data *= mul; // Trivially invertible. state ^= data; state *= mul; We know we can trivially invert the operations done on data regardless of what the current state is, so we might as well have had the following body: state ^= data; state *= mul; Now the hash starts looking rather weak indeed. The clever trick they employ is by creating two strings simultaneously, such that they differ precisely in the top bit in each 8-byte word. Why the top bit? >>> 1 << 63 9223372036854775808 >>> (1 << 63) * mul % 2**64 9223372036854775808 Since mul is odd, its least significant bit is set. Multiplying 1 << 63 by it is equivalent to shifting that bit 63 places to the left, which is once again 1 << 63. That is, 1 << 63 is a fixed point for the state *= mul operation. We also note that for the top bit XOR is equivalent to addition, as the overflow from addition is removed mod $2^{64}$. So if we have two input strings, one starting with the 8 bytes data, and the other starting with data ^ (1 << 63) == data + (1 << 63) (after doing the trivial inversions). We then find that the two states, regardless of seed, differ exactly in the top bit after state ^= data. After multiplication we find we have two states x * mul and (x + (1 << 63)) * mul == x * mul + (1 << 63)… which again differ exactly in the top bit! We are now back to state ^= data in our iteration, for the next 8 bytes. We can now use this moment to cancel our top bit difference, by again feeding two 8-byte strings that differ in the top bit (after inverting). In fact, we only have to find one pair of such strings that differ in the top bit, which we can then repeat twice (in either order) to cancel our difference again. When represented as a uint64_t if we choose the first string as x we can derive the second string as x *= mul; // Forward transformation... x ^= x >> 47; // ... x *= mul; // ... x ^= 1 << 63; // Difference in top bit. x *= mulinv; // Backwards transformation... x ^= x >> 47; // ... x *= mulinv; // ... I was unable to find a printable ASCII string that has another printable ASCII string as its partner. But I was able to find the following pair of 8-byte UTF-8 strings that differ in exactly the top bit after the Murmurhash64A input transformation: xx0rlpx! xxsXъВ Combining them as such gives two 16-byte strings that when fed through the hash algorithm manipulate the state in the same way: a collision. xx0rlpx!xxsXъВ xxsXъВxx0rlpx! But it doesn’t stop there. By concatenating these two strings we can create $2^n$ different colliding strings each $16n$ bytes long. With the current libstdc++ implementation the following prints the same number eight times: std::hash<std::u8string> h; std::u8string a = u8"xx0rlpx!xxsXъВ"; std::u8string b = u8"xxsXъВxx0rlpx!"; std::cout << h(a + a + a) << "\n"; std::cout << h(a + a + b) << "\n"; std::cout << h(a + b + a) << "\n"; std::cout << h(a + b + b) << "\n"; std::cout << h(b + a + a) << "\n"; std::cout << h(b + a + b) << "\n"; std::cout << h(b + b + a) << "\n"; std::cout << h(b + b + b) << "\n"; Even if the libstdc++ would randomize the seed used by MurmurHash64a, the strings would still collide. Breaking MurmurHash3 Nim uses used to use MurmurHash3_x86_32, so let’s try to break that. If we once again simplify to strings whose lengths are a multiple of 4 we get the following code: uint32_t rotl32(uint32_t x, int r) { return (x << r) | (x >> (32 - r)); } uint32_t murmurhash3_x86_32(const char* s, int len, uint32_t seed) { const uint32_t c1 = 0xcc9e2d51; const uint32_t c2 = 0x1b873593; const uint32_t c3 = 0x85ebca6b; const uint32_t c4 = 0xc2b2ae35; uint32_t h = seed; for (const char* p = s; p != s + len; p += 4) { uint32_t k = loadword<uint32_t>(p); k *= c1; k = rotl32(k, 15); k *= c2; h ^= k; h = rotl32(h, 13); h = h * 5 + 0xe6546b64; } h ^= len; h ^= h >> 16; h *= c3; h ^= h >> 13; h *= c4; h ^= h >> 16; return h; } I think by now you should be able to get this function to spit out any value you want if you know the seed. The inverse of rotl32(x, r) is rotl32(x, 32-r) and the inverse of h ^= h >> 16 is once again just h ^= h >> 16. Only h ^= h >> 13 is a bit different, it’s the first time we’ve seen that a xorshift’s inverse has more than one step: h ^= h >> 13 h ^= h >> 26 Compute the modular inverses of c1 through c4 as well as 5 mod $2^{32}$, and go to town. If you want to cheat or check your answer, you can check out the code I’ve used to generate the following ten strings that all hash to 1337 when fed to MurmurHash3_x86_32 with seed 0: orlp-murmurhash3_x86_32-haaaPa*+ orlp-murmurhash3_x86_32-saaaUW&< orlp-murmurhash3_x86_32-ubaa/!/" orlp-murmurhash3_x86_32-weaare]] orlp-murmurhash3_x86_32-chaa5@/} orlp-murmurhash3_x86_32-claaM[,5 orlp-murmurhash3_x86_32-fraaIx`N orlp-murmurhash3_x86_32-iwaara&< orlp-murmurhash3_x86_32-zwaa]>zd orlp-murmurhash3_x86_32-zbbaW-5G Nim uses 0 as a fixed seed. You might wonder about the ethics of publishing functions for generating arbitrary amounts of collisions for hash functions actually in use today. I did consider holding back. But HashDoS has been a known attack for almost two decades now, and the universal hash collisions I’ve shown were also published more than a decade ago now as well. At some point you’ve had enough time to, uh, fix your shit. Universal collision attack on MurmurHash3 Suppose that Nim didn’t use 0 as a fixed seed, but chose a randomly generated one. Can we do a similar attack as the one done to MurmurHash2 to still generate universal multicollisions? Yes we can. Let’s take another look at that core loop body: uint32_t k = loadword<uint32_t>(p); k *= c1; // Trivially invertable. k = rotl32(k, 15); // Trivially invertable. k *= c2; // Trivially invertable. h ^= k; h = rotl32(h, 13); h = h * 5 + 0xe6546b64; Once again we can ignore the first three trivially invertable instructions as we can simply choose our input so that we get exactly the k we want. Remember from last time that we want to introduce a difference in exactly the top bit of h, as the multiplication will leave this difference in place. But here there is a bit rotation between the XOR and the multiplication. The solution? Simply place our bit difference such that rotl32(h, 13) shifts it into the top position. Does the addition of 0xe6546b64 mess things up? No. Since only the top bit between the two states will be different, there is a difference of exactly $2^{31}$ between the two states. This difference is maintained by the addition. Since two 32-bit numbers with the same top bit can be at most $2^{31} - 1$ apart, we can conclude that the two states still differ in the top bit after the addition. So we want to find two pairs of 32-bit ints, such that after applying the first three instructions the first pair differs in bit 1 << (31 - 13) == 0x00040000 and the second pair in bit 1 << 31 == 0x80000000. After some brute-force searching I found some cool pairs (again forced to use UTF-8), which when combined give the following collision: a = "!&orlpՓ" b = "yǏglp$X" As before, any concatenation of as and bs of length n collides with all other combinations of length n. Breaking FarmHash64 Nim switched to farmhash since I started writing this post. To break it we can notice that its structure is very similar to CityHash64, so we can use those same techniques again. In fact, the only changes between the two for lengths 17-32 bytes is that a few operators were changed from subtraction/XOR to addition, a rotation operator had its constant tweaked, and some k constants are slightly tweaked in usage. The process of breaking it is so similar that it’s entirely analogous, so we can skip straight to the result. These 10 strings all hash to 1337 with FarmHash64: orlp-farmhash64-?VrJ@L7ytzwheaaa orlp-farmhash64-p3`!SQb}fmxheaaa orlp-farmhash64-pdt'cuI\gvxheaaa orlp-farmhash64-IbY`xAG&ibkieaaa orlp-farmhash64-[_LU!d1hwmkieaaa orlp-farmhash64-QiY!clz]bttieaaa orlp-farmhash64-&?J3rZ_8gsuieaaa orlp-farmhash64-LOBWtm5Szyuieaaa orlp-farmhash64-Mptaa^g^ytvieaaa orlp-farmhash64-B?&l::hxqmfjeaaa Trivial fixed-seed wyhash multicollisions Zig uses wyhash with a fixed seed of zero. While I was unable to do seed-independent attacks against wyhash, using it with a fixed seed makes generating collisions trivial. Wyhash is built upon the folded multiply, which takes two 64-bit inputs, multiplies them to a 128-bit product before XORing together the two halves: uint64_t folded_multiply(uint64_t a, uint64_t b) { __uint128_t full = __uint128_t(a) * __uint128_t(b); return uint64_t(full) ^ uint64_t(full >> 64); } It’s easy to immediately see a critical flaw with this: if one of the two sides is zero, the output will also always be zero. To protect against this, wyhash always uses a folded multiply in the following form: out = folded_multiply(input_a ^ secret_a, input_b ^ secret_b); where secret_a and secret_b are determined by the seed, or outputs of previous iterations which are influenced by the seed. However, when your seed is constant… With a bit of creativity we can use the start of our string to prepare a ‘secret’ value which we can perfectly cancel with another ASCII string later in the input. So, without further ado, every 32-byte string of the form orlp-wyhash-oGf_________tWJbzMJR hashes to the same value with Zig’s default hasher. Zig uses a different set of parameters than the defaults found in the wyhash repository, so for good measure, this pattern provides arbitrary multicollisions for the default parameters found in wyhash when using seed == 0: orlp-wyhash-EUv_________NLXyytkp Conclusion We’ve seen that a lot of the hash functions in common use in hash tables today are very weak, allowing fairly trivial attacks to produce arbitrary amounts of collisions if not randomly initialized. Using a randomly seeded hash table is paramount if you don’t wish to become a victim of a hash flooding attack. We’ve also seen that some hash functions are vulnerable to attack even if randomly seeded. These are completely broken and should not be used if attacks are a concern at all. Luckily I was unable to find such attacks against most hashes, but the possibility of such an attack existing is quite unnerving. With universal hashing it’s possible to construct hash functions for which such an attack is provably impossible, last year I published a hash function called polymur-hash that has this property. Your HTTPS connection to this website also likely uses a universal hash function for authenticity of the transferred data, both Poly1305 and GCM are based on universal hashing for their security proofs. Well, such attacks are provably impossible against non-interactive attackers, everything goes out of the window again when an attacker is allowed to inspect the output hashes and use that to try and guess your secret key. Of course, if your data is not user-controlled, or there is no reasonable security model where your application would face attacks, you can get away with faster and insecure hashes. More to come on the subject of hashing and hash tables and how it can go right or wrong, but for now this article is long enough as-is…

Suppose you have an array of floating-point numbers, and wish to sum them. You might naively think you can simply add them, e.g. in Rust: fn naive_sum(arr: &[f32]) -> f32 { let mut out = 0.0; for x in arr { out += *x; } out } This however can easily result in an arbitrarily large accumulated error. Let’s try it out: naive_sum(&vec![1.0; 1_000_000]) = 1000000.0 naive_sum(&vec![1.0; 10_000_000]) = 10000000.0 naive_sum(&vec![1.0; 100_000_000]) = 16777216.0 naive_sum(&vec![1.0; 1_000_000_000]) = 16777216.0 Uh-oh… What happened? When you compute $a + b$ the result must be rounded to the nearest representable floating-point number, breaking ties towards the number with an even mantissa. The problem is that the next 32-bit floating-point number after 16777216 is 16777218. In this case that means 16777216 + 1 rounds back to 16777216 again. We’re stuck. Luckily, there are better ways to sum an array. Pairwise summation A method that’s a bit more clever is to use pairwise summation. Instead of a completely linear sum with a single accumulator it recursively sums an array by splitting the array in half, summing the halves, and then adding the sums. fn pairwise_sum(arr: &[f32]) -> f32 { if arr.len() == 0 { return 0.0; } if arr.len() == 1 { return arr[0]; } let (first, second) = arr.split_at(arr.len() / 2); pairwise_sum(first) + pairwise_sum(second) } This is more accurate: pairwise_sum(&vec![1.0; 1_000_000]) = 1000000.0 pairwise_sum(&vec![1.0; 10_000_000]) = 10000000.0 pairwise_sum(&vec![1.0; 100_000_000]) = 100000000.0 pairwise_sum(&vec![1.0; 1_000_000_000]) = 1000000000.0 However, this is rather slow. To get a summation routine that goes as fast as possible while still being reasonably accurate we should not recurse down all the way to length-1 arrays, as this gives too much call overhead. We can still use our naive sum for small sizes, and only recurse on large sizes. This does make our worst-case error worse by a constant factor, but in turn makes the pairwise sum almost as fast as a naive sum. By choosing the splitpoint as a multiple of 256 we ensure that the base case in the recursion always has exactly 256 elements except on the very last block. This makes sure we use the most optimal reduction and always correctly predict the loop condition. This small detail ended up improving the throughput by 40% for large arrays! fn block_pairwise_sum(arr: &[f32]) -> f32 { if arr.len() > 256 { let split = (arr.len() / 2).next_multiple_of(256); let (first, second) = arr.split_at(split); block_pairwise_sum(first) + block_pairwise_sum(second) } else { naive_sum(arr) } } Kahan summation The worst-case round-off error of naive summation scales with $O(n \epsilon)$ when summing $n$ elements, where $\epsilon$ is the machine epsilon of your floating-point type (here $2^{-24}$). Pairwise summation improves this to $O((\log n) \epsilon + n\epsilon^2)$. However, Kahan summation improves this further to $O(n\epsilon^2)$, eliminating the $\epsilon$ term entirely, leaving only the $\epsilon^2$ term which is negligible unless you sum a very large amount of numbers. All of these bounds scale with $\sum_i |x_i|$, so the worst-case absolute error bound is still quadratic in terms of $n$ even for Kahan summation. In practice all summation algorithms do significantly better than their worst-case bounds, as in most scenarios the errors do not exclusively round up or down, but cancel each other out on average. pub fn kahan_sum(arr: &[f32]) -> f32 { let mut sum = 0.0; let mut c = 0.0; for x in arr { let y = *x - c; let t = sum + y; c = (t - sum) - y; sum = t; } sum } The Kahan summation works by maintaining the sum in two registers, the actual bulk sum and a small error correcting term $c$. If you were using infinitely precise arithmetic $c$ would always be zero, but with floating-point it might not be. The downside is that each number now takes four operations to add to the sum instead of just one. To mitigate this we can do something similar to what we did with the pairwise summation. We can first accumulate blocks into sums naively before combining the block sums with Kaham summation to reduce overhead at the cost of accuracy: pub fn block_kahan_sum(arr: &[f32]) -> f32 { let mut sum = 0.0; let mut c = 0.0; for chunk in arr.chunks(256) { let x = naive_sum(chunk); let y = x - c; let t = sum + y; c = (t - sum) - y; sum = t; } sum } Exact summation I know of at least two general methods to produce the correctly-rounded sum of a sequence of floating-point numbers. That is, it logically computes the sum with infinite precision before rounding it back to a floating-point value at the end. The first method is based on the 2Sum primitive which is an error-free transform from two numbers $x, y$ to $s, t$ such that $x + y = s + t$, where $t$ is a small error. By applying this repeatedly until the errors vanish you can get a correctly-rounded sum. Keeping track of what to add in what order can be tricky, and the worst-case requires $O(n^2)$ additions to make all the terms vanish. This is what’s implemented in Python’s math.fsum and in the Rust crate fsum which use extra memory to keep the partial sums around. The accurate crate also implements this using in-place mutation in i_fast_sum_in_place. Another method is to keep a large buffer of integers around, one per exponent. Then when adding a floating-point number you decompose it into a an exponent and mantissa, and add the mantissa to the corresponding integer in the buffer. If the integer buf[i] overflows you increment the integer in buf[i + w], where w is the width of your integer. This can actually compute a completely exact sum, without any rounding at all, and is effectively just an overly permissive representation of a fixed-point number optimized for accumulating floats. This latter method is $O(n)$ time, but uses a large but constant amount of memory ($\approx$ 1 KB for f32, $\approx$ 16 KB for f64). An advantage of this method is that it’s also an online algorithm - both adding a number to the sum and getting the current total are amortized $O(1)$. A variant of this method is implemented in the accurate crate as OnlineExactSum crate which uses floats instead of integers for the buffer. Unleashing the compiler Besides accuracy, there is another problem with naive_sum. The Rust compiler is not allowed to reorder floating-point additions, because floating-point addition is not associative. So it cannot autovectorize the naive_sum to use SIMD instructions to compute the sum, nor use instruction-level parallelism. To solve this there are compiler intrinsics in Rust that do float sums while allowing associativity, such as std::intrinsics::fadd_fast. However, these instructions are incredibly dangerous, as they assume that both the input and output are finite numbers (no infinities, no NaNs), or otherwise they are undefined behavior. This functionally makes them unusable, as only in the most restricted scenarios when computing a sum do you know that all inputs are finite numbers, and that their sum cannot overflow. I recently uttered my annoyance with these operators to Ben Kimock, and together we proposed (and he implemented) a new set of operators: std::intrinsics::fadd_algebraic and friends. I proposed we call the operators algebraic, as they allow (in theory) any transformation that is justified by real algebra. For example, substituting ${x - x \to 0}$, ${cx + cy \to c(x + y)}$, or ${x^6 \to (x^2)^3.}$ In general these operators are treated as-if they are done using real numbers, and can map to any set of floating-point instructions that would be equivalent to the original expression, assuming the floating-point instructions would be exact. Note that the real numbers do not contain NaNs or infinities, so these operators assume those do not exist for the validity of transformations, however it is not undefined behavior when you do encounter those values. They also allow fused multiply-add instructions to be generated, as under real arithmetic $\operatorname{fma}(a, b, c) = ab + c.$ Using those new instructions it is trivial to generate an autovectorized sum: #![allow(internal_features)] #![feature(core_intrinsics)] use std::intrinsics::fadd_algebraic; fn naive_sum_autovec(arr: &[f32]) -> f32 { let mut out = 0.0; for x in arr { out = fadd_algebraic(out, *x); } out } If we compile with -C target-cpu=broadwell we see that the compiler automatically generated the following tight loop for us, using 4 accumulators and AVX2 instructions: .LBB0_5: vaddps ymm0, ymm0, ymmword ptr [rdi + 4*r8] vaddps ymm1, ymm1, ymmword ptr [rdi + 4*r8 + 32] vaddps ymm2, ymm2, ymmword ptr [rdi + 4*r8 + 64] vaddps ymm3, ymm3, ymmword ptr [rdi + 4*r8 + 96] add r8, 32 cmp rdx, r8 jne .LBB0_5 This will process 128 bytes of floating-point data (so 32 elements) in 7 instructions. Additionally, all the vaddps instructions are independent of each other as they accumulate to different registers. If we analyze this with uiCA we see that it estimates the above loop to take 4 cycles to complete, processing 32 bytes / cycle. At 4GHz that’s up to 128GB/s! Note that that’s way above what my machine’s RAM bandwidth is, so you will only achieve that speed when summing data that is already in cache. With this in mind we can also easily define block_pairwise_sum_autovec and block_kahan_sum_autovec by replacing their calls to naive_sum with naive_sum_autovec. Accuracy and speed Let’s take a look at how the different summation methods compare. As a relatively arbitrary benchmark, let’s sum 100,000 random floats ranging from -100,000 to +100,000. This is 400 KB worth of data, so it still fits in cache on my AMD Threadripper 2950x. All the code is available on Github. Compiled with RUSTFLAGS=-C target-cpu=native and --release I get the following results: AlgorithmThroughputMean absolute error naive5.5 GB/s71.796 pairwise0.9 GB/s1.5528 kahan1.4 GB/s0.2229 block_pairwise5.8 GB/s3.8597 block_kahan5.9 GB/s4.2184 naive_autovec118.6 GB/s14.538 block_pairwise_autovec71.7 GB/s1.6132 block_kahan_autovec98.0 GB/s1.2306 crate_accurate_buffer1.1 GB/s0.0015 crate_accurate_inplace1.9 GB/s0.0015 crate_fsum1.2 GB/s0.0000 The reason the accurate crate has a non-zero absolute error is because it currently does not implement rounding to nearest correctly, so it can be off by one unit in the last place for the final result. First I’d like to note that there’s more than a 100x performance difference between the fastest and slowest method. For summing an array! Now this might not be entirely fair as the slowest methods are computing something significantly harder, but there’s still a 20x performance difference between a seemingly reasonable naive implementation and the fastest one. We find that in general the _autovec methods that use fadd_algebraic are faster and more accurate than the ones using regular floating-point addition. The reason they’re more accurate as well is the same reason a pairwise sum is more accurate: any reordering of the additions is better as the default long-chain-of-additions is already the worst case for accuracy in a sum. Limiting ourselves to Pareto-optimal choices we get the following four implementations: AlgorithmThroughputMean absolute error naive_autovec118.6 GB/s14.538 block_kahan_autovec98.0 GB/s1.2306 crate_accurate_inplace1.9 GB/s0.0015 crate_fsum1.2 GB/s0.0000 Note that implementation differences can be quite impactful, and there are likely dozens more methods of compensated summing I did not compare here. For most cases I think block_kahan_autovec wins here, having good accuracy (that doesn’t degenerate with larger inputs) at nearly the maximum speed. For most applications the extra accuracy from the correctly-rounded sums is unnecessary, and they are 50-100x slower. By splitting the loop up into an explicit remainder plus a tight loop of 256-element sums we can squeeze out a bit more performance, and avoid a couple floating-point ops for the last chunk: #![allow(internal_features)] #![feature(core_intrinsics)] use std::intrinsics::fadd_algebraic; fn sum_block(arr: &[f32]) -> f32 { arr.iter().fold(0.0, |x, y| fadd_algebraic(x, *y)) } pub fn sum_orlp(arr: &[f32]) -> f32 { let mut chunks = arr.chunks_exact(256); let mut sum = 0.0; let mut c = 0.0; for chunk in &mut chunks { let y = sum_block(chunk) - c; let t = sum + y; c = (t - sum) - y; sum = t; } sum + (sum_block(chunks.remainder()) - c) } AlgorithmThroughputMean absolute error sum_orlp112.2 GB/s1.2306 You can of course tweak the number 256, I found that using 128 was $\approx$ 20% slower, and that 512 didn’t really improve performance but did cost accuracy. Conclusion I think the fadd_algebraic and similar algebraic intrinsics are very useful for achieving high-speed floating-point routines, and that other languages should add them as well. A global -ffast-math is not good enough, as we’ve seen above the best implementation was a hybrid between automatically optimized math for speed, and manually implemented non-associative compensated operations. Finally, if you are using LLVM, beware of -ffast-math. It is undefined behavior to produce a NaN or infinity while that flag is set in LLVM. I have no idea why they chose this hardcore stance which makes virtually every program that uses it unsound. If you are targetting LLVM with your language, avoid the nnan and ninf fast-math flags.