Tetrahedral analog of the Pythagorean theorem

Posted on by John

A tetrahedron has four triangular faces. Suppose three of those faces come together like the corner of a cube, each perpendicular to the other. Let A1, A2, and A3 be the areas of the three triangles that meet at this corner and let A0 be the area of remaining face, the one opposite the right corner.

De Gua’s theorem, published in 1783, says

A0² = A1² + A2² + A3².

We will illustrate De Gua’s theorem using Python.

from numpy import *

def area(p1, p2, p3):
    return 0.5*abs(linalg.norm(cross(p1 - p3, p2 - p3)))

p0 = array([ 0, 0,  0])
p1 = array([11, 0,  0])
p2 = array([ 0, 3,  0])
p3 = array([ 0, 0, 25])

A0 = area(p1, p2, p3)
A1 = area(p0, p2, p3)
A2 = area(p0, p1, p3)
A3 = area(p0, p1, p2)

print(A0**2 - A1**2 - A2**2 - A3**2)

This prints 0, as expected.

Higher dimensions

A natural question is whether De Gua’s theorem generalizes to higher dimensions, and indeed it does.

Suppose you have a 4-simplex where one corner fits in the corner of a hypercube. The 4-simplex has 5 vertices. If you leave out any vertex, the remaining 4 points determine a tetrahedron. Let V0 be the volume of the tetrahedron formed by leaving out the vertex at the corner of the hypercube, and let V1, V2, V3, and V4 be the volumes of the tetrahedra formed by dropping out each of the other vertices. Then

V0² = V1² + V2² + V3² + V4².

You can extend the theorem to even higher dimensions analogously.

Let’s illustrate the theorem for the 4-simplex in Python. The volume of a tetrahedron can be computed as

V = det(G)1/2/6

where G is the Gram matrix computed in the code below.

def volume(p1, p2, p3, p4):
    u = p1 - p4
    v = p2 - p4
    w = p3 - p4

    # construct the Gram matrix
    G = array( [
            [dot(u, u), dot(u, v), dot(u, w)],
            [dot(v, u), dot(v, v), dot(v, w)],
            [dot(w, u), dot(w, v), dot(w, w)] ])

    return sqrt(linalg.det(G))/6

p0 = array([ 0, 0,  0, 0])
p1 = array([11, 0,  0, 0])
p2 = array([ 0, 3,  0, 0])
p3 = array([ 0, 0, 25, 0])
p4 = array([ 0, 0,  0, 7])

V0 = volume(p1, p2, p3, p4)
V1 = volume(p0, p2, p3, p4)
V2 = volume(p0, p1, p3, p4)
V3 = volume(p0, p1, p2, p4)
V4 = volume(p0, p1, p2, p3)

print(V0**2 - V1**2 - V2**2 - V3**2 - V4**2)

This prints -9.458744898438454e-11. The result is 0, modulo the limitations of floating point arithmetic.

Numerical analysis

Floating point arithmetic is generally accurate to 15 decimal places. Why is the numerical error above relatively large? The loss of precision is the usual suspect: subtracting nearly equal numbers. We have

V0 = 130978.7777777778

and

V1**2 + V2**2 + V3**2 + V4**2 = 130978.7777777779

Both results are correct to 16 decimal places, but when we subtract them we lose all precision.

The anti-Smith chart

Posted on by John

As I’ve written about several times lately, the Smith chart is the image of a rectangular grid in the right half-plane under the function

f(z) = (z − 1)/(z + 1).

What would the image of a grid in the left half-plane look like?

For starters, since f maps the right half-plane to the interior of the unit circle, it must map the left-half plane to the exterior of the unit circle.

As we said before, the function f is a Möbius transformation, and so it takes generalized circles, i.e. either a circle or a line, to generalized circles. So the grid lines in the left half-plane are either mapped to lines or circles. Which is it?

The function f has a singularity at −1 and so the image of any line (or circle) through z = −1 is unbounded, i.e. a line, not a circle. Any line not passing through −1 has a bounded image, which must be a circle.

The line Re(z) = −1 in the z plane is mapped to the line Re(w) = 1 in the w plane. Otherwise a vertical line crossing the real axis at x is mapped to a circle passing through w = (x − 1)/(x + 1). The circle also passes through w = 1 because f(∞) = 1. The circle is symmetric about the real axis, and so this is enough information to uniquely determine the circle.

Note that (x − 1)/(x + 1) > 1 when x < −1 and so vertical lines with real part less than −1 are mapped to circles to the right of w = 1. When −1 < x < 0, vertical lines are mapped to circles to the left of w = 1.

The images of horizontal lines we’ve looked at before. These are all circles passing through w = 1 and tangent to the circles that are images of vertical lines. But this time instead of taking the portion of the circles inside the unit circle, we take the portion outside the unit circle.

And without further ado, we present the anti-Smith chart, the image of a grid in the left half plane.

 

Impedance and Triangular Numbers

Posted on by John

A few days ago I wrote two posts about how to create a Smith chart, a graphical device used for impedance calculations. Then someone emailed me to point out the connection between the Smith chart and triangular numbers.

The Smith chart is the image of a rectangular grid in the right half-plane under the function

f(z) = (z − 1)/(z + 1).

If you subtract the values of f at consecutive integers, you get the reciprocal of a triangular number.

f(n) − f(n − 1) = 2/(n(n + 1)) = 1 / Tn

Or to put it another way,

f(n) − f(n − 1) = 1 / (1 + 2 + 3 + … + n).

In the first post on the Smith chart we showed that the function f maps vertical lines

in the z plane to circles in the w plane all touching at w = 1.

The circles are symmetric about the real axis and the diameter runs from f(n) to 1. The separation between the circles on the left side is thus

f(n) − f(n − 1) = 1 / Tn.

Number the circles starting from the outermost as 0, 1, 2, …. Then the maximum distance between circle n and circle n − 1 is 1 / Tn. You can see in the graph above that the distance between circle 0 and circle 1 is 1. It’s a little harder to see that the distance between circle 1 and circle 2 is 1/3. It looks like the distance between circles 2 and 3 is about half of that between circles 1 and 2, so it would be 1/6.

Related posts

Cross ratio

Posted on by John

The cross ratio of four points A, B, C, D is defined by

(A, B; C, D) = \frac{AC \cdot BD}{BC \cdot AD}

where XY denotes the length of the line segment from X to Y.

The idea of a cross ratio goes back at least as far as Pappus of Alexandria (c. 290 – c. 350 AD). Numerous theorems from geometry are stated in terms of the cross ratio. For example, the cross ratio of four points is unchanged under a projective transformation.

Complex numbers

The cross ratio of four (extended [1]) complex numbers is defined by

(z_1, z_2; z_3, z_4) = \frac{(z_3 - z_1)(z_4 - z_2)}{(z_3 - z_2)(z_4 - z_1)}

The absolute value of the complex cross ratio is the cross ratio of the four numbers as points in a plane.

The cross ratio is invariant under Möbius transformations, i.e. if T is any Möbius transformation, then

(T(z_1), T(z_2); T(z_3), T(z_4)) = (z_1, z_2; z_3, z_4)

This is connected to the invariance of the cross ratio in geometry: Möbius transformations are projective transformations on a complex projective line. (More on that here.)

If we fix the first three arguments but leave the last argument variable, then

T(z) = (z_1, z_2; z_3, z) = \frac{(z_3 - z_1)(z - z_2)}{(z_3 - z_2)(z - z_1)}

is the unique Möbius transformation mapping z1, z2, and z3 to ∞, 0, and 1 respectively.

The anharmonic group

Suppose (a, b; c, d) = λ ≠ 1. Then there are 4! = 24 permutations of the arguments and 6 corresponding cross ratios:

\lambda, \frac{1}{\lambda}, 1 - \lambda, \frac{1}{1 - \lambda}, \frac{\lambda - 1}{\lambda}, \frac{\lambda}{\lambda - 1}

Viewed as functions of λ, these six functions form a group, generated by

\begin{align*} f(\lambda) &= \frac{1}{\lambda} \\ g(\lambda) &= 1 - \lambda \end{align*}

This group is called the anharmonic group. Four numbers are said to be in harmonic relation if their cross ratio is 1, so the requirement that λ ≠ 1 says that the four numbers are anharmonic.

The six elements of the group can be written as

\begin{align*} f(\lambda) &= \frac{1}{\lambda} \\ g(\lambda) &= 1 - \lambda \\ f(f(\lambda)) &= g(g(\lambda) = z \\ f(g(\lambda)) &= \frac{1}{\lambda - 1} \\ g(f(\lambda)) &= \frac{\lambda - 1}{\lambda} \\ f(g(f(\lambda))) &= g(f(g(\lambda))) = \frac{\lambda}{\lambda - 1} \end{align*}

Hypergeometric transformations

When I was looking at the six possible cross ratios for permutations of the arguments, I thought about where I’d seen them before: the linear transformation formulas for hypergeometric functions. These are, for example, equations 15.3.3 through 15.3.9 in A&S. They relate the hypergeometric function F(a, b; c; z) to similar functions where the argument z is replaced with one of the elements of the anharmonic group.

I’ve written about these transformations before here. For example,

F(a, b; c; z) = (1-z)^{-a} F\left(a, c-b; c; \frac{z}{z-1} \right)

There are deep relationships between hypergeometric functions and projective geometry, so I assume there’s an elegant explanation for the similarity between the transformation formulas and the anharmonic group, though I can’t say right now what it is.

Related posts

[1] For completeness we need to include a point at infinity. If one of the z equals ∞ then the terms involving ∞ are dropped from the definition of the cross ratio.

Text case changes the size of QR codes

Posted on by John

Let’s make a QR code out of a sentence two ways: mixed case and upper case. We’ll use Python with the qrcode library.

>>> import qrcode
>>> s = "The quick brown fox jumps over the lazy dog."
>>> qrcode.make(s).save("mixed.png")
>>> qrcode.make(s.upper()).save("upper.png")

Here are the mixed case and upper case QR codes.

The QR code creation algorithm interprets the mixed case sentence as binary data but it interprets the upper case sentence as alphanumeric data.

Alphanumeric data, in the context of QR codes, comes from the following alphabet of 44 characters:

0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ $%*+-.:

Since 44² = 1936 < 2048 = 211 two alphanumeric characters can be encoded in 11 bits. If text contains a single character outside this alphabet, such as a lower case letter, then the text is encoded as ISO/IEC 8859-1 using 8 bits per character.

Switching from mixed-case text to upper case text reduces the bits per character from 8 to 5.5, and so we should expect the resulting QR code to require about 30% fewer pixels. In the example above we go from a 33 × 33 grid down to a 29 × 29 grid, from 1089 pixels to 841.

Application to Bitcoin addressess

Bech32 encoding uses an alphabet of 32 characters while Base58 encoding uses an alphabet of 58 characters, and so the former needs about 17% more characters to represent the same data. But Bech32 uses a monocase alphabet, and base 58 does not, and so Bech32 encoding requires fewer QR code pixels to represent the same data as Base58 encoding.

(Bech32 encoding uses a lower case alphabet, but the letters are converted to upper case before creating QR codes.)

Related posts

An ancient generalization of the Pythagorean theorem

Posted on by John

Apollonius of Perga (c. 262 BC – c. 190 BC) discovered a theorem that generalizes the Pythagorean theorem but isn’t nearly as well known.

Let ABC be a general triangle, and let D be the midpoint of the segment AB. Let a be the length of the side opposite A and b the length of the side opposite B. Let m be the length of AD and h the length of the mediant, the line CD.

Apollonius’s theorem says

a² + b² = 2(m² + h²).

To see that this is a generalization of the Pythagorean theorem, apply Apollonius’ theorem to an isosceles triangle. Now a = b and ACD is a right triangle.

Apollonius’ theorem says

2b² = 2m² + 2h²

which is the Pythagorean theorem applied to ACD with each term doubled.

Mentally compute logs base 2

Posted on by John

The previous post required computing

\frac{128}{\log_2 5}

After writing the post, I thought about how you would mentally approximate log2 5. The most crude approximation would round 5 down to 4 and use log2 4 = 2 to approximate log2 5. That would be good enough for an order of magnitude guess, but we can do much better without too much more work.

Simple approximation

I’ve written before about the approximation

\log_2 x \approx 3\frac{x - 1}{x + 1}

for x between 1/√2 and √2. We can write 5 as 4 (5/4) and so

\begin{align*} \log_2 5 &= \log_2 4 (5/4) \\ &= \log_2 4 + \log_2 (5/4) \\ &\approx 2 + 3\frac{5/4 - 1}{5/4 + 1} \\ &= 2 + 3 \frac{1/4}{9/4} \\ &= 7/3 \end{align*}

How accurate is this? The exact value of log2 5 is 2.3219…. Approximating this number by 7/3 is much better than approximating it by just 2, reducing the relative error from 16% down to 0.5%.

Origin story

Where did the approximation

\log_2 x \approx 3\frac{x - 1}{x + 1}

come from?

I don’t remember where I found it. I wouldn’t be surprised if it was from something Ron Doerfler wrote. But how might someone have derived it?

You’d like an approximation that works on the interval from 1/√2 to √2 because you can always multiply or divide by a power of 2 to reduce the problem to this interval. Rational approximations are the usual way to approximate functions over an interval [1], and for mental calculation you’d want to use the lowest order possible, i.e. degree 1 in the numerator and denominator.

Here’s how we could ask Mathematica to find a rational approximation for us [2].

Simplify[
    N[
        ResourceFunction["EconomizedRationalApproximation"][
            Log[2, x], { x, {1/Sqrt[2], Sqrt[2]}, 1, 1}]]]

This returns

(2.97035 x − 2.97155) / (1.04593 + x)

which we round off to

(3 x − 3) / (1 + x).

The N function turns a symbolic result into one with floating point numbers. Without this call we get a complicated expression involving square roots and logs of rational numbers.

The Simplify function returns an algebraically equivalent but simpler expression for its argument. In our case the function finishes the calculation by removing some parentheses.

Related posts

[1] Power series approximations are easier to compute, but power series approximations don’t give the best accuracy over an interval. Power series are excellent at the point where they’re centered, but degrade as you move away from the center. Rational approximations spread the error more uniformly.

[2] I first tried using Mathematica’s MiniMaxApproximation function, but it ran into numerical problems, so I switched to EconomizedRationalApproximation.

Physical Keys and Encryption Keys

Posted on by John

A physical key, such as a house key, is a piece of metal with cuts of differing depths. Typically there may be around 6 cuts, with five different possible depths for each cut. This allows 56 = 15,625 possible keys.

Encryption keys, such as AES keys, are a string of bits, often 128 bits, for a total of 2128 possible keys.

How long would a physical key have to be to have the same level of security as an encryption key? We’d need to solve

5n = 2128

which means

n = 128 / log25 = 55.12.

So we’d need a key with around 55 notches.

metal key with 55 notches

This only takes into account combinatorial possibilities, not the difficulty of attacking a physical key or a binary key. There are incomparably more possibilities for binary keys, but encryption attacks can be automated and carried out remotely (unless a computer is air gapped). A physical lock can only be attacked in person. It takes a lock picker orders of magnitude more time to try a key than a password cracking program. On the other hand, locks aren’t picked by trying thousands of keys.

Related post: Measuring cryptographic strength in liters of boiling water

Freshman’s dream

Posted on by John

The “Freshman’s dream” is the statement

(x + y)p = xp + yp

It’s not true in general, but it is true mod p if p is a prime. It’s a cute result, but it’s also useful in applications, such as finite field computations in cryptography.

Here’s a demonstration of the Freshman’s dream in Python.

>>> p = 5
>>> [((x + y)**p - x**p - y**p) % p for x in range(p) for y in range(p)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Here’s an example using a prime too large for verify the results by looking at the output.

>>> import numpy as np
>>> p = 103
>>> v = [((x + y)**p - x**p - y**p) % p for x in range(p) for y in range(p)]
>>> np.all( np.array(v) == 0 )
True

You can use the same code to show that the Freshman’s dream is not true in general if p is not a prime, and it’s not true in general if p is a prime but the exponent is less than p.

Why and how Bitcoin uses Merkle trees

Posted on by John

Yesterday’s post looked at a recently mined Bitcoin block, the 920,994th block in the blockchain, and verified that it contains the hash of the previous block. This post will look at the same block and verify its Merkle tree root.

Before getting down to the bytes, we’ll back up and say what a Merkle tree is and why Bitcoin uses one.

What is a Merkle tree?

A Merkle tree is a binary tree. The leaves of the tree are labeled with a hash of the data associated with each leaf. The label of each interior node is the hash of the concatenation of the two child hashes. The label of the root of the tree is the Merkle tree root.

The smallest change to any data in the tree will change the Merkle tree root [1]. If you weren’t concerned with efficiency, you wouldn’t need a Merkle tree. You could simply concatenate all the data in the leaves and compute a hash. But the Merkle tree makes it possible to verify a node without having to verify the whole tree.

To see this, imagine starting at the bottom of the tree and working your way up to compute the Merkle root. You need the hash value of the node you’re interested in and its sibling. Then you can compute the label of the parent node. If you know the label of the parent’s sibling node (you could call it the aunt or uncle) then you can compute the label of the grandparent node, and so forth until you get to the root. If you get the expected value for the root, the data is correct.

How Bitcoin uses Merkle trees

Bitcoin makes a Merkle tree with the data for each transaction being a leaf node. Simple Payment Verification (SPV) lets you verify a transaction knowing only the hash of the transaction, the hash of its sibling, and the hashes of the siblings moving up the tree: the aunts, great aunts, great great aunts etc.

If a block has 2n transactions, you only need n hash values to verify a single transaction. The binary tree has n levels, and you need two hashes at the bottom, and one hash for each of the levels between the bottom and the top. So in a block of 1024 transactions, you need 10 hash values. Note that this saves space two ways: you don’t need transaction data per se, only hashes of transaction data, and the number of hashes you need is log2 of the number of transactions.

What if the number of transactions is not a power of 2? If there is an odd number of nodes at a given level of the Merkle tree, the last node is paired with itself.

For example, the block that we’ll examine in detail contains 1279 transactions. At the bottom of our Merlke tree there are effectively 1280 nodes, with the 1280th node being a duplicate of the 1279th. Moving up the tree we have 640 nodes, 320 nodes, 160 nodes, … 10 nodes, and 5 nodes. The 5th node is paired with itself, and no the next level up has 3 nodes. The third node is paired with itself, making 2 nodes on the next level, and then the root above that.

Down to the bytes

The command

xxd -l 80 920994.dat

gives a hex dump of the block header:

0000 002e c2df b57e f582 75c2 7a64 33c5
edfd dfc1 af6a b9ae 708c 0000 0000 0000
0000 0000 bb15 accc 8940 3811 0b9b e1cf
b125 c0fa a65e fb1b 2fa4 61ed 989f 8d6f
e41e 04cf 5522 ff68 21eb 0117 f2bc ab1a

This time we’re looking at bytes 37 through 68, highlighted in orange. This is the Merkle tree root.

A miner verifying the block would recompute the Merkle tree root and confirm that it matches the value given in the header.

A Bitcoin block does not contain a Merkle tree per se, only its root. The hash of the data for each transaction, what Bitcoin calls a txid, is redundant since it can be computed. And if it were included in the block, a miner would need to recompute it to verify that it is correct.

We can fetch a list of txids for a block by calling

https://blockstream.info/api/block/{block_hash}/txids

with the hash of the block we want, which in our case is

00000000000000000001c2f89cee9f8c2fe88b4a93c2c2b75192918fc438ca32

This returns a JSON string with 1279 hash values:

fc1d784c62600565d4b4fbfe9cd45a7abc5f1c3273e294fd2b44450c5c9b9e9a
d4bce41744156a8f0e6588e42efef3519904a19169212c885930e908af1457ec
7bc6428c7d023108910a458c99359a2664de7c02343d96b133986825297e518d
…
c9132f178830d0c7a781e246bb5c2b3f4a9686c3d2b6f046b892a073d340eaff

Note that the API call did not pull the txids from the block; it computed them. Presumably the server computed the list of txids once and cached the results rather than computing them when we called the API.

We can compute the Merkle tree root by concatenating pairs of txids and hashing the results. Since the number of transactions is odd, the last txid is paired with itself. So moving one level up the tree, the first node is labeled with the hash of

fc1…e9ad4b…7ec

and the last node is labeled with the hash of

c91…affc91…aff

We continue this process until we compute the Merkle tree root. Then we can confirm that it matches the value given in the header. As with the block hash in the previous post, the Merkle tree root is stored in the header in little endian form.

Related posts

[1] There is a vanishingly small chance that a change in the data would not change the root. With a 256-bit hash, the chances of a change in the data not changing the root are 1 in 2256.