I’m sure other people have written about this, but I couldn’t find it. Maybe lots of people have written about this in Japanese but not many in English.

Japanese kana consists of two syllabaries, hiragana and katakana, that are like phonetic alphabets. Each has 46 basic characters, and each corresponds to a block of 96 Unicode characters. I had two simple questions:

- How do the 46 characters map into the 90 characters?
- Do they map the same way for both hiragana and katakana?

I’ll start with the second question because it’s easier. Hiragana and katakana are different ways of representing the same sounds, and they correspond one to one. For example, the full name of U+3047 (ぇ) is

HIRAGANA LETTER SMALL E

and the full name of its katakana counterpart U+30A7 (ェ) is

KATAKANA LETTER SMALL E

The only difference as far as Unicode goes is that katakana has three code points whose hiragana counterpart is unused, but these are not part of the *basic* letters.

The following Python code shows that the names of all the characters are the same except for the name of the system.

from unicodedata import name unused = [0, 151, 152] # not in hiragana for i in range(0,63): if i in unused: continue h = name(chr(0x3040 + i)) k = name(chr(0x30a0 + i)) assert(h == k.replace("KATAKANA", "HIRAGANA")) print("done")

You’ll see kana written in grid with one side labeled with 5 vowels and the other labeled with 10 consonants called a gojūon (五十音). That’s 50 cells, and in fact gojūon literally means 50 sounds, so how do we get 46? Five cells are empty, and one letter doesn’t fit into the grid. The empty cells are unused or archaic, and the extra character doesn’t fit the grid structure.

In the image below, the table on the left is for hiragana and the table on the right is for katakana. HTML versions of the tables available here.

Left out of each table is ん in hiragana and ン in katakana.

So does each set of 46 characters map into its Unicode code block?

Unicode numbers the letters consecutively if you traverse the grid increasing vowels first, then consonants, and adding the straggler at the end. But the reason 46 letters expand into more code points is that each letter can have one, two, or three variations. And there are various miscellaneous other symbols in the Unicode block.

For example, there is a LETTER E as well as the SMALL LETTER E mentioned above. Other variations seem to correspond to voiced and unvoiced versions of a consonant with a phonetic marker added to the voiced version. For example, く is U+304F, HIRAGANA LETTER KU, and ぐ is U+3050, HIRAGANA LETTER GU.

Here is how hiragana maps into Unicode. Each cell should be U+3000 plus the characters show.

a i u e o 42 44 46 48 4A k 4B 4D 4F 51 53 s 55 57 59 5B 5D t 5F 61 64 66 68 n 6A 6B 6C 6D 6E h 6F 72 75 78 7B m 7E 7F 80 81 82 y 84 86 88 r 89 8A 8B 8C 8D w 8F 92

The corresponding table for katakana is the previous table plus 0x60:

a i u e o A2 A4 A6 A8 AA k AB AD AF B1 B3 s B5 B7 B9 BB BD t BF C1 C4 C6 C8 n CA CB CC CD CE h CF D2 D5 D8 DB m DE DF E0 E1 E2 y E4 E6 E8 r E9 EA EB EC ED w EF F2

In each case, the letter missing from the table is the next consecutive value after the last in the table, i.e. ン is U+30F3.

In a Latin square of size *n* you have to assign one of *n* symbols to each cell so that each symbol appears exactly once in each row and column. A Room square is sort of a Latin square with pairs.

Each cell of a Room square is either empty or contains an unordered pair of symbols. Each symbol appears exactly once in each row and column, and every possible pair of symbols appears in some cell.

The following graphic shows a 7 × 7 Room square with eight symbols, each represented by a different color.

If you have trouble seeing some of the colors, take a look at the code below that made the image.

An *n* × *n* Room square corresponds to a Round Robin tournament schedule with *n* + 1 players. Each row represents a location (or *room*), and each column represents a round. Each player plays one game in each round, and each pair of players plays in exactly one location.

Here’s the code I wrote to create the image above.

import matplotlib.pyplot as plt colors = ["red", "orange", "yellow", "green", "blue", "black", "gray", "purple"] up, dn = True, False def draw_triangle(row, col, value, pos): if pos == up: x = [col, col, col+1] y = [row, row+1, row+1] else: x = [col, col+1, col+1] y = [row, row, row+1] plt.fill(x, y, colors[value]) sq = [ (1, 3, 0, 4), (1, 5, 3, 5), (1, 6, 1, 2), (1, 7, 7, 6), (2, 2, 6, 3), (2, 4, 2, 4), (2, 5, 0, 1), (2, 6, 7, 5), (3, 1, 5, 2), (3, 3, 1, 3), (3, 4, 6, 0), (3, 5, 7, 4), (4, 2, 0, 2), (4, 3, 5, 6), (4, 4, 7, 3), (4, 7, 4, 1), (5, 1, 6, 1), (5, 2, 4, 5), (5, 3, 7, 2), (5, 6, 3, 0), (6, 1, 3, 4), (6, 2, 7, 1), (6, 5, 2, 6), (6, 7, 5, 0), (7, 1, 7, 0), (7, 4, 1, 5), (7, 6, 4, 6), (7, 7, 2, 3) ] for t in sq: draw_triangle(t[0], t[1], t[2], up) draw_triangle(t[0], t[1], t[3], dn) plt.grid() plt.gca().set_aspect("equal") plt.show()

**Related post**: Balanced tournament designs

HCPCS stands for Healthcare Common Procedure Coding System. HCPCS codes are commonly pronounced like “hick pics.” I was curious/afraid to see what DALL-E would create with the prompt “HCPCS hick pics” and was pleasantly surprised that it produced the image above. More on that latter.

I occasionally need to search for medical codes—ICD-9 codes, ICD-10 codes, SNOMED codes, etc.—including HCPCS codes. Here’s some data that is supposed to contain a certain kind of code; does it? Or, more importantly, this text is supposed to have been scrubbed of medical codes; was it?

The most accurate way to search for HCPCS codes would be to have a complete list of codes and search for each one. There are a couple reasons why this isn’t so useful. First of all, I’m doing a quick scan, not writing data validation software. So a quick-and-dirty regular expression is fine. In fact, it’s better: I may be more interested in finding things that *look like* HCPCS codes than actual HCPCS codes. The former would include some typos or codes that were added after my script was written.

Another advantage of using a regular expression that approximates HCPCS codes is that HCPCS codes are copyrighted. I don’t know whether a script using a list of HCPCS codes would be fair use, but it doesn’t matter if I’m using a regular expression.

According to Wikipedia, “CPT is a registered trademark of the American Medical Association, and its largest single source of income.” CPT is part of HCPCS, namely Level I. It is controversial that CPT codes are a commercial product whose use is required by law.

There are three levels of HCPCS codes. Level I is synonymous with CPT and consists of five-digit numbers. It’s easy to search for five-digit numbers, say with the regular expression `\d{5}`

, but of course many five-digit numbers are not medical codes. But if text that has supposedly been stripped of medical codes (and zip codes) has a lot of five-digit numbers, it could be something to look into.

Level II and Level III codes are five-character codes consisting of four digits and either an F or a T. So you might search on `\d{4}[FT]`

. While a five-digit number is likely to be a false positive when searching for HCPCS codes, four-digits following by an F or T is more likely to be a hit.

I edited the image that DALL-E produced two ways. First, I cropped it vertically; the original image had a lot more grass at the bottom.

Second, I lowered the resolution. The images DALL-E generates are very large. Even after cutting off the lower half of the image, the remaining part was over a megabyte. I used Squoosh to reduce the image to 85 kb. Curiously, when I opened the squooshed output in GIMP and lowered the resolution further, the image actually got larger, so went back to the output from Squoosh.

The image looks more natural than anything I’ve made using DALL-E. I look at the image and think it would be nice to be there. But if you look at the barn closely, it has a little of that surreal quality that DALL-E images usually have. This is not a photograph or a painting.

This discussion of AI-generated images is a bit of a rabbit trail, but it does tie in with searching for medical codes. Both are a sort of forensic criticism: is this what it’s supposed to be or is something out of place? When looking for errors or fraud, you have to look at data similarly to the way you might analyze an image.

There are codes for representing tiles horizontally or vertically. And even though, for example, the 5-3 is the same domino as the 3-5, there are separate characters for representing the orientation of the tile: one for 3 on the left and one for 5 on the left.

When you include orientation like this, a domino becomes essentially a **base 7 number**: the number of spots on one end is the number of 7s and the number of spots on the other end is the number of 1s. And the order of the characters corresponds to the order as base 7 numbers:

0-0, 0-1, 0-2, …, 1-0, 1-1, 1-2, … 6-6.

The horizontal dominoes start with the double blank at U+1F031 and the vertical dominoes start with U+1F063, a difference of 32 in hexadecimal or 50 in base 10. So you can rotate a domino tile by adding or subtracting 50 to its code point.

The following tiny Python function gives the codepoint for the domino with *a* spots on the left (or top) and *b* spots on the right (or bottom).

def code(a, b, wide): cp = 0x1f031 if wide else 0x1f063 return cp + 7*a + b

We can use this function to print a (3, 5) tile horizontally and a (6, 4) tile vertically.

print( chr(code(3, 5, True )), chr(code(6, 4, False)) )

To my surprise, my computer had the fonts installed to display the results. This isn’t guaranteed for such high Unicode values.

The post Dominoes in Unicode first appeared on John D. Cook.]]>The problem with raging against the machine is that the machine has learned to feed off rage.

Someone appropriately replied with a screenshot from an episode of Star Trek TOS, Day of the Dove, about a being that feeds off anger, like contemporary media.

The post Raging against the machine first appeared on John D. Cook.]]>In my first post on Costas arrays I mentioned in a footnote that Lempel’s algorithm works more generally over any finite field, but I stuck to integers modulo primes to make things simpler.

In this post we’ll create a 7 × 7 Costas array by working in a field with 9 elements. Along the way we’ll explain in great detail how finite field arithmetic works, using the field with 9 elements as an example.

The order of a finite field can be any power of a prime. That is, a finite field can have *q* elements if (and only if) *q* = *p ^{k}* for some prime

If *q* is a prime, then any field with *q* elements is isomorphic to the integers mod *q*. Addition and multiplication are defined just as in the integers, except you always take the remainder by *q* after you add or multiply.

But if *q* is a prime *power*, things are different. So while multiplication in a field of 7 elements is simply multiplication mod 7, multiplication in a field of 9 elements is **not** multiplication mod 9.

The field with nine elements can be defined as polynomials of the form *ax* + *b* where *a* and *b* are integers mod 3, i.e. *a* and *b* can take on the values 0, 1, or 2.

You can define addition in this little field the same way you always define polynomial addition, with the understanding that the coefficients are added mod 3. So, for example,

(2*x* + 1) + (*x* + 0) = 1

because the leading coefficient is 2 + 1, which reduces to 0 mod 3.

Multiplication is also defined as polynomial multiplication, with coefficient arithmetic being carried out mod 3, with one more step: after multiplying the polynomials, you divide by an irreducible quadratic polynomial and keep the remainder.

In general, if *q* = *p ^{k}* for some prime

The definition of multiplication brings up several questions.

- What is an irreducible polynomial?
- Can I always find one?
- How can I find one?
- Which one should I pick?

An irreducible polynomial is one that cannot be factored in the given context. Here we want a polynomial of degree 2, with coefficient in integers mod 3, that cannot be factored as the product of two linear polynomials with coefficients in integers mod 3.

There are plenty of irreducible polynomials. In our case, we have 36 to choose from. (This post explains a formula for counting how many irreducible polynomials there are of a given order over a given finite field.)

There are algorithms for finding irreducible polynomials, but this is beyond the scope of this post. Here’s one we will use:

*x*² + *x* + 2.

You could verify by brute force that this polynomial cannot be factored into

(*ax* + *b*)(*cx* + *d*)

where *a* and c are non-zero. There are 2 × 3 × 2 × 3 possibilities for *a*, *b*, *c*, and *d*, and you could verify that none of them work.

We said there are 36 irreducible polynomials we could choose. What about the other 35 we didn’t choose? They all lead to the same field, where by “same” we mean *isomorphic*. They generate different *representations* of the field, but there’s a direct correspondence (an isomorphism) between any two representations.

In one sense this is nothing to worry about, and in another sense it is. You’re not going to create something fundamentally different by choosing a different irreducible polynomial. If you and I pick different polynomials, we will agree on all the algebraic properties of our fields, but we could disagree on what the product of two particular polynomials is, like two people using two different units of measure. This can be an issue when doing finite field arithmetic with software: two software packages could give different results for a given product.

An element *g* in in group *G* is primitive if every element of *G* is some power of *g*. For example, 3 is a primitive element of the integers mod 7 because if we take successive powers of 3 mod 7 we get

3, 2, 6, 4, 5, 1

and that’s all the non-zero elements of the integers mod 7. On the other hand, 2 is *not* a primitive element mod 7 because its powers are 2, 4, and 1. There’s no way to raise 2 to any power mod 8 and get 3, 5, or 6 out.

We need a primitive element of our finite field in order to carry out Lempel’s algorithm. It turns out we can use *x* as our primitive element. That is, every non-zero element of our field is a power of *x*. There are other primitive elements we could have chosen[1].

That statement may sound wrong. How, for example, can a power of *x* equal 2x + 1? We have to remember that we’re multiplying polynomials as described above, not as polynomials over real numbers.

When we multiply *x* by *x*, we get *x*² of course, but then we have to divide by our choice of irreducible polynomial, which we said above would be *x*² + *x* + 2.

If we divide *x*² by *x*² + *x* + 2 as we would in a high school algebra class, we get a quotient of 1 with a remainder of –*x* – 2. When we reduce the coefficients mod 3, we get 2*x* + 1. So in our world,

*x*² = 2*x* + 1.

OK, that takes some getting used to, but now what is *x*³? Well, it’s *x* times *x*² of course, but what is that? As before we’ll first pretend we’re doing high school algebra, then we’ll adjust. So *x* times 2*x* + 1 equals 2*x*² + *x*. When we divide by our irreducible polynomial, we get

2*x*² + *x* = 2(*x*² + *x* + 2) + 2*x* + 2

so the remainder is 2*x* + 2. Therefore in our context,

*x*³ = 2*x* + 2.

Lempel’s algorithm says to fill in the (*i*, *j*) square of our matrix if and only if

*a*^{i} + *a*^{j} = 1

where *a* is a generator of our field. In our case *a* is simply the polynomial *x*. Here *i* and *j* run from 1 through 7 (i.e. up to *q* – 2, and for us *q* = 9.)

We don’t need a full multiplication table; we only need powers of *x*. Here they are:

*x*- 2
*x*+ 1 - 2
*x*+ 2 - 2
- 2
*x* *x*+ 2*x*+ 1- 1

Incidentally, we could use this as a sort of table of logarithms. To multiply any two elements, scan the list to find out what power of *x* each element is. Then add the exponents, and see what element that corresponds to.

For each row *i*, we compute *x*^{i} and then find which value of *j* satisfies

*x*^{j} = 1 – *x*^{i}.

So here we go.

When *i* = 1, 1 – *x*^{i} = 1 – *x* = 2*x* + 1, and so *j* = 2.

When *i* = 2, 1 – *x*^{i} = 1 – *x*² = 1 – (2*x* + 1) = *x*, and so *j* = 1.

When *i* = 3, 1 – *x*^{i} = 1 – (2*x*+ 2) = *x* + 2, so *j* = 6.

When *i* = 4, 1 – *x*^{i} = 1 – 2 = 2, so *j* = 4.

When *i* = 5, 1 – *x*^{i} = 1 – 2*x* = *x* + 1, so *j* = 7.

We can figure out the rest by symmetry. As noted in the previous post, the pairs (*i*, *j*) are symmetric. So since (3, 6) is filled in, so is (6, 3). And since (5, 7) is filled in, so is (7, 5). We could have also figured out the case of *i* = 2 by symmetry.

Here we visualize the Costas array we’ve found as in the earlier post. The line segments between filled in cells are all unique, so we draw cell positions and the connecting lines.

Here we translate all the lines to the origin to show that they’re distinct.

No two lines are the same, but some lines have the same slope with different lengths. To make it possible to see these, I added a small amount of randomness to the points before plotting them.

***

[1] The number of primitive polynomials of degree *n* in a field with *q* elements is φ(*q*^{n} – 1)/*n* where φ is Euler’s totient function. In our example, *n* = 1 and *q* = 9, so there are φ(8) = 4 primitive polynomials to choose from. What are they?

Every non-zero element of our field is a power of *x*, so a primitive polynomial *y* must itself be a power of *x*, say *y* = *x*^{k}. If *y*^{j} = *x*, then we must have *x*^{jk} = *x*. The multiplicative group of our field has 8 elements, so *jk* must be congruent to 1 mod 8. This means *k* must be 1, 3, 5, or 7. So our generators are *x*, *x*^{3} = 2*x* + 2 = *x* + 2, *x*^{5} = 2*x*, and *x*^{7} = *x* + 1.

The earlier post implemented the Lempel algorithm in Python. Here we’ll implement it in Mathematica. Lempel’s algorithm says to start with parameters *p* and *a* where *p* is a prime and *a* is a primitive root mod *p* [1]. Then you can create a Costas array of size *n* = *p*-2 by filling in the (*i*, *j*) square if and only if

*a*^{i} + *a*^{j} = 1 mod *p*.

We can implement this in Mathematica by

fill[i_, j_, p_, a_] := If[Mod[a^i + a^j, p] == 1, 1, 0] m[p_, a_] := Table[fill[i, j, p, a], {i, 1, p - 2}, {j, 1, p - 2}]

We could visualize the Costas array generated by *p* = 11 and *a* = 2 using `ArrayPlot`

. The code

ArrayPlot[m[11, 2], Mesh -> All]

Generates the following image.

Note that the image is symmetric with respect to the main diagonal. That’s because our algorithm is symmetric in *i* and *j*. But Costas arrays are not always symmetrical, so this underscores the point that there is no known scalable algorithm for finding all Costas arrays.

In this example we started with knowing that 2 was a primitive root mod 11. We could have had Mathematica pick a primitive root mod *p* for us by calling `PrimitiveRoot[p]`

Just out of curiosity, let’s redo the example above, but instead of testing whether

*a*^{i} + *a*^{j} = 1 mod *p*.

we’ll just use *a*^{i} + *a*^{j} mod *p*.

The code

ArrayPlot[Table[Mod[2^i + 2^j, 11], {i, 1, 9}, {j, 1, 9}]]

produces the following image.

[1] Lempel’s algorithm generalizes to finite fields of any order, not just integers modulo a prime.

The post Costas arrays in Mathematica first appeared on John D. Cook.]]>The 2-letter abbreviations may be familiar because it is very often (but not always [1]) also the country code top-level domain (ccTLD). For example, AU is the ISO abbreviation for Australia, and .au is the ccTLD.

I was curious about the relation between the two-letter and three-letter abbreviations. Sometimes the former is a prefix of the latter, such as US and USA. Sometimes the latter adds a letter in the middle, such as going from CN to CHN. How common are these two patterns?

I wrote a script using the iso3166 Python module to find out.

Turns out that the prefix pattern is most common, and occurs 63% of the time.

The infix pattern is the next most common, occurring 29% of the time.

Suffixes are rare. There are only for instances where the 2-letter name is the last two letters of the 3-letter name: ATF, MYT, SPM, and SGS.

The remaining possibilities are miscellaneous relations, such as IL and ISR for Israel.

Here’s a table of countries and ISO codes in plain text (org-mode) format.

[1] There are four ccTLDs that are not ISO 3-166-1 alpha-2 names: uk, su, ac, and eu.

The post Two-letter vs Three-letter Country Abbreviations first appeared on John D. Cook.]]>`CountryData`

database contains flag descriptions. So I thought I’d use the flag descriptions to see which flags Mathematica things are similar.
For example, the `FlagDescription`

attribute for Chad in Mathematica is

Three equal vertical bands of blue (hoist side), yellow, and red; similar to the flag of Romania; also similar to the flags of Andorra and Moldova, both of which have a national coat of arms centered in the yellow band; design was based on the flag of France.

I had Mathematica output a list of countries and flag descriptions, then searched the output for the word “similar.” I then made the following groupings based on the output [1].

Each flag has an emoji, so here are the groupings above using emoji icons

[1] The groupings are based on Mathematica’s output, but I did some editing. Strictly following Mathematica’s descriptions would have been complicated. For example, Mathematica’s description might say A is similar to B, but not say B is similar to A. Or it might cluster four flags together that could better be split into two pairs.

The post Finding similar world flags with Mathematica first appeared on John D. Cook.]]>In this post we’re going to look at a similar problem, weakening one requirement and adding another. First we’re going to remove the requirement that no two pieces can be on the same diagonal. This turns the *n* queens problem into the *n* rooks problem because rooks cannot move diagonally.

Next, imagine running stiff wires from every rook to every other rook. We will require that no two wires have the same length and run in the same direction. in more mathematical terms, we require that the displacement vectors between all the rooks are unique.

A solution to this problem is called a Costas array. In 1965, J. P. Costas invented what are now called Costas arrays to solve a problem in radar signal processing.

Why do we call a solution a Costas *array* rather than a Costas *matrix*? Because a matrix solution can be described by recording for each column the row number of the occupied cell. For example, we could describe the eight queens solution above as

(2, 4, 6, 8, 3, 1, 7, 5).

Here I’m numbering rows from the bottom up, like points in the first quadrant, rather than top to bottom as one does with matrices.

Note that the solution to the eight queens problem above is **not** a Costas array because some of the displacement vectors between queens are the same. For example, if we number queens from left to right, the displacement vectors between the first and second queens is the same as the displacement vector between the second and third queens.

Here’s a visualization of a Costas array of order 5.

It’s clear from the plot that no two dots are in the same row or column, but it’s not obvious that all the connecting lines are different. To see that the lines are different, we move all the tails of the connecting lines to the origin, keeping their length and direction the same.

There are 10 colored lines in the first plot, but at first you may only see 9 in the second plot. That’s because two of the lines have the same direction but different length. If you look closely, you’ll see that there’s a short purple line on top of the long red line. That is, one line runs from the origin to (1, -1) and another from the origin to (4, -4).

Here is a visualization of a Costas array of order 9.

And here are its displacement vectors translated to the origin.

Here is a visualization of a Costas array of order 15.

And here are its displacement vectors translated to the origin.

There are algorithms for generating *some* Costas arrays, but not all. Every known algorithm [1] leaves out some solutions, and it is not know whether Costas arrays exist for some values of *n*.

The Costas arrays above were generated using the Lempel construction algorithm. Given a prime *p* [2] and a primitive root mod *p* [3], the following Python code will generate a Costas array of order *p* – 2.

p = 11 # prime a = 2 # primitive element # verify a is a primitive element mod p s = {a**i % p for i in range(p)} assert( len(s) == p-1 ) n = p-2 dots = [] for i in range(1, n+1): for j in range(1, n+1): if (pow(a, i, p) + pow(a, j, p)) % p == 1: dots.append((i,j)) break

[1] Strictly speaking, no *scalable* algorithm will enumerate all Costas arrays. You could enumerate all permutation matrices of order *n* and test whether each is a Costas array, but this requires generating and testing *n*! matrices and so is completely impractical for moderately large *n*.

[2] More generally, the Lempel algorithm can generate a solution of order *q*-2 where *q* is a prime power. The code above only works for primes, not prime powers. For prime powers, you have to work with finite fields of order *q* and the code would be more complicated.

[3] A primitive root for a finite field of order *q* is a generator of the multiplicative group of the field, i.e. an element *x* such that every non-zero element of the field is some power of *x*.

Denote the number of teams as 2*n*. You’d like each team to play in each round, so you need *n* locations for the games to be played.

Each game will choose 2 teams from the set of 2*n* teams, so the number of games will be *n*(2*n* – 1). And this means you’ll need 2*n* – 1 rounds. A tournament design will be a grid with *n* rows, one for each location, and 2*n*-1 columns, one for each round.

Let’s pause there for just a second and make sure this is possible. We can certainly assign *n*(2*n* – 1) games to a grid of size *n* by 2*n*-1. But can we do it in such a way that no team needs to be in two locations at the same time? It turns out we can.

Now we’d like to introduce one more requirement: we’d like to balance the games over the locations. By that we mean we’d like a team to play no more than twice in the same location. A team has to play at least once in each location, but not every team can play twice in each location. If every team played once in each location, we’d have *n*² games, and if every team played twice in each location we’d have 2*n*² games, but

*n*² < *n*(2*n* – 1) < 2*n*²

for *n* greater than 1. If *n* = 1, this is all trivial, because in that case we would have two teams. They simply play each other and that’s the tournament!

We can’t have each team play exactly the same number of times in each location, but can we have each team play either once or twice in each location? If so, we call that a **balanced tournament design**.

Now the question becomes for which values of *n* can we have a balanced tournament design involving 2*n* teams. This is called a BTD(*n*) design.

We said above that this is possible if *n* = 1: two teams just play each other somewhere. For *n* = 2, it is not possible to have a balanced tournament design: some team will have to play all three games in the same location.

It is possible to have a balanced tournament design for *n* = 3. Here’s a solution:

In fact this is *the* solution aside from relabeling the teams. That is, given any balanced tournament design involving six teams, there is a way to label the teams so that you have the design above. Said another way, there is one design in BTD(3) up to isomorphism.

There are balanced tournament designs for all positive *n* except for *n* = 2. And in general there are a lot of them. There are 47 designs in BTD(4) up to isomorphism [1].

[1] The CRC Handbook of Combinatorial Designs. Colbourn and Dinitz editors. CRC Press, 1996.

The post Balanced tournament designs first appeared on John D. Cook.]]>

Japan is divided into 8 regions and 47 prefectures. Here is a network diagram of the prefectures in the Kanto region showing which regions border each other. (In this post, “border” will be regions share a large enough border that I was able to see the border region on the map I was using. Some regions may share a very small border that I left out.)

This is a good example of why it is convenient in GraphViz to use variable names that are different from labels. I created my graphs using English versions of prefecture names, and checked my work using the English names. Then after debugging my work I changed the label names (but not the connectivity data) to use Japanese names.

To show what this looks like, my GraphViz started out like this

graph G { layout=sfdp AI [label="Aichi"] AK [label="Akita"] AO [label="Aomori"] ... AO -- AK AO -- IW AK -- IW ...

and ended up like this

graph G { layout=sfdp AI [label="愛知県"] AK [label="秋田県"] AO [label="青森県"] ... AO -- AK AO -- IW AK -- IW ...

Here’s a graph only showing which prefectures border each other within a region.

This image is an SVG, so you can rescale it without losing any resolution. Here’s the same image as a PDF.

Because this network is effectively several small networks, it would be easy to look at a map and figure out which nodes correspond to which prefectures. (It would be even easier if you could read the labels!)

Note that there are two islands—literal islands, as well as figurative islands in the image above—Hokkaido, which is its own region, and Okinawa, which a prefecture in the Kyushu region.

Here’s the graph with all bordering relations, including across regions.

The image above is also an SVG. And here’s the same image as a PDF.

The post Graphing Japanese Prefectures first appeared on John D. Cook.]]>It’s an interesting exercise to recover the geographic regions from the network. For example, take a look at the graph for the continental United States.

It’s easy to identify Alaska in the graph. The node on the left represents Maine because Maine is the only state to border exactly one other state. From there you can bootstrap your way to identifying the rest of the states.

This could make a fun classroom exercise in a math class. Students will naturally come up with the idea of the degree of a node, the number of edges that meet that node, because that’s a handy way to solve the puzzle: the only possibilities for a node of degree *n* are states that border *n* other states.

This also illustrates that networks preserve topology, not geometry. That is, the connectivity information is retained, but the shape is dramatically different.

Someone asked me on Twitter to make a corresponding graph for Brazil. Mathematica, or at least my version of Mathematica, doesn’t have data on Brazilian states, so I made an adjacency graph using GraphViz.

Labeling the blank nodes is much easier for Brazil than for the US because Brazil has about half as many states, and the topology of the graph gives you more to work with. Three nodes connect to only one other node, for example.

Here the exercise doesn’t involve as much logic, but the geography is less familiar, unless of course you’re more familiar with Brazil than the US. Labeling the graph will require staring at a map of Brazil and you might accidentally learn a little about Brazil.

The labeled version of the graph above is available here. And here are the GraphViz source files that make the labeled and unlabeled versions.

The layout of a GraphViz file is very simple. The file looks like this:

graph G { layout=sfdp AC [label="Acre"] AL [label="Alagoas"] ... AC -- AM AC -- RO ... }

There are three parts: a layout, node labels, and connections.

GraphViz has several layout engines, and the `sfdp`

one matched what I was looking for in this case. Other layout options lead to overlapping edges that were confusing.

The node names AC, AL, etc. do not appear in the output. They’re just variable names for your convenience. The text inside the label is what appears in the final output. I’ll give an example in the next post in which it’s very convenient for the variables to be different from the labels. The order of the labels doesn’t matter, only which variables are associated with which labels.

Finally, the lines with variables separated by dashes are the connection data. Here we’re telling GraphViz to connect node AC to nodes AM and RO. The order of these lines doesn’t matter.

You can reduce the problem to coloring the nodes in a graph. Each node corresponds to a region, and there is an edge between two nodes if and only if their corresponding regions share a border.

Here is a sort of topologists’s or graph theorist’s view of the continental United States.

This was created using the following sample code from the Mathematica documentation.

RelationGraph[MemberQ[#2["BorderingStates"], #1] &, EntityList[ EntityClass["AdministrativeDivision", "ContinentalUSStates"]]]

You can recognize Maine in the graph because it’s the only state that only borders one other state. Alaska is also easy to locate. Exercise for the reader: mentally add Hawaii to the graph.

The analogous graph for Texas counties took much longer to draw: there are 49 continental US states but 254 Texas counties.

This was created with the following code.

RelationGraph[MemberQ[#2["BorderingCounties"], #1] &, EntityList[EntityClass["AdministrativeDivision", "USCountiesTexas"]]]

You can find El Paso county in the top left; it only borders one county just as Maine only borders one state.

The first was a tour of Africa. Actually two tours, one for the continent and one for islands. See this post for the Mathematica code used to create the tours.

The second was about the Americas: one tour for the North American continent, one for islands, and one for South America.

This post will look at Eurasia and Oceania. As before, I limit the tours to sovereign states, though there are disputes over which regions are independent nations. I first tried to do separate tours of Europe and Asia, but this would require arbitrarily categorizing some countries as European or Asian. The distinction between Asia and Oceania is a little fuzzy too, but not as complicated.

Here’s a map of the tour of Oceania.

Here’s the order of the tour:

- Australia
- East Timor
- Indonesia
- Palau
- Papua New Guinea
- Micronesia
- Marshall Islands
- Nauru
- Solomon Islands
- Vanuatu
- Fiji
- Tuvalu
- Kiribati
- Samoa
- Tonga
- New Zealand

The total length of the tour is 28,528 kilometers or 17,727 miles.

Here’s a map of the the Eurasian tour.

Here’s the order of the tour:

- Iceland
- Norway
- Sweden
- Finland
- Estonia
- Latvia
- Lithuania
- Belarus
- Poland
- Czech Republic
- Slovakia
- Hungary
- Romania
- Moldova
- Ukraine
- Georgia
- Armenia
- Azerbaijan
- Turkmenistan
- Uzbekistan
- Afghanistan
- Pakistan
- Tajikistan
- Kyrgyzstan
- Kazakhstan
- Russia
- Mongolia
- China
- North Korea
- South Korea
- Japan
- Taiwan
- Philippines
- East Timor
- Indonesia
- Brunei
- Malaysia
- Singapore
- Cambodia
- Vietnam
- Laos
- Thailand
- Myanmar
- Bangladesh
- Bhutan
- Nepal
- India
- Sri Lanka
- Maldives
- Yemen
- Oman
- United Arab Emirates
- Qatar
- Bahrain
- Saudi Arabia
- Kuwait
- Iran
- Iraq
- Syria
- Lebanon
- Jordan
- Israel
- Cyprus
- Turkey
- Bulgaria
- North Macedonia
- Serbia
- Bosnia and Herzegovina
- Montenegro
- Albania
- Greece
- Malta
- Italy
- San Marino
- Croatia
- Slovenia
- Austria
- Liechtenstein
- Switzerland
- Monaco
- Andorra
- Spain
- Portugal
- France
- Belgium
- Luxembourg
- Germany
- Netherlands
- Denmark
- United Kingdom
- Algeria

The total length of the tour is 61,783 kilometers or 38,390 miles.

The post Shortest tours of Eurasia and Oceania first appeared on John D. Cook.]]>Here’s the North American continental tour.

The order of the tour is as follows.

- Canada
- United States
- Mexico
- Guatemala
- El Salvador
- Costa Rica
- Panama
- Nicaragua
- Honduras
- Belize

Here’s a tour of the North American islands.

Trinidad and Tabago is about ten miles from the South American continent, but the country is classified as being part of North America, at least for some purposes.

Here is the order of the tour.

- Cuba
- Jamaica
- Haiti
- Dominican Republic
- Grenada
- Trinidad and Tobago
- Barbados
- Saint Vincent and the Grenadines
- Saint Lucia
- Dominica
- Antigua and Barbuda
- Saint Kitts and Nevis
- Bahamas

Here’s the tour of South America.

Here’s the order of the tour:

- Venezuela
- Guyana
- Suriname
- French Guiana
- Brazil
- Paraguay
- Uruguay
- Falkland Islands
- Argentina
- Chile
- Bolivia
- Peru
- Ecuador
- Colombia

Here’s my first attempt at a solution using Mathematica, based on an example in the documentation for `FindShortestTour`

.

africa = CountryData["Africa"] FindShortestTour[africa] GeoGraphics[{Thick, Red, GeoPath[africa[[%[[2]]]]]}]

This produced the following map:

Hmm. Maybe I should have been more specific about what I mean by “Africa.” My intention was to find a tour of *continental* Africa, i.e. not including islands. This means I needed to remove several items from Mathematica’s list of African countries. Also, I had in mind sovereign states, not territories of overseas states and not disputed territories.

After doing this, the map is more like what I’d expect.

The tour is then

- Algeria
- Tunisia
- Libya
- Egypt
- Chad
- Central African Republic
- Democratic Republic of the Congo
- Burundi
- Rwanda
- Uganda
- South Sudan
- Sudan
- Eritrea
- Djibouti
- Somalia
- Ethiopia
- Kenya
- Tanzania
- Malawi
- Zambia
- Mozambique
- Zimbabwe
- Eswatini
- Lesotho
- South Africa
- Botswana
- Namibia
- Angola
- Republic of the Congo
- Gabon
- Equatorial Guinea
- Cameroon
- Nigeria
- Niger
- Mali
- Burkina Faso
- Benin
- Togo
- Ghana
- Ivory Coast
- Liberia
- Sierra Leone
- Guinea
- Guinea-Bissau
- Gambia
- Senegal
- Mauritania
- Morocco

The initial tour, including islands, foreign territories, and Western Sahara, was 23,744 miles or 38,213 kilometers. The second tour was 21,074 miles or 33915 kilometers.

Here’s a tour of just the islands, excluding foreign territories.

The order of the tour is

- Cape Verde
- Seychelles
- Mauritius
- Madagascar
- Comoros
- São Tomé and Príncipe

This tour is 13,034 miles or 20,976 kilometers.

**Update**: See the next two posts for tours of the Americas and Eurasia and Oceania.

Triangles on a sphere or on a hyperbolic space like a pseudosphere have two kinds of angles: the sides come together at an angle, but the sides themselves are angles. By convention, the former are denoted with upper case letters and the latter with lower case letters.

The translation rule from spherical to hyperbolic geometry is to change functions of sides from circular to hyperbolic, but to leave functions of intersection angles alone. Or in typographical terms, put an *h* on the end of functions of a lower case letter but not functions of a upper case letter.

You can find these formulas, for example, in [1].

The Pythagorean theorem on a sphere

cos(*c*) = cos(*a*) cos(*b*).

becomes

cosh(*c*) = cosh(*a*) cosh(*b*)

in hyperbolic geometry. Here you simply change cos to cosh.

The law of sines on a sphere

sin(*a*) / sin(*A*) = sin(*b*) / sin(*B*) = sin(*c*) / sin(*C*).

becomes

sinh(*a*) / sin(*A*) = sinh(*b*) / sin(*B*) = sinh(*c*) / sin(*C*).

in hyperbolic geometry. Here sin becomes sinh, but only before a lower case letter, i.e. when applied to a side.

The law of cosines on a sphere

cos(*c*) = cos(*a*) cos(*b*) + sin(*a*) sin(*b*) cos(*C*).

becomes

cosh(*c*) = cosh(*a*) cosh(*b*) – sinh(*a*) sinh(*b*) cos(*C*).

Note the negative sign, a small exception to our conversion rule.

We could rewrite the law of cosines on a sphere to be

cos(*c*) = cos(*a*) cos(*b*) + κ sin(*a*) sin(*b*) cos(*C*)

where κ stands for curvature, which equals 1 for a unit sphere. Then our theorem translation rule holds exactly:

cosh(*c*) = cosh(*a*) cosh(*b*) + κ sinh(*a*) sinh(*b*) cos(*C*)

in a hyperbolic space with curvature κ = -1.

See this post on the Unified Pythagorean Theorem for a version of the Pythagorean theorem that holds in spherical, plane, and hyperbolic geometry.

Maybe there are analogous unified laws of sines and cosines. This is left as an exercise for the reader.

[1] William P. Thurston. Three-Dimensional Geometry and Topology, Volume 1. Princeton University Press, 1997.

The post Trig in hyperbolic geometry first appeared on John D. Cook.]]>arccos(*z*) / √(2 – 2*z*)

is analytic in a neighborhood of 1. We cancel out the bad behavior of inverse cosine at 1 by dividing by another function with matching bad behavior.

We can do something similar with Bessel functions. In general, the Bessel function *J*_{ν} has a branch cut from -∞ to 0. But the function

*J*_{ν} / *z*^{ν}

is entire, i.e. analytic everywhere in the complex plane. From a certain perspective—namely the perspective advocated by Bille Carlson (1924–2013)—we got the definition of Bessel functions (and other functions) wrong. Carlson’s book Special Functions of Applied Mathematics shows how many special functions can be defined in terms of a smaller number of functions with fewer branch cuts. This simplification revealed symmetries that hadn’t been noticed before.

Here’s a plot of *J*_{1/3}(*z*). Height represents magnitude and color represents phase.

This was produced using the Mathematica command

ComplexPlot3D[BesselJ[1/3, z], {z, -1 - I, 1 + I}]

The abrupt change from purple to yellow along the negative real axis shows the discontinuity in the phase at the branch cut. You can also see a hint of the cusp at the origin.

Here’s a plot of *J*_{1/3}(*z*) / *z*^{1/3}.

This was made using

ComplexPlot3D[BesselJ[1/3, z]/z^(1/3), {z, -1 - I, 1 + I}]

The white streak is an artifact of the branch cuts for *J*_{1/3}(*z*) and for *z*^{1/3}. The branch cut in their ratio is unnecessary.

I’ve verified that the following branch cuts are used by Mathematica, Common Lisp, and SciPy. If you know of any software that follows other conventions, please let me know in a comment.

The conventional branch cuts are as follows.

- sqrt: [-∞, 0)
- log: [-∞, 0]
- arcsin: [-∞, -1] and [1, ∞]
- arccos: [-∞, -1] and [1, ∞]
- arctan: [-∞
*i*, –*i*] and [*i*, ∞*i*] - arcsinh: [-∞
*i*, –*i*] and [*i*, ∞*i*] - arccosh: [-∞, 1]
- arctanh: [-∞, -1] and [1, ∞]

- Branch points in Common Lisp
- Numbering Lambert W function branches
- Trig functions in various programming languages

[1] W. Kahan. Branch Cuts for Complex Elementary Functions or Much Ado About Nothing’s Sign Bit. *The State of the Art in Numerical Analysis*. Clarendon Preess (1987).

You can find this series, for example, here.

This comes in handy, for example, when working with the analog of the Pythagorean theorem on a sphere.

You could just use the series and be on your way. But there’s a lot more going on than immediately meets the eye.

We’re looking at a series for inverse cosine centered at 1, and yet inverse cosine is multivalued in a neighborhood of 1: for an argument slightly less than 1, there are two possible angles that have such a cosine. That matters because in order to have a power series at a point, a function has to be well behaved in a disk around that point in the complex plane. Not only is our function not well behaved, it’s not even well *defined* until we consider branch cuts. Also, the function arccos(*z*) isn’t differentiable at 1; the derivative has a singularity at 1.

In a nutshell, we’re trying to expand a function in a series at a point where the function is badly behaved. That sounds impossible, or at least ill advised.

The series above is not a series for arccos per se. If we divide both sides of the series by √(2-2*z*) we see that we actually have a series for

Although arccos is badly behaved at 1, so is √(2-2*z*), and their ratio is well behaved at 1. In fact, it’s an analytic function and so it has a power series.

Inverse cosine has a **series expansion** but it does not have a **power series**. The series at the top is not a power series because it does not consist solely of powers of *z*. There’s a square root function on the right side as well, and this function is crucial to making things work.

Why would anyone think to find a series this way, i.e. why divide by √(2-2*z*)?

For small values of *z*,

and so

One might hope, correctly it turns out, that by dividing arccos(*z*) by a good approximation might yield a function nice enough to expand in a power series.

Let’s plot

to see whether it looks like a reasonable function. There’s usually no middle ground in complex variables: functions are either analytic or badly behaved. Both the numerator and denominator have branch cuts, but hopefully the cuts coincide and the ratio can be extended smoothly across the cuts.

The plot below suggests that is the case.

This plot was produced with

f[z_] := ArcCos[z] / Sqrt[2 - 2 z] ComplexPlot3D[f[z], {z, 0 - I, 2 + I}]

The white streak across the plot is not an accidental artifact of plotting but illustrates something important.

It is not possible to extend arccos(*z*) to a function that is analytic for all *z*. You have to exclude some values of *z* from the domain, i.e. you have to make branch cuts, and Mathematica makes these cuts along the real axis for *z* ≤ -1 and *z* ≥ 1.

The square root function also requires a branch cut, and Mathematica chooses that branch cut to be along the negative real axis, with means we have to exclude *z* ≥ 1. So the branch cuts of our numerator and denominator do coincide. (Inverse cosine has an additional branch cut, but it’s not near 1 so it doesn’t matter for our purposes.)

In summary, the series at the top of the post expands arccos at a point where the function is badly behaved, by dividing it by another function that is badly defined in the same way, making a function that is well behaved.

I rewrote the example in this post in using org-mode. My org file looked like this:

#+begin_src python :session :exports none lax_lat = 33.94 lax_lon = -118.41 iah_lat = 29.98 iah_lon = -95.34 #+end_src Suppose we want to find how far a plane would travel between Los Angeles (LAX) and Houston (IAH), ... #+begin_src python :session :exports none a = round(90 - lax_lat, 2) b = round(90 - iah_lat, 2) #+end_src /a/ = src_python[:session :results raw]{f"90° – {lax_lat}° = {a}°"} and /b/ = src_python[:session :results raw]{f"90° – {iah_lat}° = {b}°"} ...

Here are some observations about the experience.

First of all, writing the post in org-mode was **more work** than writing it directly, pasting in computed values by hand, but presumably **less error-prone**. It would also be easier to update. If, for example, I realized that I had the wrong coordinates for one of the airports, I could update the coordinates in one location and everything else would be updated when I regenerated the page.

I don’t think this was the best application of org-mode. It’s easier to use org-mode like a notebook, in which case you’re not trying to hide the fact that you’re mixing code and prose. I wanted to insert computed values into the text without calling attention to the fact that the values were computed. This is fine when you mostly have a text document and you only want to insert a few computed values. When you’re doing more computing it becomes tedious to repeatedly write

src_python[:session :results raw]{...}

to insert values. It might have been easier in this case to simply write a Python program that printed out the HTML source of the example.

There are **a couple advantages** to org-mode that weren’t relevant here. One is that the same org-mode file can be exported to **multiple formats**: HTML, LaTeX, ODT, etc. Here, however, I was only interested in exporting to HTML.

Another advantage of org-mode is the ability to **mix multiple programming languages**. Here I was using Python for everything, but org-mode will let you mix dozens of languages. You could compute one result in R, another result in Haskell, pass both results as arguments into some Java code, etc. You could also include data tables and diagrams in your org-mode file with your prose and code.

In general, keeping code and documentation together reduces errors. Literate programming may be more or less work, depending on the problem, but it reduces certain kinds of errors.

The example above is sorta bargain-basement literate programming. The code being developed was very simple, and not of interest beyond the article it was used in. Literate programming really shines when used to develop complex code, as in the book Physically Based Rendering. (**Update**: The third edition of this book is available online.)

**When code requires a lot of explanation**, literate programming can be very helpful. I did a project in psychoacoustics with literate programming a few years ago that would have been hard to do any other way. The project required a lot of reverse engineering and research. A line of code could easily require a couple paragraphs of explanation. Literate programming made the code development much easier because we could develop the documentation and code together, explaining things in the order most suited to the **human reader**, not to the compiler.

This post captures the algorithm in Python code. See the earlier post for documentation.

import re def char_to_num(ch): "Assumes all characters are digits or capital letters." n = ord(ch) if n <= ord('9'): # digits return n - ord('0') if n < ord('I'): # A-I return n - ord('A') + 1 if n <= ord('R'): # J-R return n - ord('J') + 1 return n - ord('S') + 2 # S-Z def checksum(vin): w = [8, 7, 6, 5, 4, 3, 2, 10, 0, 9, 8, 7, 6, 5, 4, 3, 2] t = 0 for i, c in enumerate(vin): t += char_to_num(c)*w[i] t %= 11 check = 'X' if t == 10 else str(t) return check def validate(vin): vinregex = "^[A-HJ-NPR-Z0-9]{17}$" r = re.match(vinregex, vin) return r and checksum(vin) == vin[8]

This code assumes the VIN number is given as ASCII or Unicode text. In particular, digits come before letters, and the numeric values of letters increase with alphabetical order.

The code could seem circular: the input is the full VIN, including the checksum. But the checksum goes in the 9th position, which has weight 0. So the checksum doesn't contribute to its own calculation.

**Update** I added a regular expression to check that the VIN contains only valid characters.

I wish I knew

The root of two.

O charmed was he

To know root three.

So we now strive

To find root five.

The beginning of each stanza is a mnemonic for the number mentioned in the following line.

√ 2 = 1.414

√ 3 = 1.732

√ 5 = 2.236

I found this in Twenty Years Before the Blackboard by Michael Stueben. Steuben sites the sources as Dictionary of Mnemonics, London 1972. He doesn’t give any other information such as author or editor.

**Update**: Additional verse from the comments.

It rhymes with heaven

The root of seven.

**Update**: Here’s Python code to validate the mnemonics above.

from math import sqrt def validate(line, num): digits = [str(len(w)) for w in line.split()] x = int("".join(digits)) / 1000 assert(x == round(sqrt(num),3)) validate("I wish I knew", 2) validate("O charmed was he", 3) validate("So we now strive", 5) validate("It rhymes with heaven", 7)

**Related post**: Numbers worth memorizing.

This page shows how to compute the expected number of runs in a sequence of binary outcomes when the two outcomes have probabilities *p* and 1-*p*. The number we’re after is one less than this, because the first run does not requiring changing input modes. The expected number of input mode changes is

μ = 2(*n* – 1)*p*(1 – *p*).

In our case we can take *p* = 26/36. It doesn’t matter whether *p* represents the probability of a letter or a digit because μ is symmetric in *p* and 1-*p*.

For moderately large *n*, the number of mode changes is approximately proportional to *n*. So for a password of *n* characters we should expect approximately

2 × 26/36 × 10/36 *n* ≈ 0.4 *n*

mode changes.

The variance in the number of mode changes is

σ² = 2*pq*(2*n *– 3 – 2*pq *(3*n *– 5))

which for large *n* is approximately

4*npq*(1 – 3*pq*).

In our case this reduces to

4 × 26/36 × 10/36 × (1 – 3 × 26/36 × 10/36) *n* ≈ 0.32 *n*.

The standard deviation is the square root of variance, which is approximately

σ ≈ 0.57 √*n*

in our example.

So typing a password of *n* characters will often require around

0.4 *n* ± 0.57 √*n*

mode changes.

For example, typing a 20-character password would often require between 5 and 11 input mode changes.

Variations on this problem come up fairly often. For example, yesterday I needed to create a spreadsheet with a permutation of the numbers 1 through 26. I had a spreadsheet with a column containing a permutation of the numbers 1 through 100, and so I copied it and deleted the rows with numbers greater than 26. The math above shows that on average you’d expect to change from a run of rows to keep to a run of rows to delete about 38 times.

The post Runs of letters and digits in passwords first appeared on John D. Cook.]]>- Typesetting modal logic
- Modal and temporal logic for security
- Word problems, modal logic, and regular expressions
- Axioms for S1 through S5
- Naming and numbering modal logic systems
- Modal logic and science fiction
- Dual axioms
- Temporal and polymodal logics

I created the image above by asking DALL-E 2 for an image of diamonds and boxes, a sort of pun on the ◇ (“diamond”) and □ (“box”) operators in modal logic.

The post Modal logic posts first appeared on John D. Cook.]]>This post will look at how circumference and area vary as a function of radius in three spaces: a surface with constant curvature 1 (i.e. a unit sphere), a surface of constant curvature 0 (i.e. a plane), and a surface of constant curvature -1 (a hyperbolic surface). The radius in each case is the distance from the center to the circle as measured on the surface.

In the spherical case, a circle of radius *r* has circumference

*C*(*r*) = 2π sin(*r*)

and area

*A*(*r*) = 2π (1 – cos(*r*)).

The circumference formula is valid for 0 ≤ *r* ≤ π and the area formula is valid for 0 ≤ *r* ≤ π/2.

In the hyperbolic case, a circle of radius *r* has circumference

*C*(*r*) = 2π sinh(*r*)

and area

*A*(*r*) = 2π (cosh(*r*) – 1).

These formulas are valid for 0 ≤ *r* < ∞.

Let’s make a couple plots to illustrate the equations above. First, circumference as a function of radius.

The top subplot looks surprising at first. Can circumference decrease when radius increases? Yes, on a sphere. Imagine a circle around the north pole. As we pull that circle down from the pole, the circumference increases until the circle becomes the equator. Past that point, the circumference of the circle decreases as the circle descends further south.

Using the same range of *r* for all three subplots obscures the fact that the circumference of a circle in the hyperbolic case grows exponentially with the radius.

Next let’s look at area as a function of radius.

In each case, area grows approximately quadratically with radius. But again the growth in the hyperbolic case is exponential as *r* increases further than is possible in the spherical case.

One last note: In the plane, the ratio of area to circumference is proportional to *r*. In the hyperbolic case, the same ratio is asymptotically constant, independent of *r*.

Here’s a plot of C(*r*) / *r*.

And here’s a plot of *A*(*r*) / *r*².

cos(*c*) = cos(*a*) cos(*b*)

where sides *a* and *b* are measured in radians. If we’re on a sphere of radius *R* and we measure the sides in terms of arc length rather than in radians, the formula becomes

cos(*c*/*R*) = cos(*a*/*R*) cos(*b*/*R*)

because an of length *x* has angular measure *x*/*R*.

How does this relate to the more familiar Pythagorean theorem on the plane? If *a*, *b*, and *c* are small relative to *R*, then the plane Pythagorean theorem holds approximately:

*c*² ≈ *a*² + *b*².

In this post I’ll present a version of the Pythagorean theorem that holds *exactly* on the sphere and the plane, and on a pseudosphere (hyperbolic space) as well. This is the Unified Pythagorean Theorem [1].

A sphere of radius *R* has curvature 1/*R*². A plane has curvature 0. A hyperbolic plane can have curvature *k* for any negative value of *k*.

Let *A*(*r*) be the area of a circle of radius *r* as measured on a surface of curvature *k*. Here area and radius are measured intrinsic to the surface. Then the Unified Pythagorean Theorem says

*A*(*c*) = *A*(*a*) + *A*(*b*) – *k* *A*(*a*) *A*(*b*) / 2π.

If *k* = 0, the final term on the right drops out, and we’re left with the ordinary Pythagorean theorem with both sides of the equation multiplied by π.

On a sphere of radius *R*, the area of a circle of radius *r* is

*A*(*r*) = 2π*R*²(1 – cos(*r*/*R*)).

Note that for small *x*,

1 – cos(*x*) ≈ *x*²/2,

and so *A*(*r*) ≈ π*r*² when *R* ≫ *r*. (Notation explained here.)

When you substitute the above definition for *A* in the unified theorem and plug in *k* = 1/*R*² you get

cos(*c*/*R*) = cos(*a*/*R*) cos(*b*/*R*)

as before.

In a hyperbolic space of curvature *k* < 0, let *R* = 1/√(-*k*). Then the area of a circle of radius *r* is

*A*(*r*) = 2π*R²*(cosh(*r*/*R*) – 1)

As with the spherical case, this is approximately the plane area when *R* ≫ *r* because

cosh(*x*) – 1 ≈ *x*²/2

for small *x*. Substituting the definition of *A* for hyperbolic space into the Universal Pythagorean Theorem reduces to

cosh(*c*/*R*) = cosh(*a*/*R*) cosh(*b*/*R*),

which is the hyperbolic analog of the Pythagorean theorem. Note that this is the spherical Pythagorean theorem with cosines replaced with hyperbolic cosines.

[1] Michael P. Hitchman. Geometry with an Introduction to Cosmic Topology. Theorem 7.4.7. Available here.

The post Unified Pythagorean Theorem first appeared on John D. Cook.]]>

Unfortunately this will be the last episode of Eclectic Tech. This last episode had several interesting stories. In addition to the story above, the episode discussed synchronizing clocks by observing cosmic ray events, a new bioinspired metamaterial, and NASA’s Inspire project.

As before we denote the sides of a triangle with *a*, *b*, and *c* and denote the angle opposite a side by the capital version of that letter. So *A* is the angle opposite side *a* and so forth.

The law of sines on a sphere says

sin *a* / sin *A* = sin *b* / sin *B* = sin *c* / sin *C*.

Suppose for a moment that we’re looking at a relatively small triangle, one for which the length of the sides is small relative to the radius of the sphere. Then we can apply the small angle approximation

sin *x* ≈ *x*

to each side [1] and find

*a* / sin *A* ≈ *b* / sin *B* ≈ *c* / sin *C*.

In other words, the law of sines from plane trigonometry holds approximately in spherical trigonometry for small triangles. This is not surprising since we know we can ignore the curvature of the earth when working over relatively small areas, but it is reassuring.

In the previous post we found the air distance between LAX (Los Angeles) and IAH (Houston) airports by making a spherical triangle with vertices at LAX, IAH, and the north pole. We labeled the side from LAX to the north pole *a* and the side from IAH to the north pole *b*. We found the length of side *c* using the law of cosines.

Now suppose we’re at LAX and about to fly to IAH. What direction should we fly? Since Los Angeles and Houston are approximately at the same latitude, and they’re not terribly far apart relative to the size of the earth, we know that we’ll roughly fly due east. Let’s find the flight direction more precisely.

The angle *B* in our triangle gives us the angle between our heading and a heading of due north. We know *b* from the latitude of IAH, we know *C* from the difference in longitude of the two airports, and we found *c* in the previous post. Therefore we can find *B* from the law of sines:

sin *b* / sin *B* = sin *c* / sin *C*

and we know everything but *B*. So

sin *B* = sin *b* sin *C* / sin *c* = sin 60.02° sin 23.07° / sin 19.92°.

This tells us sin *B* = 0.9960. So what is *B*? If we ask software for the arcsine of 0.9960 we will get 84.88°. This would say we head 5.12° north of due east to get from Los Angeles to Houston. But Houston is at a lower latitude than Los Angeles and so that doesn’t seem right.

In the first post in this series I mentioned that a flight between Houston and Wuhan would fly over Alaska even though the two cities are at essentially the same latitude. Wuhan is at a slightly higher latitude than Houston, but a flight leaving Wuhan for Houston would head somewhat north. So it is *possible* that a flight from a higher latitude to a lower latitude would head north [2].

But Los Angeles is much closer to Houston than Wuhan is, and the difference in latitude is greater. A flight from LA to Houston would head mostly east and a little south. So where did we go wrong?

The arcsine function in any software package returns **a** solution to the equation

sin(*x*) = *y*,

but it might not return the solution we need. The angle 84.88° is a valid solution to the equation

sin(*B*) = 0.9960

but it’s not the appropriate solution in our context. We need the solution 95.12°. Our flight heada 5.12° **south** of east. This assumes our flight takes a great circle path; airlines may have reasons not to take the shortest possible path between two airports.

***

[1] The approximation assumes *x* is measured in radians. If we measure angles in degrees or any other units, the right side of the approximation will need to be multiplied by a constant. But the same constant will appear in all the numerators and denominators of the law of sines and so the approximation holds in any units.

[2] Here’s another illustration of the same point. Draw a circle around the north pole, say a circle a mile in diameter, and suppose you want to walk from longitude 0 to longitude 180°. Following a path of constant latitude is walking around the circle. Heading north instead is walking across the diameter of the circle. If your destination is a little bit outside the circle on the other side, it’s still shorter to walk across the circle, i.e. head north, even though your destination is at a lower latitude, i.e. south.

The post Direction between two cities first appeared on John D. Cook.]]>