MoreRSS

site iconJohn D. CookModify

I have decades of consulting experience helping companies solve complex problems involving applied math, statistics, and data privacy.
Please copy the RSS to your reader, or quickly subscribe to:

Inoreader Feedly Follow Feedbin Local Reader

Rss preview of Blog of John D. Cook

Computing polar triangles with linear algebra

2026-10-25 09:01:55

The previous post defined the polar triangle of a spherical triangle in terms of its sides and angles, but it doesn’t construct the triangle. It only specifies one of a family of triangles located anywhere on the sphere.

How do we know that polar triangles exist? There is a geometric construction that starts with a spherical triangle and constructs a unique new spherical triangle, and one can show that this new triangle has the angles and sides described in the previous post.

There is also a linear algebraic construction of polar triangles, and that’s the topic of this post. In a nutshell, take your sphere to be the unit sphere in 3D. Write the coordinates of each vertex as the columns of a matrix. Invert that matrix and normalize each row. Then the rows are the coordinates of the polar vertices.

This process is an involution: the polar of the polar of a triangle is the original triangle.

The following Python code illustrates the process described above. It begins with a triangle on the globe with vertices A, B, and C at Chicago, Quito, and Paris respectively. It computes the polar triangle, then computes the polar of that, and shows that the result is the original triangle.

In case you’re curious,

  • A′ is at 40.77° N, 168.33° W
  • B′ is at 34.53° S, 35.70° W
  • C′ is at 10.00° N, 11.43° E
import numpy as np

def latlon_to_cartesian(lat_deg, lon_deg):
    lat = np.radians(lat_deg)
    lon = np.radians(lon_deg)
    x = np.cos(lat) * np.cos(lon)
    y = np.cos(lat) * np.sin(lon)
    z = np.sin(lat)
    return np.array([x, y, z])

def cartesian_to_latlon(v):
    v = v / np.linalg.norm(v)
    lat = np.degrees(np.arcsin(np.clip(v[2], -1, 1)))
    lon = np.degrees(np.arctan2(v[1], v[0]))
    return lat, lon

def polar_triangle(M):
    """
    Polar triangle via matrix inversion.
    M: 3×3 matrix with columns = unit vectors of the original vertices.
    Returns 3×3 matrix whose columns are the polar unit vectors.
    """
    Minv = np.linalg.inv(M)
    # Normalize each row of the inverse (rows are parallel to the polar directions)
    polar_rows = Minv / np.linalg.norm(Minv, axis=1, keepdims=True)
    return polar_rows.T          # columns = polar vertices

# Original vertices (Chicago, Quito, Paris)
A = latlon_to_cartesian(41.85003, -87.65005)
B = latlon_to_cartesian(-0.22985, -78.52495)
C = latlon_to_cartesian(48.85341,  2.3488)

M = np.column_stack([A, B, C])
print("M (columns = original vertices):\n", M)
print("det(M) =", np.linalg.det(M))

# Polar triangle
M_prime = polar_triangle(M)
print("\nM' (columns = polar vertices):\n", M_prime)

print("\nPolar geographic coordinates (lat, lon):")
for name, col in zip(["A'", "B'", "C'"], range(3)):
    print(name, cartesian_to_latlon(M_prime[:, col]))

# Polar of the polar
M_double = polar_triangle(M_prime)
print("\nPolar-of-polar matrix:\n", M_double)
print("max |difference from original M| =",
      np.max(np.abs(M_double - M)))
The post Computing polar triangles with linear algebra first appeared on John D. Cook.

Duality in spherical trigonometry

2026-10-25 08:28:20

I’ve written before about the law of cosines on a sphere:

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

But that’s just one of two laws of cosines on a sphere. The equation above is the law of cosines for sides. The law of cosines for angles is typographically very similar, just inverting the case of the letters used as arguments.

cos(C) = cos(A) cos(B) + sin(A) sin(B) cos(c).

Incidentally, the law of cosines for sides shows how the sides of a spherical triangle determine the angles. The law of cosines for angles shows how, unlike in the case of plane geometry, angles determine the sides.

Two kinds of angles

Let’s back up and say what the letters refer to. As is conventional in plane geometry, the sides of a triangle are a, b, and c, and the angles are A, B, and C, with a being the side opposite the angle A, and similarly for the other letters.

However, unlike plane geometry, the sides are also angles, but of a different sort. The sides a, b, and c are central angles, and the “angles” A, B, and C are dihedral angles.

The sides are arcs of a great circle. The side a, for example, is a portion of a great circle, running between the vertices at B and C [1]. Think of a plane that contains the center of the sphere O and the vertices B and C. The length of a is the angle at O.

The angles A, B, and C are the angles at which great circles meet. In the formalism of differential geometry, these angles exist in the tangent plane at the vertices. These angles are intrinsic, measured along the (tangent plane of) the surface of the sphere. The central angles are extrinsic, measured in a plane that cuts through the sphere to its center.

Polar triangles

Now because sides are angles and vertices are angles, we can interchange the two kinds of angles to form the polar triangle of the original triangle. This is a new triangle whose central angles are the complements of the original dihedral angles, and whose dihedral angles are the complements of the original central angles. In symbols,

a′ = π − A
b′ = π − B
c′ = π − C
A′ = π − a
B′ = π − b
C′ = π − c

Now we return to the laws of cosines at the top of the post: the law of cosines for angles is the law of cosines for sides applied to the polar triangle. Every equation in spherical trigonometry has a polar counterpart.

 

[1] There’s some overloaded terminology here, just as in plane geometry. The symbol a is used to denote a side of a triangle and also the length of that side. The symbol A is used to denote both a vertex and the dihedral angle at that vertex.

The post Duality in spherical trigonometry first appeared on John D. Cook.

Working around a strange Python problem

2026-10-22 23:12:10

A Python script that worked a few days ago wouldn’t work today. I suspect it might have something to do with my updating everything I could think of in response to recent security breaches.

I got an error message that said

error: The script contains multiple PEP 723 metadata blocks

I searched on what that means, then took out the blank line in my header

#!/usr/bin/env -S uv run --script

# /// script
# dependencies = ["pyperclip"]
# ///

but I still got the same error message.

The script in question copies miscellaneous snippets of text to the clipboard for me. One of those snippets is the header block that uv requires. The script contains the header block marker # ///, and that was the source of the problem.

According to Grok

The PEP 723 parser is purely textual (it uses a regex that scans the whole file). It does not understand Python syntax, so the # /// sequences inside your triple-quoted string are treated as a second metadata block.

It suggested changing my quoted string to use # / / / and then take out the extra spaces:

"""#!/usr/bin/env -S uv run --script
# / / / script
# dependencies = ["numpy", "matplotlib"]
# / / /
import numpy as np
import matplotlib.pyplot as plt
""".replace("# / / /", "# ///")

This is a dirty hack, but it works.

The post Working around a strange Python problem first appeared on John D. Cook.

How the eMudrha certificate debacle might have happened

2026-10-22 22:01:55

You may have seen references to a password being the number 4. There are countless allusions to this online, but they don’t give the context.

As I mentioned in an earlier post, eMudrha was using the numbers 4, 56, 115, 155, 156, 230, and 237 as secret keys. I haven’t seen details of how this happened, but in this post I’ll give a plausible account of how it might have happened, what someone could have done to produce those results. I’m not trying to do forensic analysis on this specific incident, just exploring what someone might do.

eMudrha was using the NIST elliptic curve P-384, and a secret key should be a random number chosen between 1 and 2384. That upper bound isn’t exactly correct, but it’s close. See this post if you want to know the exact number.

Someone might have chosen numbers up to 384 rather than 2384 out of ignorance. Or they might done so by a typo. Maybe they generated a random value that might have been larger than 2384 and then intended to reduce the number mod 2384, but instead of typing

% 2**384

they typed

% 384

It seems more likely, however, that somehow someone took one byte of a random value because all the keys are less than 28 = 256.

The probability of choosing seven numbers from 1 to 384 and all being less than or equal to 237 is around 0.034. But the same probability for numbers chosen between 1 and 256 is about 0.60.

The post How the eMudrha certificate debacle might have happened first appeared on John D. Cook.

The number of points on an elliptic curve

2026-10-22 19:08:51

An elliptic curve E over a finite field F consists of the set of points (x, y) satisfying the equation

y² = x³ + ax + b

with x and y in F. Typically in cryptography the field F is the integers mod p where p is a very large prime number.

The number of points on E is roughly the number of elements in F. So for a curve over the integers mod p, the curve has roughly p points.

The next section will give an intuitive argument for why the number of points on E should be roughly the number of elements in F. Then the next section will give a rigorous bound. Then after that we’ll work through an example.

Hand-wavy estimate

For a prime p, half the non-zero integers mod p have a square root and half do not. The ones that do have a square root have two square roots. For example, look at the integers mod 7. The numbers 1, 2, and 4 have square roots but the numbers 3, 5, and 6 do not. Here’s a little bit of Python to verify this.

>>> [x**2 % 7 for x in range(1,7)]
[1, 4, 2, 2, 4, 1]

Now if we wanted to count the number of solutions to

y² = x (mod p)

we’d find it equals p. That’s because (0, 0) is a point on the curve, and for (p − 1)/2 values of x there are two values of y, and for the remaining values of x there are no values of y.

To put it another way, if we pick a value of x at random, it will have two square roots or no square roots with roughly equal probability.

Now suppose we pick an x and random and compute x³ + ax + b. We might guess that this process scatters x values in an unbiased way, being no more likely to produce a square than a non-square. Our guess turns out to be approximately correct, though there may be a small bias in favor of squares or non-squares.

Precise bound

Hasse’s theorem gives a precise bound on the number of points on an elliptic curve E. It says

#E = p + 1 − t

where

|t| ≤ 2√p.

So if p is on the order of 2n, the number of points on E is also on the order of 2n because √p is far smaller than p. That’s all you need to know at a high level. The rest of this post will look at a particular example in full detail.

Detailed example

The NIST curve P-384 has been in the news lately. For this curve, the based field is integers modulo a 384-bit prime p and so the number of points on the curve is also 384-bit number. Specifically

p = 2384 − 2128 − 296 + 232 − 1,

and so we know a priori that the number of points on the curve equals 2384 plus or minus a factor less than 2193.

The number of points on the elliptic curve P-384 is known to be

 n = 39402006196394479212279040100143613805079739270465446667946905279627659399113263569398956308152294913554433653942643

See NIST SP 800-186, §3.2.1.4⁠.

We can calculate t at a Python REPL.

>>> n = 39402006196394479212279040100143613805079739270465446667946905279627659399113263569398956308152294913554433653942643
>>>  p = 2**384 - 2**128 - 2**96 + 2**32 - 1
>>> t = p + 1 - n
>>> t
1388124618062372383606759648309780106643088307173319169677

The fact that t is positive means the number of points n on our curve is slightly less than p. Here “slightly” is relative. The value of p − n is a huge number, but it’s small relative to p.

We can verify that |t| is less than our a priori bound of 2√p.

>>> from math import sqrt
>>> t/(2*sqrt(p))
0.11057050503394951

So the actual value of t is about nine times smaller than its bound.

The post The number of points on an elliptic curve first appeared on John D. Cook.

Backdoors, Incompetence, and Elliptic Curves

2026-10-21 23:16:51

The NIST elliptic curves have been in the news twice lately. A few days ago Anthropic posted a paper on the NIST curve seeds that I wrote about here. Now the curves are in the news again, this time over broken certificates.

My post a few days ago was about the possibility that the NSA had added a backdoor to the NIST curves. Anthropic researchers showed that this wasn’t the case. The NSA had not been able to defend how the curve seeds were chosen because they forgot, but researchers from Anthropic were able to reverse engineer the source of the seeds.

So apparently the NIST curves do not have a backdoor, but as with all encryption, the Achilles heel is human error. eMudrha was using the NIST curve P-384  with extremely small keys: 4, 56, 115, 155, 156, 230, and 237. To understand why this is a problem, I’ll briefly describe how elliptic curve cryptography works.

P-384 is an elliptic curve designed by NIST, with input from the NSA, for use in cryptography. The curve has on the order of 2384 points because the base field is the integers mod

p = 2384 − 2128 − 296 + 232 − 1,

a 384-bit prime number. Given a basepoint b on the curve and an exponent e, it is practically impossible to recover e from be. That’s the discrete logarithm problem [1], the basis for the security of elliptic curve cryptography. If and when cryptographically relevant quantum computers arrive, the discrete logarithm problem will be solvable and elliptic curve cryptography will be broken.

Until the advent of quantum computing, elliptic curve cryptography is presumably safe, unless someone does something negligent, such as using a tiny value of e like 4.

If e is chosen at random from a group with on the order of 2384 elements, you cannot recover e from be without brute force. But apparently someone said

Let’s see. Is e = 1? Is it 2? Is it 3? Is it 4? Why yes: e = 4.

Ordinarily this would be hopeless. You can’t just try 2384 exponents. But someone thought to test whether the keys were horrifically bad, and they were right.

 

[1] The exponentiation here is not over the integers, or over the integers mod p. It’s over the group structure of the elliptic curve, and that is what makes it hard to compute the discrete logarithm.

The post Backdoors, Incompetence, and Elliptic Curves first appeared on John D. Cook.