John D. Cook

Our team started working within Nvidia in early 2009 at the beginning of the ORNL Titan project. Our Nvidia contacts dealt with applications, libraries, programming environment and performance optimization. First impressions were that their technical stance on issues was very reasonable. One obscure example: in a C++ CUDA kernel were you allowed to use enums," [...]The post Experiences with Nvidia first appeared on John D. Cook.

The previous post mentioned Legendre polynomials. This post will give a brief introduction to these polynomials and a couple hints of how they are used in applications. One way to define the Legendre polynomials is as follows. P0(x) = 1 Pk are orthogonal on [-1, 1]. Pk(1) = 1 for all k >= 0. The [...]The post Legendre polynomials first appeared on John D. Cook.

To first approximation, a satellite orbiting the earth moves in an elliptical orbit. That's what would get from solving the two-body problem: two point masses orbiting their common center of mass, subject to no forces other than their gravitational attraction to each other. But the earth is not a point mass. Neither is a satellite, [...]The post The biggest perturbation of satellite orbits first appeared on John D. Cook.

Suppose you have a radio transmitter T and a receiver R with a clear line of sight between them. Some portion the signal received atR will come straight fromT. But some portion will have bounced off some obstacle, such as the ground. The reflected radio waves will take a longer path than the waves that [...]The post Transmission Obstacles and Ellipsoids first appeared on John D. Cook.

The exponential sum of the day page on my site draws an image every day by plugging the month, day, and year into a formula. Details here. Today's image looks almost solid blue in the middle. The default plotting line width works well for most days. For example, see what the sum of the day [...]The post Zooming in on a fractalish plot first appeared on John D. Cook.

All integers are real numbers, but most computer representations of integers do not equal computer representations of real numbers. To make the statement above precise, we have to be more specific about what we mean by computer integers and floating point numbers. I'll use int32 and int64 to refer to 32-bit and 64-bit signed integers. [...]The post Most ints are not floats first appeared on John D. Cook.

Years ago I wrote about a fast way to test whether an integer n is a square. The algorithm rules out a few numbers that cannot be squares based on their last (hexadecimal) digit. If the the integer passes through this initial screening, the algorithm takes the square root of n as a floating point [...]The post Test whether a large integer is a square first appeared on John D. Cook.

What does an eighteenth century French mathematician have to do with the Ethereum cryptocurrency? A pseudorandom number generator based on Legendre symbols, known as Legendre PRF, has been proposed as part of a zero knowledge proof mechanism to demonstrate proof of custody in Ethereum 2.0. I'll explain what each of these terms mean and include [...]The post Legendre and Ethereum first appeared on John D. Cook.

The previous post looked at a function formed by patching together the function f(x) = log(1 +x) for positivex and f(x) =x for negativex. The functions have a mediocre fit, which may be adequate for some applications, such as plotting data points, but not for others, such as visual design. Here's a plot. There's something [...]The post Patching functions together first appeared on John D. Cook.

I saw a post online this morning that recommended the transformation I could see how this could be very handy. Often you want something like a logarithmic scale, not for the exact properties of the logarithm but because it brings big numbers closer in. And for big values ofx there's little difference between log(x) and [...]The post Log-ish first appeared on John D. Cook.

Suppose you have a system withn possible states. The entropy of the system is maximized when all states are equally likely to occur. The entropy is minimized when one outcome is certain to occur. You can say more. Starting from any set of probabilities, as you move in the direction of more uniformity, you increase [...]The post Uniformity increases entropy first appeared on John D. Cook.

The sinc function sinc(x) = sin(x) /x comes up continually in signal processing. Ifx is moderately small, the approximation sinc(x) (2 + cos(x))/3 is remarkably good, with an error on the order of x4/180. This could be useful in situations where you're working with the sinc function and thex in the denominator is awkward [...]The post Sinc function approximation first appeared on John D. Cook.

Four years ago I wrote a blog post about simple solutions to client problems. The post opens by recounting a conversation with a friend that ended with my friend saying So, basically you're recommending division." That conversation came back to me recently when I had a similar conversation during a deposition. I had to explain [...]The post Arithmetic for fun and profit first appeared on John D. Cook.

A couple days ago I wrote about the the problem that Bitcoin requires to be solved as proof-of-work. In a nutshell, you need to tweak a block of transactions until the SHA256 double hash of its header is below a target value [1]. Not all cryptocurrencies use proof of work, but those that do mostly [...]The post Why use hash puzzles for proof-of-work? first appeared on John D. Cook.

Suppose you have a list of positive data points y1, y2, ..., yn and you wanted to find a value that minimizes the squared distances to each of the ys. Then the solution is to take to be the mean of the ys: This result is well known [1]. The following variation is [...]The post Minimize squared relative error first appeared on John D. Cook.

In order to prevent fraud, anyone wanting to add a block to the Bitcoin blockchain must prove that they've put in a certain amount of computational work. This post will focus on what problem must be solved in order produce proof of work. You'll see the proof of work function described as finding strings whose [...]The post What is the Bitcoin proof-of-work problem? first appeared on John D. Cook.

This post looks at whether you should delete names or replace names when deidentifying personal data. With structured data, generating synthetic names does not increase or decrease privacy. But with unstructured data, replacing real names with randomly generated names increases privacy protection. Structured data If you want to deidentify structured data (i.e. data separated into [...]The post Deleting vs Replacing Names first appeared on John D. Cook.

I went down a rabbit hole this week using two symbols in LaTeX. The first was the Russian letter Sha (, U+0248), and the second was the currency symbol for Bitcoin (, U+20BF). Sha I thought there would be a LaTeX package that would include as a symbol rather than as a Russian letter, [...]The post Typesetting Sha and Bitcoin first appeared on John D. Cook.

Years ago I wrote a post Golden powers are nearly integers. The post was picked up by Hacker News and got a lot of traffic. The post was commenting on a line from Terry Tao: The powers , 2, 3, ... of the golden ratio lie unexpectedly close to integers: for instance, 11= 199.005... is [...]The post Golden powers revisited first appeared on 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.