The Riemann Hypothesis (RH) is the conjecture that all the zeros of the Riemann zeta function ζ(_s_) in the critical strip, i.e. the region of the complex plane with real part between 0 and 1, have real part equal to ½. The Quasi Riemann Hypothesis (QRH) says that there exists a constant θ < 1 such that no zeros of ζ(_s_) have real part greater than θ. OpenAI has published a paper claiming QRH with θ = 7/8. The RH is so important to number theory that even partial results can have big consequences. This post will focus on one consequence: the error term in the Prime Number Theorem. The Prime Number Theorem says that π(_x_), the number of primes less than _x_ , is asymptotically equal to Li(_x_). We’d like to know more specifically at what rate π(_x_) approaches Li(_x_). The best known result before the QRH announcement was If the QRH holds for some θ, such as OpenAI’s assertion that θ = 7/8, If RH holds, θ = ½. Incidentally, you may have seen the Prime Number Theorem stated with _x_ /log(_x_) rather than Li(_x_). These two functions are asymptotically equal, so they give the same theorem, if you’re not interested in quantifying the rate of convergence. The function Li(_x_) gives better error bounds.

The Fast Fourier Transform (FFT) algorithm can compute the discrete Fourier transform of a sequence of length _n_ in time _O_(_n_ log _n_). OpenAI recently posted a paper saying there is an algorithm that could compute the discrete Fourier transform in _O_(_n_ (log _n_)1 − ε) time for ε = 10−13. This result is amazing. It seemed that _O_(_n_ log _n_) was as good as you could do, which it provably is for sorting algorithms. The result is also of absolutely no practical value, for now. But since the theorem shows that our assumptions were wrong regarding what we thought was possible, however slightly, maybe we’re in for further surprises. Maybe the ε crack will grow. It wouldn’t be the first time.

For a real number _x_ , the irrationality index μ(_x_) is a way of measuring how well _x_ can be approximated by rational numbers. If _x_ is rational, μ(_x_) = 1. If _x_ is irrational, μ(_x_) ≥ 2. OpenAI recently published a proof that μ(π) = 2. Almost all real numbers have irrationality exponent 2, so the new result says π is typical in this regard. There are numbers proven to have irrationality index greater than 2 (more on that below), but π isn’t one of them. The irrationality exponent μ(_x_) is defined as the supremum of the set of values ν such that for infinitely many coprime integers _p_ and _q_ with _q_ > 0. This means that the approximation error for approximating π with a rational number _p_ /_q_ is typically on the order of 1/_q_ ², just like most irrational numbers. There are numbers with higher irrationality exponents. For example, Cahen’s constant _C_ has irrationality exponent 3. This means _C_ is an irrational number that has infinitely many rational approximations _p_ /_q_ with error less than 1/_q_ ³.

Gödel’s incompleteness theorem illustrated the need to distinguish between what is true and what is provable. There are true statements that cannot be proven. Let □ _p_ denote the assertion that _p_ is provable in Peano arithmetic. The logic with this interpretation for the □ operator is the Gödel-Löb logic, also called provability logic. This is a normal modal logic with the additional axiom □(□ _p_ → _p_) → □ _p_ , known as Löb’s axiom. A couple days ago I wrote about topological models for modal logic. Is there a topological model for Gödel-Löb logic? There is, but it’s not quite the same construction as in the previous post. A topological model of Gödel-Löb logic associates _p_ with a set _P_ and ◇ _p_ with the **derived set** of _P_ rather than its closure. The difference between the closure of _P_ and the derived set of _P_ is subtle, but important to this discussion. The closure of a set _P_ is the union of _P_ and all of its limit points. The derived set of _P_ is the set of limit points of _P_. The distinction is that not every point of _P_ is necessarily a limit point of _P_. A point _x_ is a limit point of _P_ if every open set containing _x_ contains a point of _P_ _in addition to x itself_. A topological space _X_ that models Gödel-Löb logic must be **scattered** , meaning that every open set must contain an isolated point, a point with no limit points. For example, consider _X_ = {0} ∪ {1, ½, ⅓, ¼, …} with the topology inherited from the ordinary topology on the real line. Then every point except 0 is isolated, and every open set contains isolated points. A statement in Gödel-Löb logic is true if its topological interpretation holds for **all** scattered spaces.

The previous post discussed a superficial connection between modal logic and topology, that both give use the terms _regular_ and _normal_ to indicate added sets of axioms. McKinsey and Tarski developed a deeper connection between modal logic and topology that we’ll discuss here. Starting with a topological space _X_ and a proposition _p_ , define [[_p_]] as the set of points in _X_ at which _p_ is true. Define □ _p_ to be true at points in the interior of [[_p_]] and define ◇ _p_ to be true on the closure of [[_p_]]. You could think of □ _p_ as the points where _p_ is robustly true. Not only is _p_ true at _x_ , there’s some wiggle room around _x_ , i.e. an open set, in which _p_ remains true. You could think of ◇ _p_ as there points where we cannot rule out the possibility of _p_ being true using open sets. If ◇ _p_ includes _x_ , any open set containing _x_ also contains part of ◇ _p_ , though it may also contain points outside of ◇ _p._ ## Regularity For any topology on _X_ , the logic constructed above is regular. The axiom ◇ _p_ ⇔ ¬ (□ ¬ _p_) holds because the interior of a set is the complement of the interior of its complement [1]. Note that this is a regularity result for the modal logic, not the topology. The topology could be arbitrary, and not necessarily regular or normal in the topological sense. ## S4 The logic constructed above also satisfies a couple more axioms. We have □ _p_ → _p_ because the interior of a set is a subset of the set, and □ _p_ → □□ _p_ because the interior of the interior of a set is simply the interior. This means the modal logic corresponding to a topology satisfies the S4 axioms. You could say S4 is the logic that corresponds to the McKinsey and Tarski logic of all topological spaces. ## More logics and more topologies So S4 is the logic that corresponds to _all_ topologies. We could look at more restricted topologies and ask what are their corresponding logics. Or we could start with a modal logic and ask whether there’s a topology that models that logic. Interesting logics correspond to badly behaved topological spaces. Familiar topological spaces like the real line correspond to S4. ### Trivial modal logic The discrete topology corresponds to the trivial modal logic. All sets are open, and closed, so the any set is the same as its interior and its closure. So □ _p_ and ◇ _p_ reduce to just _p_. ### S5 For the indiscrete topology, □ _p_ corresponds to a proposition holding everywhere and ◇ _p_ corresponds to it holding somewhere. If the topological space has infinitely many points, the corresponding modal logic is S5. [2] ### Between S4 and S5 The cofinite topology on an infinite set _X_ defines a set _U_ to be open if the complement of _U_ is finite. The McKinsey-Tarski logic of the cofinite topology is somewhere between S4 and S5. You can show that the formula _p_ ∧ ◇□ _p_ → □ _p_ holds, which doesn’t hold in S4, and the formula ◇ _p_ → □◇ _p_ does not hold, though it must hold in S5. ## Related posts * Modal logic and science fiction * Naming and numbering modal logics * Modal logic and cybersecurity [1] We should also verify that if _A_ ∩ _B_ ⊂ _C_ , then Interior(_A_) ∩ Interior(_B_) ⊂ Interior(_C_). [2] Propositions can only have a finite number of terms. Having infinite points in the topological space prevents the corresponding logic from proving theorems that don’t necessarily hold in S5.

Suppose you want to solve a triangle. You know two sides and the angle between them. Then you can solve for the third side using the law of cosines. Now suppose you want to solve a **big** triangle, a triangle on the surface of the earth so large that the curvature of the earth matters. You can still use the law of cosines, but you’ll need the spherical law of cosines: cos(_c_) = cos(_a_) cos(_b_) + sin(_a_) sin(_b_) cos(_C_). If you know the (angular) lengths of sides _a_ and _b_ , and (tangential) angle _C_ between the two sides, you can solve for _c_ by taking the inverse cosine of the right hand side above. Now suppose you want to solve this big triangle because you’re a **navigator** on a ship a couple centuries ago, doing calculations by looking up trig functions and inverse trig functions in a table. You’re interested in triangles that are so big that you have to account for the fact that you’re living on a sphere. But at the same time, you’re triangles are still fairly small relative to the size of the globe. ## The problem with the law of cosines The numbers _a_ and _b_ will often be fairly small, and so their cosines will be near 1 and their sines are near zero. So the calculation cos(_a_) cos(_b_) + sin(_a_) sin(_b_) cos(_C_) will add a number near 1 and a number near zero. That’s a problem. Say you’re working with five decimal place arithmetic. Then if the second term above is less than 10−5, its contribution to the sum gets completely lost in the addition to the first term. If the second term is larger than 10−5 but still small, its contribution to the sum will be partially lost. ## Law of haversines Enter the haversine, defined by hav(θ) = (1 − cos(θ))/2. The expression 1 − cos θ was called the versine, and so half of the versine is the haversine. In terms of the haversine, the law of cosines above becomes the law of haversines: hav(_c_) = hav(_a_ − _b_) + sin(_a_) sin(_b_) hav(_C_). Now suppose you have a table of haversines and inverse haversines. The law of haversines requires a little less work: you have one less table lookup, and you trade a product for a subtraction. But the primary advantage is numerical accuracy: the terms on the right side have roughly the same size. ## Tables Note that we’re assuming the values in your table of haversines have been calculated correctly to the given precision. If you calculated your own values of haversines from the definition above, you’d lose precision in the subtraction 1 − cos θ, defeating the advantage of the law of haversines [1]. ## History According to Wikipedia. > The first table of haversines in English was published by James Andrew in 1805, but Florian Cajori credits an earlier use by José de Mendoza y Ríos in 1801. The term _haversine_ was coined in 1835 by James Inman. ## Experiments I ran some experiments that carried out arithmetic in float16 (11 bits of precision) to approximate what someone might have done by hand. When the difference between _a_ and _b_ was on the order of 1° or 0.1°, the law of cosine method often overflowed: the right-hand side evaluated to something larger than 1 even though theoretically it should be less than 1. The haversine method never overflowed. The median error for the haversine method was a couple orders of magnitude less than that of the cosine method. ## Related posts * Lews & Clark navigation * How many trig functions are there? * Analog of Heron’s formula for a sphere [1] hav(θ) = (1 − cos(θ))/2 = sin²(θ/2). If you calculated hav θ by looking up sin(θ/2) and squaring it, you’d be doing extra work, but you wouldn’t have numerical problems.

Nothing grows exponentially forever. What appears to be an exponential curve often turns out to be some sort of S curve, such as a logistic curve. Suppose you’re collecting data on the left side of the curve. If there’s even a small amount of error in your data, you won’t be able to predict the asymptotic value with any accuracy. But if you have data on both sides of the inflection point, you can make a good prediction of the limiting value. I’ve written about this before, explaining that the problem is hard, but I didn’t say _why_ it’s hard. Here I’d like to give an idea why it’s hard. Suppose you want to fit a logistic equation to three distinct values of _t_ and the corresponding values of _y_. There is a unique solution, but in general you cannot find a solution in closed form. However, if the values of _t_ are evenly spaced there is a method [1] to solve for the parameters _L_ , _k_ , and _t_ 0. For this post we’re only interested in the limiting value _L_ , and it can be found by independent of _h_. To find out how small changes in the _y_ ‘s change the estimate of _L_ , we take the partial derivatives of _L_ with respect to the _y_ ‘s and find and All three derivatives have the same expression in the denominator: _y_ 1² − _y_ 0 _y_ 2. If the function _y_(_t_) were an exponential, this expression would be exactly zero [2]. The function _y_(_t_) is not exactly exponential, but it is _approximately_ exponential when the _t_ ‘s are in the left or right tail of the logistic curve. The further out in either tail the _t_ ‘s are, the closer the expression is to zero. So when all the _t_ ‘s come from the same side of the inflection point, _y_(_t_) is nearly exponential the partial derivatives are huge and so the fitted value of _L_ is extremely sensitive to changes in the _y_ ‘s. [1] Raymond Pearl and Lowell J. Reed. On the Rate of Growth of the Population of the United States Since 1790 and its Mathematical Representation. Proceedings of the National Academy of Sciences of the United States of America, Vol. 6, No. 6 (Jun. 15, 1920), pp. 275-288 [2] exp(_x_ + _h_)² = exp(_x_)² exp(_h_)² = exp(_x_) exp(_x_ + 2 _h_)

The product of four consecutive Fibonacci numbers equals the product of two consecutive integers. For example, 3 × 5 × 8 × 13 = 39 × 40. I ran across this theorem in a note [1] that says “The product of any four consecutive Fibonacci numbers is twice a triangular number.” Since triangular numbers have the form _n_(_n_ + 1)/2, twice a triangular number is the product of two consecutive integers. The note also gives a way to find the numbers on the right hand side. We have _F_ _n_ _F_ _n_ +1 _F_ _n_ +2 _F_ _n_ +3 = _m_(_m_ + 1) where _m_ equals _F_ _n_ +1 _F_ _n_ +2 if _n_ is odd and _F_ _n_ _F_ _n_ +3 if _n_ is even. In the example at the top, 3 is the 4th Fibonacci number, so _n_ = 4. Since 4 is even, _m_ is the product of the 4th and 7th Fibonacci numbers, i.e. _m_ = 3 × 13 = 39. ## More Fibonacci posts * Fibonacci meets Pythagoras * Certified Fibonacci numbers * Turning trig identities into Fibonacci identities [1] K. B. Subramaniam. On a link between Triangular and Fibonacci numbers. The Mathematical Gazette, Vol. 103, No. 558 (November 2019), p. 489.

Comments

Suppose I give you an _n_ -digit number and ask you what it represents. This seems impossible, and in theory it _is_ impossible. But in practice it’s often possible. Apps on a phone may automatically interpret a 10-digit number as a phone number or a 16-digit number as a package tracking number. And very often these interpretations are correct, given the kinds of things most people use their phones for. It’s not surprising that a 10-digit number _on a phone_ is a _phone number_. It’s more interesting that a 16-digit number is likely a tracking number. It could be other things, such as a credit card number. But people don’t usually write out credit card numbers in a text note; credit card numbers likely saved in some more opaque way. I run into a variation of this problem routinely, trying to infer what a number represents inside medical notes. A five-digit number could be a US postal code, or it could be a medical procedure code. A six-digit number could be a date in MMDDYY format, or it could be a medical record number. A ten-digit number could be a phone number, or it could be an NPI (National Provider Identifier) number. It’s interesting that it’s possible make a good guess at what a number means inside unstructured text. Context has been lost, but not all context: you know you’re looking at medical notes. And that meager bit of context can be surprisingly useful.

Yesterday OpenAI announced a proof that settled a long-standing question about the Navier-Stokes equations from fluid dynamics. The announcement has created a lot of buzz, as one would expect. But there’s an aspect of OpenAI’s work that I haven’t seen anyone talk about: they posted a Lean 4 formal proof at the same time as their conventional human-readable proof. Quite a few other mathematical conjectures have been settled recently using AI, and these have also been accompanied with formal proofs, using Lean 4 in particular. Until very recently, generating machine-verifiable formal proofs has been **excruciatingly tedious**. In 2005, Henk Barendregt and Freek Wiedijk wrote > To give an indication of how much work is needed for formalisation, we estimate that it takes approximately one work-week (five work-days of eight work-hours) to formalise one page from an undergraduate mathematics textbook. That was the rule of thumb: **forty hours per page**. And this in the context of undergraduate textbooks. Research publications are much denser than textbooks. Furthermore, page 100 of a textbook probably depends mostly on material on pages 1 through 99. A sentence in a research article could cite anything that has been published before. Say a research article takes 20 times more effort to formalize than page in an undergraduate textbook. Then formalizing the 166-page paper from OpenAI would take 132,800 person-hours. It took OpenAI 17 hours to verify their proof in Lean. I hesitate to use the word “revolutionary,” but lowering the cost of anything by **four orders of magnitude** is revolutionary. I’ve used AI to generate formal proofs to check my work just for a little blog post. I wouldn’t dream of doing that if I had to pay someone a week’s salary to check my work. Formal verification doesn’t just apply to mathematics. You could, for example, formally verify that a set of security policies are consistent and that, given certain assumptions, they accomplish their purpose. You could formally verify that a smart contract imposes a certain maximum liability. You could verify the correctness of mission-critical algorithms. These problems are easier than formalizing mathematics research, and it is easier to quantify the return on investment. ## Related posts * Automation and validation * Formal methods let you explore the corners * When are formal methods worth the effort?

There are rumors that a long-standing math problem, one of the Millennium Prize problems, has been solved. The problem concerns technical properties of solutions to the **Navier-Stokes equations** [1], an equation that describe the dynamics of fluid flow. Popular accounts of the problem are often oversimplified and misleading. Some reports will speak of the problem as “solving the Navier-Stokes equations.” The task is not to write down a closed-form solution, which can’t be done, or solve the equations numerically, which has been done for decades. The problem is to prove theoretical properties of solutions which are of little interest in practice. There has been progress toward settling the Navier-Stokes problem. Terence Tao wrote a post on this yesterday. What’s also interesting is the intrigue around the possible solution. A post this morning says > If I am reading this correctly, Tristan Buckmaster is alleging OAI has a resolution of Navier-Stokes … which maybe used info from Buckmaster and Levent Alpöge’s private Codex sessions. Buckmaster asked OpenAI whether they used his private sessions and they have not responded. ## Personal note This topic connects parts of my career spanning decades. My graduate work was in PDEs and I had some interest in the Navier-Stokes equations. Here are some notes I wrote back in the day. Now I work more with privacy than with PDEs. The question of whether OpenAI uses private data, contradicting their stated policy, is more relevant to my current work than whether the Navier-Stokes equations have global regular solutions. ## Related posts * Engineering a waterpark * Euler-Lagrange equations * Data privacy [1] I never know whether to say equation or equations. You’ll hear both. You could think of Navier-Stokes as one vector-valued PDE or three scalar-valued equations.

The previous post discussed the motivation for and application of the rank-trace theorem. This post will give a proof. Suppose _A_ is a real symmetric matrix. The rank-trace inequality says where tr is the trace operator, the sum of the elements along the diagonal of the matrix. ## Terse proof Here’s the proof in a nutshell: diagonalize _A_ and use the Cauchy-Schwarz inequality. ## Detailed proof Now let’s unpack that. Any real symmetric matrix _A_ is similar to a matrix _D_ with the eigenvalues of _A_ along the diagonal. The trace of a matrix stays the same under a similarity transformation, i.e. multiplying by _P_ on one side and its inverse on the other side. So without loss of generality we may as well assume _A_ is diagonal. The rank of a matrix equals the number of non-zero eigenvalues, so a vector containing the non-zero eigenvalues of _A_ has length _r_ where _r_ is the rank of _A_. Define _w_ to be the vector of dimension _r_ consisting of all 1’s. Then by the Cauchy-Schwarz inequality we have ## Cyclic trace property Why should a matrix _A_ and its diagonalization _D_ have the same trace? The trace of a matrix product _AB_ equals the trace of the product _BA_. To prove this, write out matrix products and the traces, then note that the two expressions are equal. Therefore More generally, trace has the cyclic property However, not all permutations preserve the trace. For example, let Then but

Suppose you want to know the rank of an _n_ × _n_ matrix _A_ , the number of linearly independent rows of _A_ , or equivalently the number of linearly independent columns. There are at least three difficulties. ## Difficulties in computing rank First of all, rank is not a continuous function of a matrix. Since rank is an integer, an arbitrarily small change in the matrix could cause a discrete change in the rank [1]. A small error in computing _A_ could produce a matrix with a different rank. Second, finding the rank takes _O_(_n_ ³) operations, which may or may not be an issue depending on context. Third, you may not have the matrix _A_ in an explicit form. Maybe you’re able to compute products _Av_ for vectors _v_ but it’s not practical to form the entire matrix _A_. ## Rank-trace inequality If you don’t need to know the rank of _A_ per se, but only need to know whether it is above a certain size, a lower bound on the rank may enough. Suppose _A_ is a Hermitian matrix. If _A_ is real, this means _A_ is symmetric. If _A_ is complex, this means _A_ equals its conjugate transpose. Then the rank-trace inequality says \operatorname{rank}(A)\ge\frac{(\operatorname{tr} A)^2}{\operatorname{tr}(A^2)} The quantity on the right hand side is known as the **stable rank** of _A_. It’s not a rank in any algebraic sense, but it gives a lower bound on rank. And it solves the three problems listed above. ### Stability First of all, trace _is_ a continuous function of a matrix, and so stable rank is also a continuous function of a matrix, provided the denominator isn’t zero. A small change to a matrix only makes a small change to its stable rank. That’s why stable rank is called stable. ### Efficiency Second, although computing rank takes _O_(_n_ ³) operations, computing stable rank takes only _O_(_n_ ²) operations, though this isn’t immediately obvious. The trace of _A_ takes _n_ operations: simply sum the elements on the diagonal of _A_. But how do you take the trace of _A_ ²? Squaring _A_ takes _n_ ³ operations, and so if you had to square _A_ to find the trace of _A_ ² the rank-trace inequality would have no efficiency advantage over finding the rank of _A_. But you can compute the trace of _A_ ² via \operatorname{tr}(A^2) = \sum_{i=1}^n \sum_{j=1}^n |a_{ij}|^2 ### Formation Now suppose you don’t have the matrix _A_ per se but you do have a way of probing _A_ , computing the product of vectors with _A_. Maybe _A_ is too large to fit into memory, or explicitly computing the elements of _A_ would take too long. There are Monte Carlo algorithms for estimating the traces of _A_ and _A_ ² that could be used together to estimate the stable rank of _A_. ## Demonstration The following Python code illustrates the discussion above. import numpy as np np.random.seed(20260904) n = 5 B = np.random.randn(n, n) A = B.T @ B + 1e-8 * np.eye(n) # Gram matrix plus a tiny shift => SPD rank_A = np.linalg.matrix_rank(A) tr_A = np.trace(A) tr_A2 = np.trace(A @ A) # matrix product sum_sq = np.sum(A * A) # element-by-element product stable_rank = (tr_A ** 2) / tr_A2 print(f"A =\n{A}\n") print(f"rank(A) = {rank_A}") print(f"tr(A) = {tr_A:.12f}") print(f"tr(A^2) direct = {tr_A2:.12f}") print(f"tr(A^2) indirect = {sum_sq:.12f}") print(f"stable rank = {stable_rank:.12f}") The code above produces the output below. A = [[ 1.09945682 0.4899665 0.98901845 0.66983113 -1.35006341] [ 0.4899665 0.98531254 0.35067791 0.89757603 -0.72037507] [ 0.98901845 0.35067791 4.31233926 0.94556225 -0.54819048] [ 0.66983113 0.89757603 0.94556225 1.3494295 -1.33840786] [-1.35006341 -0.72037507 -0.54819048 -1.33840786 3.54858332]] rank(A) = 5 tr(A) = 11.295121449420 tr(A^2) direct = 51.035447533673 tr(A^2) indirect = 51.035447533673 stable rank = 2.499826585688 [1] Topological argument: A map from a connected space (such as ℝ _n_ × _n_) onto a discrete space (such as ℤ) cannot be continuous, otherwise the inverse images of the points in the range would partition the connected space into disjoint open sets, violating the definition of a connected space.

Eric Lu announced on X today that he has factored RSA-260, a number _N_ with 260 digits (862 bits) that is the product of two large primes [1]. RSA numbers are challenge problems posed to gauge the security of RSA encryption, which rests on the difficulty of factoring large numbers [2]. The naming scheme is confusing because RSA-_n_ might have _n_ digits or _n_ bits. For example, RSA-768 is smaller than RSA-260 because the former has 768 bits and the latter has 260 digits. RSA-260 is the largest RSA number factored so far. What does the news of its factorization say about the security of RSA? Based on equations here, an RSA key with 862 bits would have a security level of 74 bits, i.e. the same security level as symmetric encryption with a 74-bit key. The minimum recommended RSA key size now is 2048 bits, which has a security level of 107 bits. Security levels are on a logarithmic scale: each additional bit of security doubles the effort required to break the encryption by brute force. So breaking a 2048-bit RSA key would take 234, roughly 1010, times more effort than factoring RSA-260. All this depends on numerous assumptions, such as the state of factorization algorithms and the non-existence of CRQC [3]. ## Related posts * RSA implementation flaws * Generating and inspecting an RSA key * Martin Gardner’s RSA article * RSA munitions T-shirt [1] _N_ = _pq_ = 22112825529529666435281085255026230927612089502470015394413748319128822941402001986512729726569746599085900330031400051170742204560859276357953757185954298838958709229238491006703034124620545784566413664540684214361293017694020846391065875914794251435144458199 _p_ = 4397328654844826923795068102505872571721883526553349659561256924505973939597593482272505698004801207988043088656411102133523080581 _q_ = 5028695206842569864686141618253083416610081090075366674776775706538324961364412200138116378509733307971876652984898985905923678379 [2] The ability to efficiently factor large primes would break RSA. It’s possible that there’s a way to break RSA without being able to factor large numbers. More on that here. [3] Cryptographically-relevant quantum computer. Quantum computers exist, but so far they’re cryptographically irrelevant. So far quantum computers cannot factor 21 without cheating.

This post provides a couple examples to go along with two earlier posts. The pattern we’re illustrating is families of polynomials _p_ _n_(_x_) that each satisfy a differential equation and a three-term recurrence. The differential equations have a second solution _q_ _n_(_x_) that is the larger solution with respect to _x_ but the smaller solution with respect to _n_. In both the examples below _p_ _n_(_x_) is a polynomial, and so bounded on the interval [−1, 1], and _q_ _n_(_x_) is not a polynomial, with singularities at ±1. This is analogous to the previous examples with Bessel functions _J_ _n_(_x_) and _Q_ _n_(_x_) that satisfy the same differential equation but have contrasting behavior with respect to _x_ versus _n_. ## Legendre polynomials The differential equation has two solutions for each _n_ , _P_ _n_(_x_) and _Q_ _n_(_x_). The solutions _P_ _n_(_x_) are the Legendre polynomials. The solutions _Q_ _n_(_x_) are not polynomials but involve a term log((1 + _x_)/(1 − _x_)) that blows up at 1 and −1. But for fixed _x_ and increasing _n_ , _P_ _n_(_x_) grows exponentially and _Q_ _n_(_x_) decays exponentially. ## Chebyshev polynomials The differential equation has two solutions for each _n_ , _T_ _n_(_x_) and _V_ _n_(_x_). The solutions _T_ _n_(_x_) are the Chebyshev polynomials. The solutions _V_ _n_(_x_) are not polynomials but involve a term √(1 — _x_ ²) that become vertical up at 1 and −1. But for fixed _x_ and increasing _n_ , _T_ _n_(_x_) grows exponentially and _V_ _n_(_x_) decays exponentially.

When you’re interested in studying a family of functions, it can be useful to look at a differential equation that the functions solve. This is a theme I’ve written about several times, most recently here and here, but also three years ago here. Orthogonal polynomials are mathematically elegant as well as very useful in applications [1]. Various families of orthogonal polynomials satisfy various differential equations. These equations have a polynomial and non-polynomial solutions. What use are the latter? If the differential equation modeled something physical, then the second solution would be necessary to have a complete basis of solutions. But if the differential equation is only instrumental in studying the orthogonal polynomials, what use is a non-polynomial solution? These non-polynomial solutions turn out to be useful. Just as “junk” DNA turned out not to be junk, these “junk” solutions are important. Junk DNA doesn’t directly code for proteins, but it regulates DNA that does code for proteins and serves other purposes. Similarly, these non-polynomial solutions carry information related to the polynomial solutions. For example, orthogonal polynomials are used to construct numerical integration methods, such as Gaussian quadrature, and the associated non-polynomial solutions describe the error in these integration methods. Incidentally, Gaussian quadrature is based on Legendre polynomials, mentioned in the previous post. For every family of orthogonal polynomials there is a corresponding integration method. See these notes. Another tie-in to recent posts is that these non-polynomial solutions are the minimal solution to the polynomial family’s three-term recurrence, the solution that takes extra care to compute numerically. This post has been very high-level, alluding to ideas without going into details. I’d like to write future posts that go into more depth regarding the ideas introduced here. [1] “Real analysts cannot do without Fourier, complex analysts cannot do without Laurent, and numerical analysts cannot do without Chebyshev [polynomials].” — Lloyd N. Trefethen”

When I hear the term _ultraspherical_ I think of something extremely spherical. For example, a baseball is spherical, but a billiard ball is more spherical. Maybe a highly polished billiard ball is ultraspherical. Using this line of thought, the term **ultraspherical polynomial** is inexplicable. This is an example of the arcane terminology I wrote about recently. In this post I’ll explain what it conveys. A spherical polynomial is a polynomial that naturally falls out of solving Laplace’s equation in spherical coordinates, using separation of variables. Legendre polynomials are spherical polynomials. **Gegenbauer polynomials** are so called because a man named Gegenbauer studied them, just as Legendre polynomials take their name from Legendre. Gegenbauer polynomials are also called ultraspherical polynomials. Why is that? There are two possible reasons. I’m not sure which is the historical reason, but both are plausible and are useful mnemonics. Ultraspherical polynomials are not extremely spherical, they’re _beyond_ spherical in some sense. More modern terminology uses the hyper- prefix rather than ultra-, which helps a bit. Ultraspherical polynomials are beyond spherical in two ways. Gegenbauer polynomials are a generalization of Legendre polynomials, so they’re beyond Legendre polynomials in this sense. More importantly, Gegenbauer polynomials fall out of solving Laplace’s equation on a hypersphere, i.e. a sphere in ℝ _n_ for _n_ > 3, just as Legendre polynomials fall out of the case _n_ = 3. It makes sense to call these polynomials **hyperspherical** because they fall out of solving an equation on a hypersphere. Unfortunately the classical term is _ultraspherical_ rather than _hyperspherical_. I think, but I’m not sure, that at one time higher dimensional spheres were called hyperspheres, but the the higher dimensional analog of spherical coordinates was called ultraspherical coordinates. If so, it would be understandable that the adjective modifying _coordinates_ would be applied to the polynomials that result from solving equations in these coordinates.

The previous post gave several examples of three-term recurrence relations for special functions. These relations can be computationally useful, but they have to be applied carefully. Several years ago I wrote a post on stable and unstable recurrences. In that post I show that the stability of the recurrence relation for Bessel functions produces depends on which kind of Bessel function and which direction the recurrence is applied. In the forward direction, computing higher order values from lower order values, works well for Bessel functions of the second kind _Y_ _n_ but not for Bessel functions of the first kind _J_ _n_. In the reverse direction, the recurrence is stable for _J_ _n_ but not for _Y_ _n_. I didn’t explain in that post why this is. In this post I will. Second order linear difference equations have two independent solutions, just like second order linear differential equations. For both kinds of equations, all solutions are linear combinations of the two solutions. Suppose one solution grows with _n_ and the other decays. You may want to compute the decaying solution, but in doing so you might pick up a small component of the growing solution due to rounding error. This post illustrates this phenomena for differential equations, and this post illustrates it for difference equations. When you look at a plot of Bessel functions in a text book, you’ll probably see a few plots of _J_ _n_(_x_) and _Y_ _n_(_x_) for a few small values of _n_. The functions seem to behave roughly the same way, like sine and cosine. And that’s true, **as functions of _x_**. But it’s not true for _J_ _n_(_x_) and _Y_ _n_(_x_) as functions of _n_ for fixed _x_. As _n_ increases, _J_ _n_(_x_) decays to zero and _Y_ _n_(_x_) goes off to −∞. That’s the source of numerical instability. And there will be similar instability problems for other recurrences where the ratios of the two independent solutions goes to zero or infinity as a function of _n_. There are techniques for computing the solution that does not diverse, the so-called minimal solution, such as Miller’s algorithm mentioned here.

Probability density function must integrate to 1, and so if you know a density function up to a constant, the constant is determined. When you’re looking at a probability density _f_(_x_) for the first time, it helps to ignore the normalizing constant. Concentrate on the part of the function involving _x_ and know that the normalizing constant is whatever it has to be. For example, about half of the ink that it takes to write down a beta or chi-squared density is devoted to the normalization constant; the rest of the expression is easier to understand. This post will do the opposite of the advice above and focus on normalization constants because this ties into the previous post on modified Bessel functions. The **von Mises** probability distribution on a circle has two parameters, μ and κ, and its density function is The normalizing constant is 2π _I_ 0(κ). The factor of 2π is unsurprising for anything defined on a circle. The more interesting part is _I_ 0, the modified Bessel function of order 0. The **von Mises-Fisher** distribution is the generalization of the von Mises distribution to a sphere in _p_ dimensions. The density function is where the normalization constant _C_ _p_(κ) is where _I_ _p_ /2 − 1 is the modified Bessel function of order _p_ /2 − 1. The values of **x** and **μ** are in bold face because they are now vectors, points on the unit sphere. When _p_ = 2, we have the “sphere” in two dimensions, i.e. the circle, and the von Mises-Fisher distribution reduces to the von Mises distribution. But where did the cosine go? The inner product of **x** and **μ** is the cosine of the angle between the two vectors. When _p_ = 3, obviously an important special case, the von Mises-Fisher distribution is known as the **Fisher** distribution. In that case the normalizing constant _C_ 3(κ) can be written without using modified Bessel functions because when ν = ½ + _n_ for an integer _n_ , _I_ ν(_x_) is an elementary function.

There many examples of families of functions where each function can be computed as a linear combination of the two previous terms where _a_ and _b_ are functions of _x_ but not on _n_. This is called a three-term recurrence formula. It’s amazing how often you can run into three-term recurrence formulas. There are theorems that give conditions for such recurrences to hold, but I haven’t reached the bottom of that rabbit hole [1]. For this post I just want to give examples. **Bessel functions** of the first and second kind: **Modified Bessel functions** of the first and second kind: **Chebyshev polynomials** of the first and second kind: **Hermite polynomials** (physicists’ convention): **Legendre polynomials** : [1] See Bochner’s theorem for orthogonal polynomials, the Nikiforov–Uvarov method, and Infeld-Hull factorization.

Special functions often have arcane names that not very helpful without some context. The previous post goes into some reasons for this. This post will expand on a point at the end of the post about “modified” functions. Things are given their names for a reason. Discovering that reason helps you understand their motivation and use. ## Pure math perspective For each integer _n_ , the modified Bessel function _I n_ is essentially the Bessel function _J n_ evaluated along the imaginary axis. Specifically, From a certain shallow perspective, that’s the end of the story: modified Bessel functions are modified in the sense that the argument is multiplied by _i_. And there’s a fiddly constant term up front for no apparent reason. But of course that’s not the end of the story or else this wouldn’t be worth an entire post. The equation above is analogous to the relationships between circular and hyperbolic functions These relationships are interesting because the circular and hyperbolic functions are independently meaningful. If you view these equations merely as definitions you lose their significance. Circular and hyperbolic functions were widely used before Euler discovered the connection between them. Similarly, there’s a reason the modified Bessel functions were given a name their own. If you were led to Bessel functions and modified Bessel functions separately by different applications, you would regard the equation as a **discovery** rather than just a definition. The following section explains why someone would be interested in modified Bessel functions. Before we move on, I’d like to explain the reason for the term _i_ − _n_ term. In general for all real ν. The reason for the exp(νπ _i_ /2) term is that it makes _I_ ν(_x_) real for all real _x_. ## Applied math perspective Bessel functions often arise from solving problems with **radial symmetry**. Solving the **wave equation** in cylindrical coordinates using separation of variables leads to Bessel’s differential equation and its solutions _J n_ and _Y n_, Bessel functions of the first and second kind. Solving the **heat equation** in cylindrical coordinates with separation of variables leads to the _modified_ Bessel equation and its solutions _I n_ and _K n_, the _modified_ Bessel functions of the first and second kind. This is the reason behind the complex analysis perspective above: the change of variables sending _x_ to _ix_ changes the sign of the _x_ ² term in Bessel’s equation. Bessel functions describe radially symmetric **oscillations** , such as the vibrations of a drum head. Modified Bessel functions describe radially symmetric **exponential decay** , such as in the heat in a cylinder. ## Other modified functions **Struve functions** are closely related to Bessel functions. The (modified) Struve functions also satisfy Bessel’s (modified) differential equation, but with a non-zero right hand side. The modified Struve functions are proportional to the unmodified Struve functions evaluated along the imaginary axis, with a proportionality constant that makes the modified Struve functions real for real arguments. There’s a similar relationship between the **Mathieu functions** and modified Mathieu functions. The general pattern is that “modified” in the context of special functions means “evaluated at _ix_ and multiplied by a constant to make the function real for real arguments.”

Special functions are special because they’re useful. They can also be shrouded in arcane terminology. These two facts are related. The more widely useful a function is, the more likely it is that the function will be discovered independently multiple times. Independent discoveries lead to varying definitions and notations. For example, there are two widely used definitions of Hermite polynomials, one used in probability and another used in physics, that only differ by a scaling factor. This also explains why there are so many variations on the definitions of the Fourier transform and spherical coordinates. Special functions were discovered and applied before they were studied systematically. As with most mathematics, practice preceded theory. In hindsight, some names and conventions were less than ideal, at least from the perspective of someone seeking to organize a theory. Functions can have arcane names for several reasons, one being that their usefulness became apparent long ago. In this case the odd terminology is evidence that the topic is worth knowing about. Sometimes special functions have bland, uninformative names because the names stuck before anybody could think of something better. Bob looks into an interesting family of functions, then later he finds another interesting family of functions. These become known as “Bob’s functions of the first kind” and “Bob’s functions of the second kind.” These names are quite understandable at the time, though in the future people will want to know what distinguishes the functions, other than the fact that Bob discovered them, and what the groupings have in common other than the order in which Bob found them. I started this post intending to discuss modified Bessel functions and explain what exactly is modified about them, but my preface became its own post. “Modified” is an example of the bland terminology mentioned above. There are Bessel functions and modified Bessel functions. Without more context, the “modified” term isn’t very informative. But it does provide a clue that there’s some kind of close relationship between the modified and unmodified functions. That’ll be the topic of my next post.

There’s a pizza shop near my home with a sign out front that says “Coming Soon.” When I drove by it this morning I thought about how you would model the time until an event happens that is “coming soon.” Suppose I look at the sign one day and guess how many days until the pizza shop will open. When I drive by a week later and guess again, should my guess be smaller? You might argue that the shop will open some day, fixed in time but unknown to me, and so every day I’m one day closer to the eventual opening. You might model the pizza shop opening like radioactive decay and say that the estimated number of days until it opens is always the same until the day it actually opens. Now I think this shop has been “coming soon” for over a year. So instead of decreasing, every day I increase my estimate of the time until the shop opens. Something has gone wrong that the owners didn’t expect when they put up the sign. Maybe the reasonable thing would be for estimated days until opening to decrease over time, but only up to a point. After some point, the longer a business has been “coming soon” the less like that it is coming soon, or coming at all. This brings up an interesting point about modeling. There are two probability distributions at work: the probability that the shop will eventually open, and the time until opening assuming it eventually opens. When the sign first goes up saying the business is coming soon, there’s some change that it is in fact not coming. Maybe you’re optimistic and think this probability is small, but it would seem unreasonable to think the probability is zero. That means the _expected_ number of days until opening is always infinite. If there’s a probability ε that the shop never opens, the expected time to opening is ε × ∞ + (1 − ε) × something = ∞.

A few weeks ago I wrote about the number system used in labeling spreadsheet columns. Labels run from A through Z, then AA through AZ, etc. This looks a lot like base 26, but it’s not quite the same. It has no analog of zero. If Z were like zero, Y would be followed by AZ. The Excel labeling system is not base 26, but what’s called bijective base 26. If you found fragments of writing from an ancient culture and inferred that five symbols were used as digits, how could you distinguish base 5 from bijective base 5? Suppose you believe these five symbols were digits ★ ☂ ☘ ☗ ☢ but you don’t know in what order. You just see sequences like ☂☘☢ and ★★☂ and believe they’re numbers. If you noticed that numbers often contain ☘, but ☘ never appears at the beginning of a number, you might infer that ☘ is a zero. But this would take a fairly large sample. If you found only 20 numbers, for example, you could hardly conclude ☘ never appears at the beginning of a number just because it doesn’t come at the beginning of any number you’ve seen. Now suppose you’ve found writing with more number symbols. Say you’ve found 17 numeric symbols. You might infer that the writing used a base 20 system, because it would be hard to imagine a human culture using base 17. Now imagine you find more fragments and confirmed that indeed there are 20 numeric symbols. Approached as a purely statistical problem, you’d need a very large sample to infer what the digits correspond to and whether they use a base 20 or bijective base 20 system (or some other system). You’re best hope is to find numbers in some context where you know what number is being represented. If you knew somehow that some symbol corresponds to 20, then you’d know they didn’t use base 20 because base _b_ doesn’t have a single symbol for _b_. If you had a huge collection of numbers but no context, which is highly unlikely, you could use Benford’s law to infer the meaning of the number symbols: the most common leading digit is probably 1, the next most common is probably 2, etc. This is interesting to think about, but it seems much more realistic that a number system would be decoded by finding context, such as a list of consecutive numbers or numbers with known meaning.

Suppose you have a planet in an elliptical orbit around a star. The math is identical for any light object orbiting a heavy object, such as a moon or satellite orbiting a planet, but we’ll call the heavy object a star and the light object a planet. The center of the star is not quite the center of the orbit. The planet moves along an ellipse with the star at one focus of that ellipse. Let _a_ be the semi-major axis of planet’s orbit, the maximum distance from the center of the ellipse to a point on the ellipse. Then the distance of a focus to the center of the ellipse is _ae_ where _e_ is the eccentricity of the ellipse. This defines eccentricity. The center of earth’s orbit is between three and four solar radii away from the center of the sun [1]. The planet is farthest from the star when it is along the major axis of the ellipse on the opposite side as the star. The distance is then _a_ + _ae_ , the distance to the center plus the distance from the center to the star. The planet is closest to the star on the opposite side of its orbit. There the distance is _a_ − _ae_. In summary the maximum distance to the star is _a_(1 + _e_) and the minimum distance is _a_(1 − _e_). If you had to guess the _average_ distance between the planet and its star, _a_ would be a good guess since it’s the average of the maximum and minimum distance. And that’s a good approximation, provided _e_ is small. The mean distance over time is _a_(1 + ½ _e_ ²). See derivation. The average distance is greater than _a_ because the planet moves faster when nearest the star and slower when further from the star. The relative error in approximating the mean distance by _a_ is then ½ _e_ ². When _e_ is small, ½ _e_ ² is very small. For the earth’s orbit, _e_ = 0.01671, and so the approximation is off by around 0.014%. The eccentricity of Pluto’s orbit is 0.2488, and so in that case the approximation is off by about 3.1%. The eccentricity of a Molniya orbit, used by some Russian satellites, is 0.74 [2]. For such satellites the error in approximating the mean distance to earth as the semimajor axis is around 27%. [1] For earth’s orbit, _e_ = 0.01671, _a_ = 1.496×1011 m, and the sun’s radius is _r_ = 6.957×108 m. And so _ea_ / _r_ = 3.59. [2] An object in such a highly elliptical orbit will spend a long time at the far side of its orbit, i.e. over Russia. Sort of a poor man’s geostationary orbit.

The first post in the recent series of posts on Hadamard matrices describes a way of constructing new Hadamard matrices from two other Hadamard matrices by taking their Kronecker product. Starting with a Hadamard matrix _H_ 0 and a Hadamard matrix _G_ , you can construct a sequence of Hadamard matrices by _H_ _n_ +1 = _G_ ⊗ _H_ _n_ for positive integers _n_. This is known as the generalized Sylvester method. Let _p_ _n_ be the proportion of 1s in _H_ _n_ and let _q_ be the proportion of 1s in _G_. Then you can show that the recurrence holds _p_ _n_ +1 = _q_ _p_ _n_ + (1 − _q_)(1 − _p n_). You can solve the recurrence to show that lim _n_ → ∞ _p_ _n_ = ½ and so as the iterations proceed, the ratio of number of 1s to the number of −1s approaches 1. This doesn’t say anything Hadamard matrices in general, but it does apply to all Hadamard matrices created by repeatedly applying the generalized Sylvester method. If you set _G_ and _H_ equal to the matrix then _p_ 0 = _q_ = ¾. Then for _n_ = 1, 2, 3, …, 8 the values of _p_ _n_ are 0.625 0.5625 0.53125 0.515625 0.5078125 0.50390625 0.501953125 0.5009765625.

Error correcting codes are most simply described in terms of the errors they can certainly correct. For example, the Hadamard code used for the Mariner 9 probe to Mars encoded each 6-bit pixel to a 32-bit codeword in such a way that the original pixel could be recovered if no more than 7 bits were corrupted in transit. What is the _probability_ that a pixel could be repaired if corrupted? That depends on your probability model. We will assume that the probability of each bit being flipped is _p_ and that errors are independent. (Are errors independent, i.e. if a bit flips, is the next bit more or less likely to flip? That would depend on context.) It’s straight-forward to calculate the probability that 7 or fewer or fewer bits out of 32 flip; this is the cumulative distribution of a binomial random variable. The following Python code will return the probability of _k_ or fewer successes out of _n_ trials, each with probability of success _p_ : from scipy.stats import binom print(binom.cdf(k, n, p)) For example, if there is a 10% chance that each bit will flip, there’s a 98.8% chance that 7 or fewer bits out of 32 will flip. However this only gives a **lower bound** on the probability of correcting an error. If eight bits flip in transit, we cannot tell with certainty which codeword was sent, but that doesn’t mean all possibilities are equally likely. Here things get messier. For the Hadamard code mentioned above, there’s a 50-50 chance of being able to recover a pixel transmitted with 8 flipped bits in the corresponding code word. The probability of correct recovery gets smaller with more corruption, but it doesn’t go to zero. Now suppose you’re given a desired error recovery rate and have to determine what value of _p_ it can sustain. For example, someone might say they want a 98.8% chance of recovering a pixel correctly, and you could come back and say _p_ must be less than or equal to 0.1. This would be a conservative answer because as discussed above, _p_ = 0.1 gives a pixel recovery probability of something more than 98.8, though it’s messy to calculate how much more. You could solve for _p_ by trial and error, or you could use some more sophisticated math to compute _p_ directly. Given a probability _F_ , you can solve for _p_ such that the probability of up to _k_ successes out of _n_ trials using the inverse of the regularized incomplete beta function. from scipy.special import betaincinv p = 1 - betaincinv(n - k, k + 1, F) Calculating _F_ given _n_ , _k_ , and _p_ could be a homework exercise in an introductory probability course. Solving for _p_ given _F_ , _n_ , and _k_ either requires some numerical programming or special functions and so would be a more challenging problem. ## Related posts * Golay code used in Voyager * Vehicle Identification Number checksum * Credit card checksum

Hadamard matrices are in the news following the recent announcement of a newly discovered Hadamard matrix. I’ve written three posts on Hadamard matrices recently, one as a sort of introduction and two on applications: the error correcting code used in the Mariner 9 probe and constructing sphere packings. A Hadamard matrix is an orthogonal matrix with all entries equal to ±1. Jacques Hadamard conjectured that there exist Hadamard matrices of order 4 _n_ for all positive integers _n_. It’s necessary that the order be divisible by 4, and Hadamard conjectured that this is sufficient [1]. How could you compactly represent a Hadamard matrix? Since the entries are all either 1 or − 1 each entry could be represented by a single bit, and _n_ ² bits could store an _n_ × _n_ Hadamard matrix. But we can do better. ## Methodical matrices If the matrix can be produced by an algorithm, you only need to store the name of the algorithm and the argument to the algorithm. So, for a 1024 × 1024 matrix applied by iterating Sylvester’s algorithm could be stored by saying “Apply Sylvester’s algorithm 10 times” rather than storing a megabyte of data. Paley’s method can create a Hadamard matrix corresponding to every prime power. So you could determine a Paley type matrix by storing the prime and the exponent. Next in complexity would be hybrid algorithms, such as start with the Paley method applied to 376 and then apply Sylvester’s method 3 times. There are more methods of creating Hadamard matrices than Sylvester’s method and Paley’s method, though they’re harder to describe and parameterize. ## Sporadic matrices If a Hadamard matrix cannot be constructed using an algorithm, you can still store the matrix in fewer than _n_ ² bits. Since the rows are orthogonal, the last row of the matrix is determined by all the previous rows, up to sign. So you could store a Hadamard matrix using _n_(_n_ − 1) + 1 bits. Some Hadamard matrices are symmetric or skew. A symmetric matrix is determined by its diagonal and the elements above the diagonal. So a symmetric Hadamard matrix could be represented by _n_(_n_ + 1)/2 bits. A skew Hadamard matrix isn’t quite skew-symmetric. A matrix _M_ is skew symmetric if _M_ T = − _M_. This implies the diagonal elements are 0, and Hadamard matrices cannot contain 0s. A Hadamard matrix _H_ is called skew if _H_ + _H_ T = 2 _I_. This implies the diagonal elements are all 1s and the elements below the diagonal have the opposite sign of the elements above the diagonal. Since the elements on the diagonal are determined, a skew Hadamard matrix can be sotred using _n_(_n_ − 1)/2 bits. Incidentally, there is a conjecture that there exist skew Hadamard matrices of order 4 _n_ for all positive _n_. [1] There are Hadamard matrices of order 1 and 2, but larger orders must be divisible by 4.

Yesterday Levent Alpöge announced that he and his colleagues had discovered a new Hadamard matrix using Claude AI. That motivated a post I wrote this morning on how to construct Hadamard matrices. I mentioned in that post that these matrices arise in applications. This evening I gave an example, describing how NASA used a Hadamard matrix of order 32 to transmit photos from the Mariner 9 spacecraft in 1971. This post will give another application: **sphere packing**. Conway and Sloane [1] give a correspondence between binary codes and sphere packings that they call Construction A. Given an (_n_ , _M_ , _d_) binary code _C_ , center a sphere on a point _x_ if and only if _x_ is a codeword in _C_. Here (_n_ , _M_ , _d_) means an error correcting code that encodes _M_ bits of data as strings of _n_ bits, with a minimum Hamming distance between code words of _d_ , i.e. all codewords differ in at least _d_ bits. The previous post described how to create a (32, 6, 16) code by stacking a Hadamard matrix _H_ of order 32 on top of − _H_ and turning −1’s into 0’s. The analogous construction for a (8, 4, 4) Hadamard code gives _E_ 8, the densest packing in ℝ8. We start with the Hadamard matrix and obtain the matrix whose centers form the sphere packing. This doesn’t look like the E8 sphere packing as it is usually presented, but it’s isomorphic. [1] J. H. Conway and N. J. A. Sloane. Sphere Packings, Lattices and Groups. Springer. 1999.