John D. Cook

Suppose you need to write software to compute the sine function. You were told in a calculus class that this is done using Taylor series-it's not, but that's another story-and so you start writing code to implement Taylor series. How many terms do you need? Uh, ..., 20? Let's go with that. from math import [...]The post Naively computing sine first appeared on John D. Cook.

The Euler-Mascheroni constant is defined as the limit So an obvious way to try to calculate would be to evaluate the right-hand side above for large n. This turns out to not be a very good approach. Convergence is slow and rounding error accumulates. A much better approach is to compute It's not obvious [...]The post Computing the Euler-Mascheroni Constant first appeared on John D. Cook.

It is possible to express every positive integer as a sum of powers of the golden ratio using each power at most once. This means it is possible to create a binary-like number system using as the base with coefficients of 0 and 1 in front of each power of . This system [...]The post Golden ratio base numbers first appeared on John D. Cook.

Sometimes a paper comes out that has the seeds of a great idea that could lead to a whole new line of pioneering research. But, instead, nothing much happens, except imitative works that do not push the core idea forward at all. For example the McCulloch Pitts paper from 1943 showed how neural networks could [...]The post Scientific papers: innovation ... or imitation? first appeared on John D. Cook.

I just stumbled across the binomial number system in Exercise 5.38 of Concrete Mathematics. The exercise asks the reader to show that every non-negative integer n can be written as and that the representation is unique if you require 0 a <b <c. The book calls this the binomial number system. I skimmed a paper [...]The post Binomial number system first appeared on John D. Cook.

Casting out nines is a well-known way of finding the remainder when a number is divided by 9. You add all the digits of a number n. And if that number is bigger than 9, add all the digits of that number. You keep this up until you get a number less than 9. This [...]The post Additive and multiplicative persistence first appeared on John D. Cook.

Logarithms give a way to describe large numbers compactly. But sometimes even the logarithm is a huge number, and it would be convenient to take the log again. And maybe again. For example, consider googolplex = 10googol where googol = 10100. (The name googol is older than Google," and in fact the former inspired the [...]The post Iterated logarithm first appeared on John D. Cook.

The inverse of an approximation is an approximation of the inverse." This statement is either obvious or really clever. Maybe both. Logs and exponents Let's start with an example. For moderately small x, That means that near 1 (i.e. near 10 raised to a small number, The inverse of a good approximation is a good [...]The post Approximation of Inverse, Inverse of Approximation first appeared on John D. Cook.

Fermat's primality test Fermat's little theoremsays that ifpis a prime andais not a multiple ofp, then ap-1= 1 (modp). The contrapositive of Fermat's little theorem says if ap-1 1 (modp) then eitherp is not prime ora is a multiple ofp. The contrapositive is used to test whether a number is prime. Pick a numbera less [...]The post False witnesses first appeared on John D. Cook.

When I hear American cities criticized for lack of long-term planning, my initial thought is that this is a good thing because the future cannot be predicted or directed very far in advance, and cities should be allowed to grow organically. Houston, in particular, is subject to a lot of criticism because of its lack [...]The post Houston's long term planning first appeared on John D. Cook.

There are several ways to define an ellipse. If you want to write software draw an ellipse, it's most convenient to have a parametric form: This gives an ellipse centered at (h, k), with semi-major axis a, semi-minor axis b, and with the major axis rotated by an angle from horizontal. But if you're [...]The post Gardener's ellipse first appeared on John D. Cook.

The previous post discussed fitting an ellipse and a parabola to the same data. Both fit well, but the ellipse fit a little better. This will often be the case because an ellipse has one more degree of freedom than a parabola. There is one way to fit a parabola to an ellipse at an [...]The post Fitting a parabola to an ellipse and vice versa first appeared on John D. Cook.

Thenth row of Pascal's triangle contains the binomial coefficientsC(n,r) forr ranging from 0 ton. For large n, if you print out the numbers in the nth row vertically in binary you can see a circular arc. Here's an example withn = 50. You don't have to use binary. Here's are the numbers in the row [...]The post The ellipse hidden inside Pascal's triangle first appeared on John D. Cook.

Imagine you are a soldier in charge of stacking cannonballs. Your fort has a new commander, a man with OCD who wants the cannonballs stacked in a very particular way. The new commander wants the balls stacked in tetrahedra. The balls on the bottom of the stack are arranged into a triangle. Then the next [...]The post Stacking positive and negative cannonballs first appeared on John D. Cook.

Mathematica includes code to draw various whimsical images. For example, if you enter the following command in Mathematica Entity["PopularCurve", "DragonflyCurve"][ EntityProperty["PopularCurve", "Image"]] you get an image of a dragonfly. It draws such images with Fourier series. You can tell by asking for the parameterization of the curve. If you enter Entity["PopularCurve", "DragonflyCurve"][ EntityProperty["PopularCurve", "ParametricEquations"]] you'll [...]The post How Mathematica Draws a Dragonfly first appeared on John D. Cook.

How often do the hour and and minute hand of an analog clock point in the same direction? The hands start pointing the same direction at midnight, then the hour hand moves clockwise (!) by 360 in 12 hours. The minute hand moves clockwise at the rate of 360 in one hour. So when does [...]The post How often do the hands of a clock line up? first appeared on John D. Cook.

Ever since 1952, the largest known primes have all had the form 2n - 1, with one exception in 1989. The reason the largest known primes have this form is that it is easier to test whether these numbers are prime than other numbers. A number of the form 2n - 1 is called a [...]The post The bad version of a good test first appeared on John D. Cook.

It's possible to represent complex numbers as a pair of real numbers or 2 * 2 matrices with real entries. And it's possible to represent quaternions as pairs of complex numbers or 2 * 2 matrices with complex entries werez* is the complex conjugate ofz. And it's also possible to represent octonions as pairs of [...]The post Representing octonions as matrices, sorta first appeared on John D. Cook.

You can multiply pairs of real numbers using the rules of complex numbers. Complex numbers have all the algebraic structure of the real numbers, i.e. they form a field. There is a general process, the Cayley-Dickson construction, that let's you bootstrap multiplication from 1 real number to 2, from 2 to 4, from 4 to [...]The post Octonions sometimes associate first appeared on John D. Cook.

There's an old joke about a drunk man looking for his keys under a lamppost. Someone stops and offers to help. He asks, So, did you lose your keys here?" The drunk replies No, I lost them over there, but here's where the light is." I routinely talk to people who have strong technical skills [...]The post Looking for keys under the lamppost first appeared on John D. Cook.

I wrote two unrelated blog posts this morning, one about the math paper H = W and one about a joke putting numbers into the D&D alignment matrix. I used Grok to look up the reference to the H = W paper, and to confirm that the alignment matrix originated with Dungeons & Dragons [1]. [...]The post Structured frameworks for complex systems first appeared on John D. Cook.

Dan Piponi posted a chart like the one below on Mastodon. At the risk of making a joke not funny by explaining it, I'd like to explain Dan's table. The alignment matrix above comes from Dungeons & Dragons and has become a kind of meme. The number neutral good number is the golden ratio, [...]The post Dungeons, Dragons, and Numbers first appeared on John D. Cook.

A paper came out in 1964 with the title H = W." The remarkably short title was not cryptic, however. The people for whom the paper was written knew exactly what it meant. There were two families of function spaces, one denoted with H and another denoted withW, that were speculated to be the same, [...]The post My favorite paper: H = W first appeared on John D. Cook.

If you change the coefficients of a polynomial a little bit, do you change the location of its zeros a little bit? In other words, do the roots of a polynomial depend continuously on its coefficients? You would think so, and you'd be right. Sorta. It's easy to see that small change to a coefficient [...]The post Wilkinson's polynomial first appeared on John D. Cook.

You would think that interpolating at more points would give you a more accurate approximation. There's a famous example by Runge that proves this is not the case. If you interpolate the function 1/(1 + x^2) over the interval [-5, 5], as you add more interpolation points the maximum interpolation error actually increases. Here's an [...]The post Interpolation instability first appeared on John D. Cook.

The most well-known generalization of the inverse of a matrix is the Moore-Penrose pseudoinverse.But there is another notion of generalized inverse, the Drazin pseudoinverse, for square matrices. If a matrix A has an inverseA-1 then it also has a Moore-Penrose pseudoinverse A+ and a Drazin pseudoinverse AD and A-1 = A+ = AD. When the [...]The post Drazin pseudoinverse first appeared on John D. Cook.

I've mentioned the Moore-Penrose pseudoinverse of a matrix a few times, most recently last week. This post will give an application of the pseudoinverse: computing effective graph resistance. Given a graphG, imagine replacing each edge with a one Ohm resistor. The effective resistance between two nodes inG is the electrical resistance between those the two [...]The post Effective graph resistance first appeared on John D. Cook.

An earlier post claimed that there practical advantages to partitioning a matrix, thinking of the matrix as a matrix of matrices. This post will give an example. LetM be a square matrix and suppose we need to multiply M by its transpose MT. We can compute this product faster than multiplying two arbitrary matrices of [...]The post Multiplying a matrix by its transpose first appeared on John D. Cook.

Pop quiz: what does the following code do? bool is_leap_year_fast(uint32_t y) { return ((y * 1073750999) & 3221352463) <= 126976; } It determines whether the year y is a leap year in the Gregorian calendar, of course. :) It's very efficient, though I don't image that would ever matter. But it's very clever. Gregorian vs [...]The post A bit-twiddling marvel first appeared on John D. Cook.

It's often useful to partition a matrix, thinking of it as a matrix whose entries are matrices. For example, you could look at the matrix 6 * 6 as a 2 * 2 matrix whose entries are 3 * 3 matrices. M is not a diagonal matrix of real numbers, but its block structure is [...]The post Matrices of Matrices first appeared on John D. Cook.

IJK (@iconjack) has calculated the probability of rolling a Yahtzee in no more than n rolls. The first few numerical values are: p(1) = 0.0007716 p(2) = 0.0126316 p(3) = 0.0460286 The probability is 0.95 after 23 rolls, and 0.99 after 32 rolls. Here's a plot.The post Probability of rolling a Yahtzee first appeared on John D. Cook.

The term trinomial coefficient is used a couple different ways. The most natural use of the term is a generalization of bionomial coefficients to three variables, i.e. numbers of the form wherei +j +k =n. These numbers are the coefficients of yj zk in the expansion of (x + y + z)n. But there's another [...]The post Trinomial Coefficients and Kings first appeared on John D. Cook.

I saw an integral posted online that came from this year's MIT Integration Bee. My thoughts on seeing this were, in order: It looks like a beta function. The answer is a small number. You can evaluate the integral using the substitution u = 1 - x2025. I imagine most students' reactions would be roughly [...]The post Riff on an integration bee integral first appeared on John D. Cook.

Here's another little chess puzzle by Martin Gardner, taken from the same paper as earlier. The task is to place two rooks, two bishops, and two knights on a 4 by 4 chessboard so that no piece attacks any other. As before, there are two basic solutions, here and here, plus symmetries.The post Another little chess puzzle first appeared on John D. Cook.

The map that takes a quaternionx to the quaternion qx is linear, so it can be represented as multiplication by a matrix. The same is true of the map that takesx toxq, but the two matrices are not the same because quaternion multiplication does not commute. Letq =a + bi + cj +dkand let qM [...]The post Multiplying by quaternions on the left and right first appeared on John D. Cook.

I just revised a post from a week ago about rotations. The revision makes explicit the process of embedding a 3D vector into the quaternions, then pulling it back out. The 3D vector is embedded in the quaternions by making it the vector part of a quaternion with zero real part: (p1,p2,p3) (0, p1,p2,p3) [...]The post Embeddings, Projections, and Inverses first appeared on John D. Cook.

The other day I stumbled on an article [1] that advocated writing abasab and loga(b) as ab. This is a special case of Knuth's up arrow and down arrow notation. Knuth introduces his arrows with the intention of repeating them to represent hyperexponentiation and iterated logarithms. But the emphasis in [1] is more on the [...]The post Alternative exp and log notation first appeared on John D. Cook.

This morning I ran across the following tip from Joost Helberg on Mastodon: TIL don't report numbers with three digits after the decimal point. People may interpret the decimal point as a thousands separator. Using 2 or 4 digits, although wrong, avoids off by a factor thousand errors. I usually report four decimal places, but [...]The post Decimal Separator and Internationalization first appeared on John D. Cook.

Here's a puzzle by Martin Garnder [1]. Can a queen, king, rook, bishop, and knight be placed on a 4^2 board so no piece attacks another? There are two solutions, plus symmetries. Note that in all non-attacking chess puzzles, the colors of the pieces are irrelevant. In the solutions I chose the piece colors to [...]The post A crowded little chess puzzle first appeared on John D. Cook.

The Boolean satisfiability problem is to determine whether there is a way to assign values to variables in a set of Boolean formulas to make the formulas hold [1]. If there is a solution, the next task would be to enumerate the solutions. You can solve the famous eight queens problem, or its generalization ton-queens, [...]The post Formulating eight queens as a SAT problem first appeared on John D. Cook.

There are 92 ways to place eight queens on a chessboard so that no queen is attacking any other. These fall into 12 equivalence classes. The 92 solutions are all rotations and reflections of these 12 basic solutions. If you think about the previous numbers a minute, you might wonder why the total number of [...]The post Special solutions to the eight queens problem first appeared on John D. Cook.

How many bishops can you place on a chessboard so that no bishop is attacking any other bishop? For a standard 8 * 8 chessboard the answer is 14. In general, for an n * n chessboard the answer is 2n - 2. Here's one way to place the maximum number of non-attacking bishops. To [...]The post The non-attacking bishops problem first appeared on John D. Cook.

This morning I wrote about the frequencies of names for popes and kings. This involved sorting strings with Roman numerals since it's common for popes and kings to have Roman numerals after their names. Something that surprised me was that sorting Roman numerals alphabetically roughly sorts them in numerical order, especially for small numbers. It's [...]The post Sorting Roman numerals first appeared on John D. Cook.

After I wrote the code to make the bar graph of papal names for the previous post, I decided to reuse the code to make a similar graph for monarchs of England. Just as there is some complication in counting papal names, there are even more complications in counting names of English monarchs. Who was [...]The post Frequency of names of English monarchs first appeared on John D. Cook.

The new pope chose the name Leo XIV. That made me curious about the distribution of names of popes and so I made the graph below. (I'm Protestant, so wasn't familiar to me.) Looks like Leo is tied with Clement for fourth place, the top three names being John, Benedict, and Gregory. There are a [...]The post Frequency of papal names first appeared on John D. Cook.

The recent post Converting between quaternions and rotation matrices describes an algorithm as a mathematically correct but numerically suboptimal method/" Let's just look at a small piece of that post. The post explains how to find a rotation matrix to describe the same rotation as pre- and post- multiplying by a unit quaternion q, and [...]The post Square root of a small number first appeared on John D. Cook.

Large language models display emergence behaviors: when the parameter count is scaled to a certain value, suddenly the LLM is capable of performing a new task not possible at a smaller size.Some say the abruptness of this change is merely a spurious artifact of how it is measured. Even so, many would like to understand, [...]The post Why do LLMs have emergent properties? first appeared on John D. Cook.

In the previous post I wrote about representing rotations with quaternions. This representation has several advantages, such as making it clear how rotations compose. Rotations are often represented as matrices, and so it's useful to be able to go between the two representations. A unit-length quaternion (q0, q1, q2, q3) represents a rotation by an [...]The post Converting between quaternions and rotation matrices first appeared on John D. Cook.

Every rotation in 3D fixes an axis [1]. This is Euler's rotation theorem from 1775. Another way to state the theorem is that no matter how you rotate a sphere about its center, two points stay fixed. The composition of two rotations is a rotation. So the first rotation fixes some axis, the second rotation [...]The post Composing rotations using quaternions first appeared on John D. Cook.

This is the 5,000th post on this blog. I started blogging in 2008, and Wayne Joubert started contributing posts last year. We've written an average of between five and six posts a week for the last 17 years. I thought about listing some of the most popular posts over the years, but that's not as [...]The post 5,000th post first appeared on John D. Cook.