Solving hard problems

We help companies solve hard problems in mathematics, statistics, and computing. Let’s explore how we might work together.

Integer odds and prime numbers

Posted on 17 October 2018 by John

For every integer m > 1, it’s possible to choose N so that the proportion of primes in the sequence 1, 2, 3, … N is 1/m. To put it another way, you can make the odds against one of the first N natural numbers being prime any integer value you’d like [1].

For example, suppose you wanted to find N so that 1/7 of the first N positive integers are prime. Then the following Python code shows you could pick N = 3059.

    from sympy import primepi

    m = 7

    N = 2*m
    while N / primepi(N) != m:
        N += m
    print(N)

[1] Solomon Golomb. On the Ratio of N to π(N). The American Mathematical Monthly, Vol. 69, No. 1 (Jan., 1962), pp. 36-37.

The proof is short, and doesn’t particularly depend on the distribution of primes. Golomb proves a more general theorem for any class of integers whose density goes to zero.

My most popular posts on Reddit

Posted on 16 October 2018 by John

There are only three posts on this top 10 list that are also on the top 10 list for Hacker News.

Comparing trig functions and Jacobi functions

Posted on 13 October 2018 by John

My previous post looked at Jacobi functions from a reference perspective: given a Jacobi function defined one way, how do I relate that to the same function defined another way, and how would you compute it?

This post explores the analogy between trigonometric functions and Jacobi elliptic functions.

Related basic Jacobi functions to trig functions

In the previous post we mentioned a connection between the argument u of a Jacobi function and the amplitude φ:

$u = \int_0^{\varphi} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$

We can use this to define the functions sn and cn. Leaving the dependence on m implicit, we have

$\begin{align*} \mathrm{sn}(u) &= \sin(\varphi) \\ \mathrm{cn}(u) &= \cos(\varphi) \end{align*}$

If m = 0, then u = φ and so sn and cn are exactly sine and cosine.

There’s a third Jacobi function we didn’t discuss much last time, and that’s dn. It would be understandable to expect dn might be analogous to tangent, but it’s not. The function dn is the derivative of φ with respect to u, or equivalently

$\mathrm{dn}(u) = \sqrt{1 - m \sin^2\varphi}$

The rest of the Jacobi functions

Just as there are more trig functions than just sine and cosine, there are more Jacobi functions than sn, cn, and dn. There are two ways to define the rest of the Jacobi functions: in terms of ratios, and in terms of zeros and poles.

Ratios

I wrote a blog post one time asking how many trig functions there are. The answer depends on your perspective, and I gave arguments for 1, 3, 6, and 12. For this post, lets say there are six: sin, cos, tan, sec, csc, and ctn.

One way to look at this is to say there are as many trig functions as there are ways to select distinct numerators and denominators from the set {sin, cos, 1}. So we have tan = sin/cos, csc = 1/sin, etc.

There are 12 Jacobi elliptic functions, one for each choice of distinct numerator and denominator from {sn, cn, dn, 1}. The name of a Jacobi function is two letters, denoting the numerator and denominator, where we have {s, c, d, n} abbreviating {sn, cn, dn, 1}.

For example, cs(u) = sn(u) / cn(u) and ns(u) = 1 / sn(u).

Note that to take the reciprocal of a Jacobi function, you just reverse the two letters in its name.

Zeros and poles

The twelve Jacobi functions can be classified [1] in terms of their zeros and poles over a rectangle whose sides have length equal to quarter periods. Let’s look at an analogous situation for trig functions before exploring Jacobi functions further.

Trig functions are periodic in one direction, while elliptic functions are periodic in two directions in the complex plane. A quarter period for the basic trig functions is π/2. The six trig functions take one value out of {0, 1, ∞} at 0 and different value at π/2. So we have one trig function for each of the six ways to chose an permutation of length 2 from a set of 3 elements.

In the previous post we defined the two quarter periods K and K‘. Imagine a rectangle whose corners are labeled

s = (0, 0)
c = (K, 0)
d = (K, K‘)
n = (0, K‘)

Each Jacobi function has a zero at one corner and a pole at another. The 12 Jacobi function correspond to the 12 ways to chose a permutation of two items from a set of four.

The name of a Jacobi function is two letters, the first letter corresponding to where the zero is, and the second letter corresponding to the pole. So, for example, sn has a zero at s and a pole at n. These give same names as the ratio convention above.

Identities

The Jacobi functions satisfy many identities analogous to trigonometric identities. For example, sn and cn satisfy a Pythagorean identity just like sine and cosine.

$\mathrm{sn}^2 u + \mathrm{cn}^2 u = 1$

Also, the Jacobi functions have addition theorems, though they’re more complicated than their trigonometric counterparts.

$\begin{align*} \mathrm{sn}(u + v) &= \frac{\mathrm{sn}\,u\, \mathrm{cn}\, v\, \mathrm{dn}\,v\, + \mathrm{sn}\,v\, \mathrm{cn}\, u\, \mathrm{dn}\,u\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \\ \\ \mathrm{\mathrm{cn}\,}(u + v) &= \frac{\mathrm{cn}\, u\, \mathrm{cn}\, v\, - \mathrm{sn}\,u\, \mathrm{dn}\,u\, \mathrm{sn}\,v\, \mathrm{dn}\,v\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \\ \\ \mathrm{dn}(u + v) &= \frac{\mathrm{dn}\,u\, \mathrm{dn}\,v\, - m\, \mathrm{sn}\,u\, \mathrm{cn}\, u\, \mathrm{sn}\,v\, \mathrm{cn}\, v\,}{1 - m\, \mathrm{sn}^2 u\, \, \mathrm{sn^2} v\,} \end{align*}$

Derivatives

The derivatives of the basic Jacobi functions are given below.

$\begin{align*} \mathrm{sn}'(u) &= \mathrm{cn}(u)\, \mathrm{dn}(u) \\ \\ \mathrm{cn}'(u) &= -\mathrm{sn}(u) \,\mathrm{dn}(u) \\ \\ \mathrm{dn}'(u) &= -m\,\mathrm{sn}(u)\, \mathrm{cn}(u) \\ \end{align*}$

Note that the implicit parameter m makes an appearance in the derivative of dn. We will also need the complementary parameter m‘ = 1 – m.

The derivatives of all Jacobi functions are summarized in the table below.

$\begin{table} \centering \begin{tabular}{l|rrrr} \multicolumn{1}{l}{} & s & n & d & c \\ \cline{2-5} s & & dn cn & nd cd & nc dc \\ n & $-$ ds cs & & $m$ sd cd & sc dc \\ d & $-$ ns cs & $-m$ sn cn & & $m'$ sc nc \\ c & $-$ ns ds & $-$ sn dn & $-m'$ sd nd & \end{tabular} \end{table}$

The derivatives of the basic Jacobi functions resemble those of trig functions. They may look more complicated at first, but they’re actually more regular. You could remember them all by observing the patterns pointed out below.

Let w, x, y, z be any permutation of {s, n, d, c}. Then the derivative of wx is proportional to yx zx. That is, the derivative of every Jacobi function f is a constant times two other Jacobi functions. The names of these two functions both end in the same letter as the name of f, and the initial letters are the two letters not in the name of f.

The proportionality constants also follow a pattern. The sign is positive if and only if the letters in the name of f appear in order in the sequence s, n, d, c. Here’s a table of just the constants.

$\begin{table} \centering \begin{tabular}{l|rrrr} \multicolumn{1}{l}{} & s & n & d & c \\ \cline{2-5} s & & 1 & 1 & 1 \\ n & $-1$ & & $m$ & 1 \\ d & $-1$ & $-m$ & & $m'$ \\ c & $-1$ & $-1$ & $-m'$ & \end{tabular} \end{table}$

Note that the table is skew symmetric, i.e. its transpose is its negative.

[1] An elliptic function is determined, up to a constant multiple, by its periods, zeros, and poles. So not only do the Jacobi functions have the pattern of zeros and poles described here, these patterns uniquely determine the Jacob functions, modulo a constant. For (singly) periodic functions, the period, zeros, and poles do not uniquely determine the function. So the discussion of zeros and poles of trig functions is included for comparison, but it does not define the trig functions.

Clearing up the confusion around Jacobi functions

Posted on 12 October 2018 by John

The Jacobi elliptic functions sn and cn are analogous to the trigonometric functions sine and cosine. The come up in applications such as nonlinear oscillations and conformal mapping. Unfortunately there are multiple conventions for defining these functions. The purpose of this post is to clear up the confusion around these different conventions.

Plot of Jacobi sn

The image above is a plot of the function sn [1].

Modulus, parameter, and modular angle

Jacobi functions take two inputs. We typically think of a Jacobi function as being a function of its first input, the second input being fixed. This second input is a “dial” you can turn that changes their behavior.

There are several ways to specify this dial. I started with the word “dial” rather than “parameter” because in this context parameter takes on a technical meaning, one way of describing the dial. In addition to the “parameter,” you could describe a Jacobi function in terms of its modulus or modular angle. This post will be a Rosetta stone of sorts, showing how each of these ways of describing a Jacobi elliptic function are related.

The parameter m is a real number in [0, 1]. The complementary parameter is m‘ = 1 – m. Abramowitz and Stegun, for example, write the Jacobi functions sn and cn as sn(u | m) and cn(u | m). They also use m₁ = rather than m‘ to denote the complementary parameter.

The modulus k is the square root of m. It would be easier to remember if m stood for modulus, but that’s not conventional. Instead, m is for parameter and k is for modulus. The complementary modulus k‘ is the square root of the complementary parameter.

The modular angle α is defined by m = sin² α.

Note that as the parameter m goes to zero, so does the modulus k and the modular angle α. As any one of these three goes to zero, the Jacobi functions sn and cn converge to their counterparts sine and cosine. So whether your dial is the parameter, modulus, or modular angle, sn converges to sine and cn converges to cosine as you turn the dial toward zero.

As the parameter m goes to 1, so does the modulus k, but the modular angle α goes to π/2. So if your dial is the parameter or the modulus, it goes to 1. But if you think of your dial as modular angle, it goes to π/2. In either case, as you turn the dial to the right as far as it will go, sn converges to hyperbolic secant, and cn converges to the constant function 1.

Quarter periods

In addition to parameter, modulus, and modular angle, you’ll also see Jacobi function described in terms of K and K‘. These are called the quarter periods for good reason. The functions sn and cn have period 4K as you move along the real axis, or move horizontally anywhere in the complex plane. They also have period 4iK‘. That is, the functions repeat when you move a distance 4K‘ vertically [2].

The quarter periods are a function of the modulus. The quarter period K along the real axis is

$K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$

and the quarter period K‘ along the imaginary axis is given by K‘(m) = K(m‘) = K(1 – m).

The function K(m) is known as “the complete elliptic integral of the first kind.”

Amplitude

So far we’ve focused on the second input to the Jacobi functions, and three conventions for specifying it.

There are two conventions for specifying the first argument, written either as φ or as u. These are related by

$u = \int_0^{\varphi} \frac{d\theta}{\sqrt{1-m\sin^2\theta}}$

The angle φ is called the amplitude. (Yes, it’s an angle, but it’s called an amplitude.)

When we said above that the Jacobi functions had period 4K, this was in terms of the variable u. Note that when φ = π/2, u = K.

Jacobi elliptic functions in Mathematica

Mathematica uses the u convention for the first argument and the parameter convention for the second argument.

The Mathematica function JacobiSN[u, m] computes the function sn with argument u and parameter m. In the notation of A&S, sn(u | m).

Similarly, JacobiCN[u, m] computes the function cn with argument u and parameter m. In the notation of A&S, cn(u | m).

We haven’t talked about the Jacobi function dn up to this point, but it is implemented in Mathematica as JacobiDN[u, m].

The function that solves for the amplitude φ as a function of u is JacobiAmplitude[um m].

The function that computes the quarter period K from the parameter m is EllipticK[m].

Jacobi elliptic functions in Python

The SciPy library has one Python function that computes four mathematical functions at once. The function scipy.special.ellipj takes two arguments, u and m, just like Mathematica, and returns sn(u | m), cn(u | m), dn(u | m), and the amplitude φ(u, m).

The function K(m) is implemented in Python as scipy.special.ellipk.

[1] The plot was made using JacobiSN[0.5, z] and the function ComplexPlot described here.

[2] Strictly speaking, 4iK‘ is a period. It’s the smallest vertical period for cn, but 2iK‘ is the smallest vertical period for sn.

Hadamard product

Posted on 10 October 2018 by John

The first time you see matrices, if someone asked you how you multiply two matrices together, your first idea might be to multiply every element of the first matrix by the element in the same position of the corresponding matrix, analogous to the way you add matrices.

But that’s not usually how we multiply matrices. That notion of multiplication hardly involves the matrix structure; it treats the matrix as an ordered container of numbers, but not as a way of representing a linear transformation. Once you have a little experience with linear algebra, the customary way of multiplying matrices seems natural, and the way that may have seemed natural at first glance seems kinda strange.

The componentwise product of matrices is called the Hadamard product or sometimes the Schur product. Given two m by n matrices A and B, the Hadamard product of A and B, written A ∘ B, is the m by n matrix C with elements given by

c_ij = a_ij b_ij.

Because the Hadamard product hardly uses the linear structure of a matrix, you wouldn’t expect it to interact nicely with operations that depend critically on the linear structure. And yet we can give a couple theorems that do show a nice interaction, at least when A and B are positive semi-definite matrices.

The first is the Schur product theorem. It says that if A and B are positive semi-definite n by n matrices, then

det(A ∘ B) ≥ det(A) det(B)

where det stands for determinant.

Also, there is the following theorem of Pólya and Szegö. Assume A and B are symmetric positive semi-definite n by n matrices. If the eigenvalues of A and B, listed in increasing order, are α_i and β_i respectively, then for every eigenvalue λ of A ∘ B, we have

α₁ β₁ ≤ λ ≤ α_n β_n.

Python implementation

If you multiply two (multidimensional) arrays in NumPy, you’ll get the componentwise product. So if you multiply two matrices as arrays you’ll get the Hadamard product, but if you multiply them as matrices you’ll get the usual matrix product. We’ll illustrate that below. Note that the function eigvalsh returns the eigenvalues of a matrix. The name may look a little strange, but the “h” on the end stands for “Hermitian.” We’re telling NumPy that the matrix is Hermitian so it can run software specialized for that case [1].

    
    from numpy import array, matrix, array_equal, all
    from numpy.linalg import det, eigvalsh
    
    A = array([
        [3, 1],
        [1, 3]
    ])
    
    B = array([
        [5, -1],
        [-1, 5]
    ])
    
    H = array([
        [15, -1],
        [-1, 15]
    ])
    
    AB = array([
        [14,  2],
        [ 2, 14]
    ])
    
    # Hadamard product
    assert(array_equal(A*B, H))
    
    # Ordinary matrix product
    assert(array_equal(A@B, AB))
    
    # Schur product theorem
    assert(det(H) >= det(A)*det(B))
    
    # Eigenvalues
    eigA = eigvalsh(A)
    eigB = eigvalsh(B)
    eigH = eigvalsh(A*B)
    
    lower = eigA[0]*eigB[0]
    upper = eigA[1]*eigB[1]
    assert(all(eigH >= lower))
    assert(all(eigH <= upper))

The code above shows that the eigenvalues of A are [2, 4], the eigenvalues of B are [4, 6], and the eigenvalues of A ∘ B are [14, 16].

[1] For complex matrices, Hermitian means conjugate symmetric, which in the real case reduces to simply symmetric. The theorem of Pólya and Szegö is actually valid for Hermitian matrices, but I simplified the statement for the case of real-valued matrices.

Prime interruption

Posted on 8 October 2018 by John

Suppose you have a number that you believe to be prime. You start reading your number aloud, and someone interrupts “Stop right there! No prime starts with the digits you’ve read so far.”

It turns out the person interrupting you shouldn’t be so sure. There are no restrictions on the digits a prime number can begin with. (Ending digits are another matter. No prime ends in 0, for example.) Said another way, given any sequence of digits, it’s possible to add more digits to the end and make a prime. [1]

For example, take today’s date in ISO format: 20181008. Obviously not a prime. Can we find digits to add to make it into a prime? Yes, we can add 03 on to the end because 2018100803 is prime.

What about my work phone number: 83242286846? Yes, just add a 9 on the end because 832422868469 is prime.

This works in any base you’d like. For example, the hexadecimal number CAFEBABE is not prime, but CAFEBABE1 is. Or if you interpret SQUEAMISH as a base 36 number, you can form a base 36 prime by sticking a T on the end. [2]

In each of these example, we haven’t had to add much on the end to form a prime. You can show that this is to be expected from the distribution of prime numbers.

[1] Source: R. S. Bird. Integers with Given Initial Digits. The American Mathematical Monthly, Vol. 79, No. 4 (Apr., 1972), pp. 367-370

[2] CAFEBABE is a magic number used at the beginning of Java bytecode files. The word “squeamish” here is an homage to “The Magic Words are Squeamish Ossifrage,” an example cleartext used by the inventors of RSA encryption.

Probability of winning the World Series

Posted on 6 October 2018 by John

Astros win 2018 World Series

Suppose you have two baseball teams, A and B, playing in the World Series. If you like, say A stands for Houston Astros and B for Milwaukee Brewers. Suppose that in each game the probability that A wins is p, and the probability of A losing is q = 1 – p. What is the probability that A will win the series?

The World Series is a best-of-seven series, so the first team to win 4 games wins the series. Once one team wins four games there’s no point in playing the rest of the games because the series winner has been determined.

At least four games will be played, so if you win the series, you win on the 4th, 5th, 6th, or 7th game.

The probability of A winning the series after the fourth game is simply p⁴.

The probability of A winning after the fifth game is 4 p⁴ q because A must have lost one game, and it could be any one of the first four games.

The probability of A winning after the sixth game is 10 p⁴ q² because A must have lost two of the first five games, and there are 10 ways to choose two items from a set of five.

Finally, the probability of A winning after the seventh game is 20 p⁴ q³ because A must have lost three of the first six games, and there are 20 ways to choose three items from a set of six.

The probability of winning the World Series is the sum of the probabilities of winning after 4, 5, 6, and 7 games which is

p⁴(1 + 4q + 10q² + 20q³)

Here’s a plot:

Probability of winning the world series as a function of the probability of winning one game

Obviously, the more likely you are to win each game, the more likely you are to win the series. But it’s not a straight line because the better team is more likely to win the series than to win any given game.

Now if you only wanted to compute the probability of winning the series, not the probability of winning after different numbers of games, you could pretend that all the games are played, even though some may be unnecessary to determine the winner. Then we compute the probability that a Binomial(7, p) random variable takes on a value greater than or equal to 4, which is

35p⁴q³ + 21p⁵q² + 7p⁶q + p⁷

While looks very different than the expression we worked out above, they’re actually the same. If you stick in (1 – p) for q and work everything out, you’ll see they’re the same.

Physical constants in Python

Posted on 4 October 2018 by John

You can find a large collection of physical constants in scipy.constants. The most frequently used constants are available directly, and hundreds more are in a dictionary physical_constants.

The fine structure constant α is defined as a function of other physical constants:

$\alpha = \frac{e^2}{4 \pi \varepsilon_0 \hbar c}$

The following code shows that the fine structure constant and the other constants that go into it are available in scipy.constants.

    import scipy.constants as sc

    a = sc.elementary_charge**2
    b = 4 * sc.pi * sc.epsilon_0 * sc.hbar * sc.c
    assert( abs(a/b - sc.fine_structure) < 1e-12 )

Eddington’s constant

In the 1930’s Arthur Eddington believed that the number of protons in the observable universe was exactly the Eddington number

$N_{\mathrm{Edd}} = \frac{2^{256}}{\alpha}$

Since at the time the fine structure constant was thought to be 1/136, this made the number of protons a nice even 136 × 2²⁵⁶. Later he revised his number when it looked like the fine structure constant was 1/137. According to the Python code above, the current estimate is more like 1/137.036.

Eddington was a very accomplished scientist, though he had some ideas that seem odd today. His number is a not a bad estimate, though nobody believes it could be exact.

The constants in scipy.constants have come up in a couple previous blog posts.

The post on Koide’s coincidence shows how to use the physical_constants dictionary, which includes not just the physical constant values but also their units and uncertainty.

The post on Benford’s law shows that the leading digits of the constants in scipy.constants follow the logarithmic distribution observed by Frank Benford (and earlier by Simon Newcomb).

Technological context

Posted on 3 October 2018 by John

As you read this blog post, you probably have a hierarchy of contexts in the back of your mind. It comes so naturally to you that you don’t consciously think of it.

If you’re reading this in a web browser, you probably know what browser you’re using. If not, you’re at least aware that you are using a browser, even if you forget momentarily which one you have open. And you probably know what operating system is hosting your browser. You understand that you are reading a page on my site, that this page is not your browser, but content hosted by your browser. If you’ve subscribed via RSS or email, you know what email or RSS client you’re using and understand how this post is organized with respect to your other content.

Some people do not have this kind of context. Anything on their computer is simply “The Computer.” They don’t really understand what an operating system, web browser, or email client are. And they don’t need to know, most of the time. They can get their work done, but then occasionally they have inexplicable problems.

I’m not saying this to criticize or make fun of the people who don’t have the kind of technological context I’m talking about. It’s a remarkable achievement that software has gotten so easy to use that people can get along without knowing much technological context. But if this sort of thing is second nature to you, you might have a hard time understanding how a large number of people work.

You probably take it for granted that you can access the same web site from different computers. Some people do not. Their desktop at work is one computer, and their iPhone is a different computer. They don’t really understand what a web site is.

I know what a web browser is because I have been using computers since before there was a web. Old timers know what various technological components are because they’ve seen them develop. And “digital natives” know to get things done because they’ve grown up with computers, though their gaps in context show occasionally. Seems like the people in the middle would have the hardest time, not having grown up with the technology but not having watched it develop either.

I’m writing this because I’m becoming increasingly aware of how difficult life can be for people who don’t have adequate mental models for technology. I imagine most of my readers are tech savvy, and may have a hard time seeing some of the same things that I’ve had a hard time seeing, that a lot of people don’t understand things we take for granted.

Source http://dilbert.com/strip/1995-06-24

It used to be that anybody who used a computer had to know some basic things. If you were a Unix user a generation ago, you might not know anything about the internals of Unix, but you at least knew that you were a Unix user. There were wizards and mere mortals, but the two groups shared more context than the most tech savvy and least tech savvy share today.

It’s good that people don’t need to know as much context, but occasionally it produces bewildering situations, both for the user and the person trying to help them.

Interlaced zeros of ODEs

Posted on 2 October 2018 by John

Sturm’s separation theorem says that the zeros of independent solutions to an equation of the form

alternate. That is, between any two consecutive zeros of one solution, there is exactly one zero of the other solution. This is an important theorem because a lot of differential equations of this form come up in applications.

If we let p(x) = 0 and q(x) = 1, then sin(x) and cos(x) are independent solutions and we know that their zeros interlace. The zeros of sin(x) are of the form nπ, and the zeros of cos(x) are multiples of (n + 1/2)π.

What’s less obvious is that if we take two different linear combinations of sine and cosine, as long as they’re not proportional, then their zeros interlace as well. For example, we could take f(x) = 3 sin(x) + 5 cos(x) and g(x) = 7 sin(x) – 11 cos(x). These are also linearly independent solutions to the same differential equation, and so the Sturm separation theorem says their roots have to interlace.

If we take p(x) = 1/x and q(x) = 1 – (ν/x)² then our differential equation becomes Bessel’s equation, and the Bessel functions J_ν and Y_ν are independent solutions. Here’s a little Python code to show how the zeros alternate when ν = 3.

    import matplotlib.pyplot as plt
    from scipy import linspace
    from scipy.special import jn, yn

    x = linspace(4, 30, 100)
    plt.plot(x, jn(3, x), "-")
    plt.plot(x, yn(3, x), "-.")
    plt.legend(["$J_3$", "$Y_3$"])
    plt.axhline(y=0,linewidth=1, color="k")
    plt.show()

Plotting Bessel functions J_3 and Y_3

Thirty days hath September …

Posted on 1 October 2018 by John

Here’s today’s exponential sum, 10/1/2018:

And here’s what it would look like if September had 31 days, 9/31/2018:

The valley of medium reliability

Posted on 27 September 2018 by John

Last evening my electricity went out and this morning it was restored. This got met thinking about systems that fail occasionally [1]. Electricity goes out often enough that we prepare for it. We have candles and flashlights, my work computer is on a UPS, etc.

A residential power outage is usually just an inconvenience, especially if the power comes back on within a few hours. A power outage to a hospital could be disastrous, and so hospitals have redundant power systems. The problem is in between, if power is reliable enough that you don’t expect it to go out, but the consequences of an outage are serious [2].

If a system fails occasionally, you prepare for that. And if it never fails, that’s great. In between is the problem, a system just reliable enough to lull you into complacency.

Dangerously reliable systems

For example, GPS used to be unreliable. It made useful suggestions, but you wouldn’t blindly trust it. Then it got a little better and became dangerous as people trusted it when they shouldn’t. Now it’s much better. Not perfect, but less dangerous.

For another example, people who live in flood planes have flood insurance. Their mortgage company requires it. And people who live on top of mountains don’t need flood insurance. The people at most risk are in the middle. They live in an area that could flood, but since it hasn’t yet flooded they don’t buy flood insurance.

flooded park

So safety is not an increasing function of reliability, not always. It might dip down before going up. There’s a valley between unreliable and highly reliable where people are tempted to take unwise risks.

Artificial intelligence risks

I expect we’ll see a lot of this with artificial intelligence. Clumsy AI is not dangerous; pretty good AI is dangerous. Moderately reliable systems in general are dangerous, but this especially applies to AI.

As in the examples above, the better AI becomes, the more we rely on it. But there’s something else going on. As AI failures become less frequent they also become weird.

Adversarial attacks

You’ll see stories of someone putting a tiny sticker on a stop sign and now a computer vision algorithm thinks the stop sign is a frog or an ice cream sundae. In this case, there was a deliberate attack: someone knew how to design a sticker to fool the algorithm. But strange failures can also happen unprompted.

Unforced errors

Amazon’s search feature, for example, is usually very good. Sometimes I’ll get every word in a book title wrong and yet it will figure out what I meant. But one time I was searching for the book Universal Principles of Design.

I thought I remembered a “25” in the title. The subtitle turns out to be “125 ways to enhance reliability …” I searched on “25 Universal Design Principles” and the top result was a massage machine that will supposedly increase the size of a woman’s breasts. I tried the same search again this morning. The top result is a book on design. The next five results are

a clip-on rear view mirror
a case of adult diapers
a ratchet adapter socket
a beverage cup warmer, and
a folding bed.

The book I was after, and whose title I remembered pretty well, was nowhere in the results.

Because AI is literally artificial, it makes mistakes no human would make. If I went to a brick-and-mortar book store and told a clerk “I’m looking for a book. I think the title is something like ’25 Universal Design Principles,” the clerk would not say “Would you like to increase your breast size? Or maybe buy a box of diapers?”

In this case, the results were harmless, even entertaining. But unexpected results in a mission-critical system would not be so entertaining. Our efforts to make systems fool-proof has been based on experience with human fools, not artificial ones.

[1] This post is an elaboration on what started as a Twitter thread.

[2] I’m told that in Norway electrical power is very reliable, but also very dependent on electricity, including for heating. Alternative sources of fuel such as propane are hard to find.

How to compute the area of a polygon

Posted on 26 September 2018 by John

If you know the vertices of a polygon, how do you compute its area? This seems like this could be a complicated, with special cases for whether the polygon is convex or maybe other considerations. But as long as the polygon is “simple,” i.e. the sides meet at vertices but otherwise do not intersect each other, then there is a general formula for the area.

The formula

List the vertices starting anywhere and moving counterclockwise around the polygon: (x₁, y₁), (x₂, y₂), …, (x_n, y_n). Then the area is given by the formula below.

$A = \frac{1}{2} \begin{vmatrix} x_1 & x_2 & \ldots & x_n & x_1\\ y_1 & y_2 & \ldots & y_n & y_1 \end{vmatrix}$

But what does that mean? The notation is meant to be suggestive of a determinant. It’s not literally a determinant because the matrix isn’t square. But you evaluate it in a way analogous to 2 by 2 determinants: add the terms going down and to the right, and subtract the terms going up and to the right. That is,

x₁ y₂ + x₂ y₃ + … x_n y₁ – y₁ x₂ – y₂ x₃ – … – y_n x₁

This formula is sometimes called the shoelace formula because the pattern of multiplications resembles lacing a shoe. It’s also called the surveyor’s formula because it’s obviously useful in surveying.

Numerical implementation

As someone pointed out in the comments, in practice you might want to subtract the minimum x value from all the x coordinates and the minimum y value from all the y coordinates before using the formula. Why’s that?

If you add a constant amount to each vertex, you move your polygon but you don’t change the area. So in theory it makes no difference whether you translate the polygon before computing its area. But in floating point, it can make a difference.

The cardinal rule of floating point arithmetic is to avoid subtracting nearly equal numbers. If you subtract two numbers that agree to k significant figures, you can lose up to k figures of precision. We’ll illustrate this by taking a right triangle with base 2 and height π. The area should remain π as we translate the triangle away from the origin, we lose precision the further out we translate it.

Here’s a Python implementation of the shoelace formula.

    def area(x, y):
        n = len(x)
        s = 0.0
        for i in range(-1, n-1):
            s += x[i]*y[i+1] - y[i]*x[i+1]
        return 0.5*s

If you’re not familiar with Python, a couple things require explanation. First, range(n-1) is a list of integers starting at 0 but less than n-1. Second, the -1 index returns the last element of the array.

Now, watch how the precision of the area degrades as we shift the triangle by powers of 10.

    import numpy as np

    x = np.array([0.0, 0.0, 2.0])
    y = np.array([np.pi, 0.0, 0.0])

    for n in range(0, 10):
        t = 10**n
        print( area(x+t, y+t) )

This produces

    3.141592653589793
    3.1415926535897825
    3.1415926535901235
    3.1415926535846666
    3.141592651605606
    3.1415929794311523
    3.1416015625
    3.140625
    3.0
    0.0

Shifting by small amounts doesn’t make a noticeable difference, but we lose between one and two significant figures each time we increase t by a multiple of 10. We only have between 15 and 16 significant figures to start with in a standard floating point number, and so eventually we completely run out of significant figures.

When implementing the shoelace formula, we want to do the opposite of this example: instead of shifting coordinates so that they’re similar in size, we want to shift them toward the origin so that they’re not similar in size.

Passwords and power laws

Posted on 25 September 2018 by John

According to this paper [1], the empirical distribution of real passwords follows a power law [2]. In the authors’ terms, a Zipf-like distribution. The frequency of the rth most common password is proportional to something like 1/r. More precisely,

f_r = C r^–s

where s is on the order of 1. The value of s that best fit the data depended on the set of passwords, but their estimates of s varied from 0.46 to 0.91.

This means that the most common passwords are very common and easy to guess.

Size of password spaces

If passwords come from an alphabet of size A and have length n, then there are Aⁿ possibilities. For example, if a password has length 10 and consists of uppercase and lowercase English letters and digits, there are

62¹⁰ = 839,299,365,868,340,224

possible such passwords. If users chose passwords randomly from this set, brute force password attacks would be impractical. But brute force attacks are practical because passwords are not chosen uniformly from this large space of possibilities, far from it.

Attackers do not randomly try passwords. They start with the most common passwords and work their way down the list. In other words, attackers use Pareto’s rule.

Rules requiring, say, one upper case letter, don’t help much because most users will respond by using exactly one upper case letter, probably the first letter. If passwords must have one special character, most people will use exactly one special character, most likely at the end of the word. Expanding the alphabet size A exponentially increases the possible passwords, but does little to increase the actual number of passwords.

What’s interesting about the power law distribution is that there’s not a dichotomy between naive and sophisticated users. If there were, there would be a lot of common passwords, and all the rest uniformly distributed. Instead, there’s a continuum between the most naive and most sophisticated. That means a lot of people are not exactly naive, but not as secure as they think they are.

Some math

Under the Zipf model [3], the number of times we’d expect to see the most common password is NC where N is the size of the data set. The constant C is what it has to be for the frequencies to sum to 1. That is, C depends on the number of data points N and the exponent s and is given by

$C_{N, s} = \frac{1}{\sum_{r=1}^n r^{-s}}$

We can bound the sum in the denominator from above and below with integrals, as in the integral test for series convergence. This gives us a way to see more easily how C depends on its parameters.

$1 + \int_1^N x^{-s} \, dx > \frac{1}{C} = \sum_{r=1}^N r^{-s} > 1 + \int_2^{N+1} x^{-r}\, dx$

When s = 1 this reduces to

$1 + \log(N) > \frac{1}{C} > 1 + \log(N+1) - \log(2)$

and otherwise to

$1 + \frac{N^{1-s} - 1}{1-s} > \frac{1}{C} > 1 + \frac{(N+1)^{1-s} - 2^{1-s}}{1-s}$

Suppose you have N = 1,000,000 passwords. The range of s values found by Wang et al varied from roughly 0.5 to 0.9. Let’s first set s = 0.5. Then C is roughly 0.0005. This mean the most common password appears about 500 times.

If we keep N the same and set s to 0.9, then C is roughly 0.033, and so the most common password would appear about 33,000 times.

How to choose passwords

If you need to come up with a password, randomly generated passwords are best, but hard to remember. So people either use weak but memorable passwords, or use strong passwords and don’t try to remember them. The latter varies in sophistication from password management software down to Post-it notes stuck on a monitor.

One compromise is to concatenate a few randomly chosen words. Something like “thebestoftimes” would be weak because they are consecutive words from a famous novel. Something like “orangemarbleplungersoap” would be better.

Another compromise, one that takes more effort than most people are willing to expend, is to use Manuel Blum’s mental hash function.

[1] In case the link rots in the future, the paper is “Zipf’s Law in Passwords” by Ding Wang, Haibo Cheng, Ping Wang, Xinyi Huang, and Gaopeng Jian. IEEE Transactions on Information Forensics and Security, vol. 12, no. 11, November 2017.

[2] Nothing follows a power law distribution exactly and everywhere. But that’s OK: nothing exactly follows any other distribution everywhere: not Gaussian, not exponential, etc. But a lot of things, like user passwords, approximately follow a power law, especially over the middle of their range. Power law’s like Zipf’s law tend to not fit as well at the beginning and the end.

[3] Here I’m using a pure Zipf model for simplicity. The paper [1] used a Zipf-like model that I’m not using here. Related to the footnote [2] above, it’s often necessary to use a minor variation on the pure Zipf model to get a better fit.

Solving hard problems

Related posts

Related basic Jacobi functions to trig functions

The rest of the Jacobi functions

Ratios

Zeros and poles

Identities

Derivatives

Modulus, parameter, and modular angle

Quarter periods

Amplitude

Jacobi elliptic functions in Mathematica

Jacobi elliptic functions in Python

Related posts

Python implementation

Related posts

Related posts

Eddington’s constant

Related posts

Related posts

Related posts

Dangerously reliable systems

Artificial intelligence risks

Adversarial attacks

Unforced errors

The formula

Numerical implementation

Size of password spaces

Some math

How to choose passwords

Related posts