Solving hard problems

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

Crinkle Crankle Calculus

Posted on 19 November 2019 by John

A crinkle crankle wall

A crinkle crankle wall, also called a serpentine wall, is a wavy wall that may seem to sacrifice some efficiency for aesthetics. The curves add visual interest, but they use more material than a straight way. Except they don’t! They use more bricks than a straight wall of the same thickness but they don’t have to be as thick.

Crinkle crankle walls resist horizontal forces, like wind, more than straight wall would. So if the alternative to a crinkle crankle wall one-brick thick is a straight wall two or more bricks thick, the former saves material. How much material?

The amount of material used in the wall is proportional to the product of its length and thickness. Suppose the wall is shaped like a sine wave and consider a section of wall 2π long. If the wall is in the shape of a sin(θ), then we need to find the arc length of this curve. This works out to the following integral.

$\int_0^{2\pi} \sqrt{1 + a^2 \cos^2(x)}\, dx$

The parameter a is the amplitude of the sine wave. If a = 0, we have a flat wave, i.e. a straight wall, as so the length of this segment is 2π = 6.2832. If a = 1, the integral is 7.6404. So a section of wall is 22% longer, but uses 50% less material per unit length as a wall two bricks thick.

The integral above cannot be computed in closed form in terms of elementary functions, so this would make a good homework exercise for a class covering numerical integration.

The integral can be computed in terms of special functions. It equals 4 E(-a²) where E is the “complete elliptic integral of the second kind.” This function is implemented as EllipticE in Mathematica and as scipy.special.ellipe in Python.

As the amplitude a increases, the arc length of a section of wall increases. You could solve for the value of a to give you whatever arc length you like. For example, if a = 1.4422 then the length is twice that of a straight line. So a crinkle crankle wall with amplitude 1.4422 uses about as many bricks as a straight wall twice as thick.

***

Photo “Crinkle crankle wall, Fulbourn” by Bob Jones is licensed under CC BY-SA 2.0

Rook graphs and Paley graphs

Posted on 19 November 2019 by John

An m by n rook graph is formed by associating a node with each square of an m by n chess board and connecting two nodes with an edge if a rook can legally move from one to the other. That is, two nodes are connected if they are either in the same row or in the same column.

A Paley graph of order q is a graph with q nodes where q is a prime power and congruent to 1 mod 4. Label each node with an element of the finite group F with q elements and connect two nodes with an edge if the difference between their labels is a quadratic residue. That is, for nodes labeled a and b, connect these nodes if there exist an x in F such that a – b = x².

We’ll look at a graph which is both a rook graph and a Paley graph: the 3 by 3 rook graph is the same as the Paley graph of order 9.

Rook graph 3 x 3

Constructing the rook graph is easy. If we number our 3 by 3 chess board as follows

    1 2 3
    4 5 6
    7 8 9

then we connect 1 to 2, 3, 4, and 7. We connect 2 to 1, 3, 5, and 8. Etc. This gives us the following graph.

Rook graph

(The graph was made with Sketchviz introduced here.)

Paley graph of order 9

They Paley graph is a little more complicated to construct. First we need to explain what a field with 9 elements is.

A field with 3 elements is just integers mod 3. We can think of the field with 9 elements as the analog of complex numbers: numbers of the form a + bi where a and b are integers mod 3. As with the complex numbers, i² = -1.

When we square each of the nine elements of our field, we get four distinct non-zero results: 1, 2, i, and 2i. For example, i is one of our squares because

(1 + 2i)(1 + 2i) = (1 – 4) + 4i = i

Here we used -3 = 0 mod 3 and 4 = 1 mod 3.

Since our squares are 1, 2, i, and 2i, that says 0 is connected to these four values. Then we connect 1 to 2, 0, 1 + i and 1 + 2i because each of these numbers minus 1 is a square. If we keep doing this for all nine nodes, we get the following graph.

Paley graph

Isomorphism

It’s clear that the two graphs above are relabelings of each other. We can see this explicitly by relabeling our chess board as follows:

    0    i    2i
    1  1+i  1+2i
    2  2+i  2+2i

Moving horizontally adds i, which is a square, and moving vertically adds 1, which is a square. So the Paley connections are rook connections.

Uniqueness

Are there more rook graphs that are isomorphic to Paley graphs? No, this is the only one.

Every pair of vertices in a rook graph that are not neighbors have two common. If the points with coordinates (a, b) and (c, d) are on different rows and columns, then both are connected to (a, d) and (c, b).

In a Paley graph of order q, every pair of nodes that aren’t neighbors have (q – 1)/4 common neighbors. For a Paley graph to be isomorphic to a rook graph, we must have (q – 1)/4 = 2, and so q = 9.

NP vs small P

Posted on 18 November 2019 by John

The P vs NP conjecture has always seemed a little odd to me. Or rather the interpretation of the conjecture that any problem in P is tractable. How reassuring is to to know a problem can be solved in time less than some polynomial function of the size of the input if that polynomial has extremely high degree?

But this evening I ran across a fascinating comment by Avi Wigderson [1] that puts the P vs NP conjecture in a different light:

Very few known polynomial time algorithms for natural problems have exponents above 3 or 4 (even though at discovery the initial exponent may have been 30 or 40).

Problems in P may be more tractable in practice than in (current) theory. Wigderson’s comment suggests that if you can find any polynomial time algorithm, your chances are improved that you can find a small-order polynomial time algorithm. It seems there’s something deep going on here that would be hard to formalize.

[1] From his new book Mathematics and Computation

Inverse trig function implementations

Posted on 18 November 2019 by John

Programming languages differ in the names they use for inverse trig functions, and in some cases differ in their return values for the same inputs.

Choice of range

If sin(θ) = x, then θ is the inverse sine of x, right? Maybe. Depends on how you define the inverse sine. If -1 ≤ x ≤ 1, then there are infinitely many θ’s whose sine is x. So you have to make a choice.

The conventional choice is for inverse sine to return a value of θ in the closed interval

[-π/2, π/2].

Inverse tangent uses the same choice. However, the end points aren’t possible outputs, so inverse tangent typically returns a value in the open interval

(-π/2, π/2).

For inverse cosine, the usual choice is to return values in the closed interval

[0, π].

Naming conventions

Aside from how to define inverse trig functions, there’s the matter of what to call them. The inverse of tangent, for example, can be written tan^-1, arctan, or atan. You might even see arctg or atn or some other variation.

Most programming languages I’m aware of use atan etc. Python uses atan, but SciPy uses arctan. Mathematica uses ArcTan.

Software implementations typically follow the conventions above, and will return the same value for the same inputs, as long as inputs and outputs are real numbers.

Arctangent with two arguments

Many programming languages have a function atan2 that takes two arguments, y and x, and returns an angle whose tangent is y/x. Note that the y argument to atan2 comes first: top-down order, numerator then denominator, not alphabetical order.

However, unlike the one-argument atan function, the return value may not be in the interval (-π/2, π/2). Instead, atan2 returns an angle in the same quadrant as x and y. For example,

    atan2(-1, 1)

returns -π/4; because the point (1, -1) is in the 4th quadrant, it returns an angle in the 4th quadrant. But

    atan2(1, -1)

returns 3π/4. Here the point (-1, 1) and the angle 3π/4 are in the 2nd quadrant.

In Common Lisp, there’s no function named atan2; you just call atan with two arguments. Mathematica is similar: you call ArcTan with two arguments. However, Mathematica uses the opposite convention and takes x as its first argument and y as its second.

Complex values

Some software libraries support complex numbers and some do not. And among those that do support complex numbers, the return values may differ. This is because, as above, you have choices to make when extending these functions to complex inputs and outputs.

For example, in Python,

    math.asin(2)

returns a domain error and

    scipy.arcsin(2)

returns 1.5707964 + 1.3169578i.

In Common Lisp,

    (asin 2)

returns 1.5707964 – 1.3169578i. Both SciPy and Common Lisp return complex numbers whose sine equals 2, but they follow different conventions for what number to chose. That is, both

    scipy.sin(1.5707964 - 1.3169578j)
    scipy.sin(1.5707964 + 1.3169578j)

in Python and

    (sin #C(1.5707964 -1.3169578))
    (sin #C(1.5707964 +1.3169578))

in Common Lisp both return 2, aside from rounding error.

In C99, the function casin (complex arc sine) from <complex.h>

    double complex z = casin(2);

returns 1.5707964 + 1.3169578i. The Mathematica call

    ArcSin[2.]

returns 1.5707964 – 1.3169578i.

Branch cuts

The Common Lisp standardization committee did a very careful job of defining math functions for complex arguments. I’ve written before about this here. You don’t need to know anything about Lisp to read that post.

The committee ultimately decided to first rigorously define the two-argument arctangent function and bootstrap everything else from there. The post above explains how.

Other programming language have made different choices and so may produce different values, as demonstrated above. I mention the Common Lisp convention because they did a great job of considering every detail, such as how to handle +0 and -0 in IEEE floating point.

arctan( k tan(x) )

Posted on 18 November 2019 by John

I recently had to work with the function

f(x; k) = arctan( k tan(x) )

The semicolon between x and k is meant to imply that we’re thinking of x as the argument and k as a parameter.

This function is an example of the coming full circle theme I’ve written about several times. Here’s how a novice, a journeyman, and an expert might approach studying our function.

Novice: arctan( k tan(x) ) = kx.
Journeyman: You can’t do that!
Expert: arctan( k tan(x) ) ≈ kx for small x.

Novices often move symbols around without thinking about their meaning, and so someone might pull the k outside (why not?) and notice that arctan( tan(x) ) = x.

Someone with a budding mathematical conscience might conclude that since our function is nonlinear in x and in k that there’s not much that can be said without more work.

Someone with more experience might see that both tan(x) and arctan(x) have the form x + O(x³) and so

arctan( k tan(x) ) ≈ kx

should be a very good approximation for small x.

Here’s a plot of our function for varying values of k.

Plot of atan( k tan(x) ) for varying k

Each is fairly flat in the middle, with slope equal to its value of k.

As k increases, f(x; k) becomes more and more like a step function, equal to -π/2 for negative x and π/2 for positive x, i.e.

arctan( k tan(x) ) ≈ sgn(x) π/2

for large k. Here again we might have discussion like above.

Novice: Set k = ±∞. Then ±∞ tan(x) = ±∞ and arctan(±∞) = ±π/2.
Journeyman: You can’t do that!
Expert: Well, you can if you interpret everything in terms of limits.

Orbital resonance in Neptune’s moons

Posted on 17 November 2019 by John

Phys.com published an article a couple days ago NASA finds Neptune moons locked in ‘dance of avoidance’. The article is based on the scholarly paper Orbits and resonances of the regular moons of Neptune.

The two moons closest to Neptune, named Naiad and Thalassa, orbit at nearly the same distance, 48,224 km for Naiad and 50,074 km for Thalassa. However, the moons don’t come as close to each other as you would expect from looking just at the radii of their orbits.

Although the radii of the orbits differ by 1850 km, the two moons do not get closer than about 2800 km apart. The reason has to do with the inclination of the two orbital planes with each other, and a resonance between their orbital periods. When the moons approach each other, one dips below the other, increasing their separation.

Assume the orbits of both moons are circular. (They very nearly are, with major and minor axes differing by less than a kilometer.) Also, choose a coordinate system so that Thalassa orbits in the xy plane. The position of Thalassa at time t is

r_T (cos 2πt/T_T, sin 2πt/T_T, 0)

where r_T is the radius of Thalassa’s orbit and T_T is its period.

The orbit of Naiad is inclined to that of Thalassa by an angle u. The position of Naiad at time t is

r_N (cos u cos 2πt/T_N, sin 2πt/T_N, -sin u cos 2πt/T_N).

I tried implementing the model above using data from Wikipedia articles on the two moons, but in my calculations the moons get closer than reported. I suspect the parameters have to be set carefully in order to demonstrate the kind of resonance we see in observation, and we don’t know these parameters accurately enough to make a direct geometric approach like this work.

from numpy import *

# Orbital radii in km 
r_T = 50074.44
r_N = 48224.41

# Period in days
T_T = 0.31148444
T_N = 0.29439580

def thalassa(t):

    # frequency
    f = 2*pi/T_T

    return r_T * array([
        cos(f*t),
        sin(f*t),
        0
    ])

def naiad(t):
    
    # inclination in radians relative to Thalassa
    i = 0.082205
    
    # frequency
    f = 2*pi/T_N

    return r_N * array([
        cos(i)*cos(f*t),
        sin(f*t),
        -sin(i)*cos(f*t)
    ])

def dist(t):
    return sum((naiad(t) - thalassa(t))**2)**0.5

d = [dist(t) for t in linspace(0, 10*73*T_N, 1000)]
print(d[0], min(d))

In this simulation, the moons are 4442 km apart when t = 0, but this is not the minimum distance. The code above shows an instance where the distance is 2021 km. I tried minor tweaks, such as adding a phase shift to one of the planets or varying the angle of inclination, but I didn’t stumble on anything that cane closer to reproducing the observational results. Maybe resonance depends on factors I’ve left out. Naiad and Thalassa are practically massless compared to Neptune, but perhaps the simulation needs to account for their mass.

The periods of the two moons are nearly in a ratio of 69 to 73. We could idealize the problem by making the periods rational numbers in exactly a ratio of 69 to 73. Then the system would be exactly periodic, and we could find the minimum distance. Studying a model system close to the real one might be useful. If we do tweak the periods, we’d need to tweak the radii as well so that Kepler’s third law applies, i.e. we’d need to keep the squares of the periods proportional to the cubes of the radii.

Related post: Average distance between planets

The hardest logarithm to compute

Posted on 16 November 2019 by John

Suppose you want to compute the natural logarithms of every floating point number, correctly truncated to a floating point result. Here by floating point number we mean an IEEE standard 64-bit float, what C calls a double. Which logarithm is hardest to compute?

We’ll get to the hardest logarithm shortly, but we’ll first start with a warm up problem. Suppose you needed to compute the first digit of a number z. You know that z equals 2 + ε where ε is either 10^-100 or -10^-100. If ε is positive, you should return 2. If ε is negative, you should return 1. You know z to very high precision a priori, but to return the first digit correctly, you need to compute z to 100 decimal places.

In this post we’re truncating floating point numbers, i.e. rounding down, but similar considerations apply when rounding to nearest. For example, if you wanted to compute 0.5 + ε rounded to the nearest integer.

In general, it may not be possible to know how accurately you need to compute a value in order to round it correctly. William Kahan called this the table maker’s dilemma.

Worst case for logarithm

Lefèvre and Muller [1] found that the floating point number whose logarithm is hardest to compute is

x = 1.0110001010101000100001100001001101100010100110110110 × 2⁶⁷⁸

where the fractional part is written in binary. The log of x, again written in binary, is

111010110.0100011110011110101110100111110010010111000100000000000000000000000000000000000000000000000000000000000000000111 …

The first 53 significant bits, all the bits a floating point number can represent [2], are followed by 65 zeros. We have to compute the log with a precision of more than 118 bits in order to know whether the last bit of the fraction part of the float representation should be a 0 or a 1. Usually you don’t need nearly that much precision, but this is the worst case.

Mathematica code

Here is Mathematica code to verify the result above. I don’t know how to give Mathematica floating point numbers in binary, but I know how to give it integers, so I’ll multiply the fraction part by 2⁵² and take 52 from the exponent.

n = FromDigits[
  "10110001010101000100001100001001101100010100110110110", 2]
x = n  2^(678 - 52)
BaseForm[N[Log[x], 40], 2]

This computes the log of x to 40 significant decimal figures, between 132 and 133 bits, and gives the result above.

A few days ago I wrote about another example by Jean-Michel Muller: A floating point oddity.

William Kahan also came up recently in my post on summing an array of floating point numbers.

***

[1] Vincent Lefèvre, Jean-Michel Muller. Worst Cases for Correct Rounding of the Elementary Functions in Double Precision. November 2000. Research Report 2000-35. Available here.

[2] There are 52 bits in the the faction of an IEEE double. The first bit of a fraction is 1, unless a number is subnormal, so the first bit is left implicit, squeezing out one extra bit of precision. That is, storing 52 bits gives us 53 bits of precision.

Hand-drawn Graphviz diagrams

Posted on 13 November 2019 by John

The Sketchviz site lets you enter GraphViz code and get out something that looks hand-drawn. I decided to try it on some of the GraphViz diagrams I’ve made for this blog.

Here’s a diagram from my post on Rock, Paper, Scissors, Lizard, Spock.

The Sketchviz site defaults to Tinos font, which does not look hand-drawn. Maybe you’d want that sometimes: neatly typeset labels and hand-drawn boxes and arrows. In my case, I thought it looked better to use Sedgwick Ave which looks closer to handwriting.

***

Some posts with GraphViz diagrams:

Some of these posts include GraphViz source.

Permutation distance

Posted on 13 November 2019 by John

If you have two sequences, and one is a permutation of the other, how do you measure how far apart the two sequences are?

This post was motivated by the previous post on anagram statistics. For a certain dictionary, the code used in that post found that the longest anagram pairs were the following:

certifications, rectifications
impressiveness, permissiveness
teaspoonsful, teaspoonfuls

Anagrams are more amusing when the two words are more dissimilar, at least in my opinion.

There are several ways to measure (dis)similarity. The words in the first two pairs above are dissimilar in meaning, but the words in the last pair are synonyms [1]. The pairs of words above are not that dissimilar as character strings.

It would be hard to quantify how dissimilar two words are in meaning, but easier to quantify how dissimilar they are as sequences of characters. We’ll look at two approaches: Hamming distance and transposition distance.

Hamming distance

The Hamming distance between two words is the number of positions in which they differ. We can see the Hamming distances between the words above by viewing them in a fixed-width font.

certifications
rectifications

impressiveness
permissiveness

teaspoonsful
teaspoonfuls

The Hamming distances above are 2, 5, and 4.

The anagrams with maximum Hamming distance in my dictionary are “imperfections” and “perfectionism.” (There’s something poetic about these words being anagrams.) The Hamming distance is 13 because these 13-letter words differ in every position.

imperfections
perfectionism

Transposition distance

Another way to measure the distance between two permutations is how many elements would have to be swapped to turn one list into the other. This is called the transposition distance. For example, the transposition distance between “certifications” and “rectifications” is 1, because you only need to swap the first and third characters to turn one into the other.

It’s not so obvious what the transposition distance is between “impressiveness” and “permissiveness”. We can find an upper bound on the distance by experimentation. The two words only differ in the first five letters, so the distance between “impressiveness” and “permissiveness” is no larger than the distance between “impre” and “permi”. The latter pair of words are no more than 4 transpositions apart as the following sequence shows:

impre
pmire
peirm
perim
permi

But is there a shorter sequence of transpositions? Would it help if we used the rest of the letters “ssiveness” as working space?

Computing transposition distance is hard, as in NP-hard. This is unfortunate since transposition has more practical applications than just word games. It is used, for example, in genetic research. It may be practical to compute for short sequences, but in general one must rely on bounds and approximations.

Note that transposition distance allows swapping any two elements. If we only allowed swapping consecutive elements, the problem is much easier, but the results are not the same. When restricted to consecutive swaps, the distance between “certifications” and “rectifications” is 3 rather than 1. We can swap the “c” and the “r” to turn “cer” into “rec” so we have to do something like

cer
ecr
erc
rec

I don’t know a name for the distance where you allow rotations. The words “teaspoonsful” and “teaspoonfuls” differ by one rotation, turning “sful” into “fuls.” While this can be done in one rotation, it would take several swaps.

[1] Although “teaspoonfuls” is far more common, I remember a schoolmarm teaching us that this is incorrect.

And I still recoil at “brother in laws” and say “brothers in law” instead; the majority is on my side on this one.

Anagram frequency

Posted on 13 November 2019 by John

An anagram of a word is another word formed by rearranging its letters. For example, “restful” and “fluster” are anagrams.

How common are anagrams? What is the largest set of words that are anagrams of each other? What are the longest words which have anagrams?

You’ll get different answers to these questions depending on what dictionary you use. I started with the american-english dictionary from the Linux package wamerican. This list contains 102,305 words. I removed all words containing an apostrophe so that possessive forms didn’t count as separate words. I then converted all words to lower case and removed duplicates. This reduced the list to 72,276 words.

Next I alphabetized the letters in each word to create a signature. Signatures that appear more than once correspond to anagrams.

Here are the statistics on anagram classes.

    |------+-------|
    | size | count |
    |------+-------|
    |    1 | 62093 |
    |    2 |  3600 |
    |    3 |   646 |
    |    4 |   160 |
    |    5 |    61 |
    |    6 |    13 |
    |    7 |     2 |
    |    8 |     1 |
    |------+-------|

This means that 62,093 words or about 86% are in an anagram class by themselves. So about 14% of words are an anagram of at least one other word.

The largest anagram class had eight members:

least
slate
stael
stale
steal
tales
teals
tesla

Stael is a proper name. Tesla is a proper name, but it is also a unit of magnetic induction. In my opinion, tesla should count as an English word and Stael should not.

My search found two anagram classes of size seven:

pares
parse
pears
rapes
reaps
spare
spear

and

carets
caster
caters
crates
reacts
recast
traces

The longest words in this dictionary that form anagrams are the following, two pair of 14-letter words and one pair of 12-letter words.

certifications, rectifications
impressiveness, permissiveness
teaspoonsful, teaspoonfuls

Dictionary of anagrams

I made a dictionary of anagrams here. Every word which has a anagram is listed, followed by its anagrams. Here are the first few lines:

abby: baby
abeam: ameba
abed: bade, bead
abel: able, bale, bela, elba
abet: bate, beat, beta
abets: baste, bates, beast, beats, betas
abetter: beretta
abhorred: harbored

There is some redundancy in this dictionary for convenience: every word in the list of anagrams will also appear as the first entry on a line.

Here’s the Python code that produced the dictionary.

from collections import defaultdict

lines = open("american-english", "r").readlines()

words = set()
for line in lines:
    if "'" not in line:
        line = line.strip().lower()
        words.add(line)

def sig(word):
    return "".join(sorted(word))

d = defaultdict(set)
for w in words:
    d[sig(w)].add(w)

for w in sorted(words):
    anas = sorted(d[sig(w)])
    if len(anas) > 1:
        anas.remove(w)
        print("{}: {}".format(w, ", ".join(anas)))

Words with the most consecutive vowels.

Solving hard problems

Crinkle Crankle Calculus

Rook graphs and Paley graphs

Rook graph 3 x 3

Paley graph of order 9

Isomorphism

Uniqueness

Related posts

NP vs small P

Related posts

Inverse trig function implementations

Choice of range

Naming conventions

Arctangent with two arguments

Complex values

Branch cuts

arctan( k tan(x) )

Related posts

Orbital resonance in Neptune’s moons

The hardest logarithm to compute

Worst case for logarithm

Mathematica code

Related posts

Hand-drawn Graphviz diagrams

Permutation distance

Hamming distance

Transposition distance

Related posts

Anagram frequency

Dictionary of anagrams

Related post