Splitting _N_! into _N_ factors is trivial. It’s _defined_ as the product of _N_ integers, the integers 1 through _N_. Can you still split _N_! into _N_ factors if you put a minimum size requirement on the factors? You can’t split 3! into three factors unless one of the factors is 1, but you can split 4! into factors ≥ 2: 4! = 24 = 2 × 2 × 2 × 3 Define _t_(_N_) as the largest number such that _N_! can be expressed as the product of integers greater than or equal to _t_(_N_). So based on the discussion above, _t_(_N_) ≥ 1 for all _N_ , _t_(3) = 1 and _t_(4) = 2. The _t_(_N_) are given in OEIS sequence A034258. Last Thursday an article came out proving an asymptotic estimate for _t_(_N_). You can find a preprint on arXiv. An upper bound on _t_(_N_) would be the _N_ th root of _N_! because that would correspond to splitting _N_! into _N_ factors of equal size. But the _N_ th root of _N_! will not be an integer for _N_ > 1, so the bound would not be tight. If we apply Stirling’s approximation for factorial we get _t_(_N_) / _N_ ≤ 1/_e_ + _O_(log(_N_) / _N_). It was conjectured that for large _N_ _t_(_N_) / _N_ ∼ 1/_e_ − _c_ 0 / log(_N_) and the new paper proves _t_(_N_) / _N_ ∼ 1/_e_ − _c_ 0 / log(_N_) + _O_(1 / log1 + _c_ _N_). .

About 14% of English words can be spelled by stringing together chemical symbols. For example, _sarcasm_ can be written S Ar Ca Sm using the symbols for sulfur, argon, calcium, and samarium. How would you write a program to determine whether a word can be written by stringing together chemical element abbreviations ignoring case? Obviously this isn’t a practical problem, but it illustrates a practical problem. You could solve the problem with a regular expression: (h|he|li|be|b|c|n|o|f|ne|na|mg|al|si| … |og)+ This solution could either be very efficient or very inefficient, depending on the inner workings of the software used to match the regular expression. A DFA (deterministic finite automata) regular expression engine would find a match, or report failure, linear time, i.e. the execution time would be proportional to the length of the string being matched. A NFA (nondeterministic finite automata) engine that uses backtracking could take exponential time. The root of the problem is the ambiguity in dividing words into single letters or pairs of letters. Many two-letter chemical symbols are made of letters that are themselves chemical symbols. For example, the symbol for cobalt, Co, consists of the symbols for carbon, C, and oxygen, O. A traditional NFA regular expression engine uses backtracking and could take exponentially long to test whether the word cococo…j matches a regular expression for concatenated chemical symbols. If there are _n_ copies of `co`, the engine would try all 2 _n_ possibilities of treating each pair as either `co` (cobalt) or `c o` (carbon oxygen). A DFA regular expression engine always runs in _O_(_n_) time, but a (traditional) NFA can take exponential time. So why ever use an NFA engine? Because it can match more advanced patterns such as Perl regular expressions (which are not “regular” expressions in the original meaning of the term). Newer NFA engines have added safeguards to eliminate or at least reduce the chance of running into catastrophic backtracking.

NIST published a set of elliptic curves for cryptography in 2000. The seeds for the curves were supplied by the NSA, and cryptographers suspected the agency may have installed backdoors in the curves. And not without reason: the NSA did insert a backdoor in the Dual_EC_DRBG random number generator. Two days ago researchers at Anthropic [1] showed that the NIST curves were not seeded with a backdoor. The rest of this post will give the background necessary to appreciate this result. * * * The security of elliptic curve cryptography depends on the difficult of solving the discrete logarithm problem. In the context of elliptic curve arithmetic, the analog of exponentiation is easy but the analog of taking logarithms is hard. **In general** it’s not practical to compute discrete logarithms over large elliptic curves. But it could be easy for someone with inside information to compute discrete logarithms if something about the curve was not “general.” Specifically, the NIST curves have the form _y_ ² = _x_ ³ − 3 _x_ + _b_ (mod _p_) for various values of _p_ [2], and the suspicion regards the choice of _b_. For nearly all values of _b_ encryption over the elliptic curve should be secure, as far as we know, but in principle someone could engineer values of _b_ to weaken the encryption. So where did the _b_ ‘s come from? Were they chosen in good faith? Concern over the choice of seeds, the _b_ ‘s, grew following the Snowden revelations. Daniel J. Bernstein emailed Jerry Solinas of the NSA on September 6, 2014 asking where the seeds had come from. Solinas replied in subsequent emails. He said the seeds were generated by hashing the ASCII representation of a humorous English message, but he could no longer remember the exact message or the precise hashing procedure. On October 5, 2023 Philippo Valsorda announced a bounty for reverse engineering the origin of the NIST elliptic curve seeds. Two days ago researchers from Anthropic claimed the bounty. Valsorda verified the results and awarded the prize. It turns out Solinas found the curve seeds in good faith, but his memory had faded over the years since the curves were published. We now know he used a date, not a phrase, as part of the process of constructing the seeds. The NIST curves have some nice properties, such as enabling some clever tricks for arithmetic in their base fields. But adoption of the curves has been held back over suspicion regarding their origin. Now the curves have been exhonorated. As Valdorda notes > … this nicely completes the “rehabilitation arc” of NIST elliptic curve cryptography, which I would now personally recommend over Curve25519/Ed25519/X25519. Just in time for them to be retired in favor of post-quantum signatures. ## Related posts * What is an elliptic curve? * Monero’s elliptic curve * Lewis Carroll and zero knowledge proofs [1] This does not imply that LLMs are necessarily a threat to encryption. The search for the NIST seeds required a lot of brute force. Some cleverness was required in order to come up with possibilities to check by brute force, but even that is sorta brute force at a higher level of abstraction. There has been speculation that labs are looking for weaknesses in cryptographic primitives and protocols, but that’s a different matter than searching for seeds. [2] The curves P-192, P-224, P-256, P-384, and P-521 are built over the field of integers mod _p_ , where the number of bits in _p_ corresponds to the number in the name. The P-192 curve, for example, has as its base field the integers mod 2192 − 264 − 1, a 192-bit prime.

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.