Math

Quality metrics

Posted on by John

Sometimes you’ll hear that a process has so many nines of reliability or that the error rate is so many sigmas. A few years ago I wrote a post on converting between nines and sigmas. See that post for details, approximations, etc.

Here I’d like to post a new graph that I believe is an improvement on the graph in the original post.

The grid lines make it easier to see how to convert between nines and sigmas. For example, 5 nines corresponds to a little more than 4 sigmas.

It’s curious that 3 nines is approximately 3 sigmas, and 9 nines is very nearly 6 sigmas. There are no more points on the red curve with near integer coordinates until you go out to 72 nines, which is nearly 18 sigmas.

Turning trig identities into Fibonacci identities

Posted on by John

In 2013, John Conway and Alex Ryba published a brief note [1] on how to convert identities involving sine and cosine into identities involving Fibonacci and Lucas numbers.

Fibonacci and Lucas

The Fibonacci numbers Fn are defined by F0 = 0, F1 = 1, and

Fn+2 = Fn + Fn+1

for n > 1. Similarly, Lucas numbers Ln are defined by the same recurrence relation

Ln+2 = Ln + Ln+1

but starting with L0 = 2 and L1 = 1.

The recipe

The recipe for turning trigonometric identities into Fibonacci and Lucas number identities is as follws.

  1. Arguments become subscripts.
  2. Replace sines by ½ in Fn.
  3. Replace cosines by ½ in Ln.
  4. Insert a factor of (−5)k in front of every term that contains 2k or 2k + 1 sine terms.

The first step is admittedly unclear. It’s unclear in [1] as well. If sine or cosine takes an argument of

pα + qβ

then this becomes the subscript

pa + qb

where the Latin letters are integers and the Greek letters are angles.

There’s no proof in [1] why the recipe should work, but an earlier paper [2] gives more details.

Examples

Double angle

For a simple example, let’s start with the double angle identity

sin 2θ = 2 sin θ cos θ.

This becomes

½ i2n F2n = 2 (½ in Fn)(½ in Ln)

which simplifies to

F2n = Fn Ln.

Sum angle

For a little more involved example, let’s look at the sum identity

cos(θ + φ) = cos θ cos φ − sin θ sin φ.

This becomes

½ ia + b La + b = ½ ia La ½ ib Lb −(−5) ½ ia Fa ½ ib Fb

which simplifies to

2La + b = La Lb  + 5 Fa Fb.

Triple angle

In the previous examples the powers of i in the recipe above cancelled out. This example will show that they are necessary.

We start with the triple angle identity

sin 3θ = 3 sin θ − 4 sin³ θ.

This becomes

½ i3n F3n = 3 (½) in Fn − 4(−5)(½ in Fn)³

which simplifies to

F3n = 3 (−1)n Fn + 5 (Fn)³.

Demonstration

The following Python code demonstrates that the Fibonacci / Lucas identities above are correct.

def F(n):
    if n == 0: return 0
    if n == 1: return 1
    return F(n-1) + F(n-2)

def L(n):
    if n == 0: return 2
    if n == 1: return 1
    return L(n-1) + L(n-2)

for n in range(10):
    assert(F(2*n) == F(n) * L(n))

for n in range(10):
    for m in range(10):
        assert(2*L(m + n) == L(m)*L(n) + 5*F(m)*F(n))

for n in range(10):
    assert(F(3*n) == 3*(-1)**n * F(n) + 5*F(n)**3)

Related posts

[1] John Conway and Alex Ryba. Fibonometry. The Mathematical Gazette, November 2013, pp. 494–495

[2] Barry Lewis, Trigonometry and Fibonacci Numbers, The Mathematical Gazette, July 2007, pp. 216–226

Ethereum’s consensus layer elliptic curve

Posted on by John

I’ve written before about Bitcoin’s elliptic curve and Monero’s elliptic curve. In the Monero post I wrote “Bitcoin and Ethereum use the elliptic curve secp256k1.” That’s true, but it’s incomplete. Ethereum does use the elliptic curve secp256k1 for digital signatures, as does Bitcoin, but Ethereum also uses a different elliptic curve for its consensus layer.

Ethereum’s consensus layer uses the elliptic curve BLS12-381. This post will say a little bit about this curve, starting with unpacking the cryptic name. I won’t go into the advanced math behind the curve’s design but instead will focus on concrete calculations in Python to try to make the curve more tangible. The same curve is used in other cryptocurrencies, including Zcash.

First of all, BLS stands for Paulo Barreto, Ben Lynn, and Michael Scott, the developers of a family of elliptic curves known as the BLS curves, including BLS12-381. Incidentally, in the context of cryptography, BLS can refer to the BLS curves or to BLS signatures. The “L” is the same person in both instances, but the “B” and “S” in BLS signatures refer to Dan Boneh and Hovav Shacham.

The 381 in BLS12-381 refers to the fact that it is defined over a finite field whose order is a 381-bit number, namely

p = 0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab

We can verify that p is prime and that it has 381 bits.

>>> from sympy import isprime
>>> isprime(p)
True
>>> 2**380 < p < 2**381
True

The value of p comes from applying a certain polynomial to a parameter z with low Hamming weight, i.e. with lots of zeros in its binary representation. Why this matters is beyond the scope of this post, but we can show that it’s true.

>>> z = -0xd201000000010000
>>> p == (z-1)**2*(z**4 - z**2 + 1)//3 + z
True

The elliptic curve BLS12-381 is the set of points satisfying

y² = x³ + 4 mod p.

The 12 in BLS12-381 refers to an embedding degree k that we’ll get to shortly.

The elliptic curve BLS12-381 is pairing friendly, which is the reason for its use in the Ethereum consensus layer. This layer uses pairing-based cryptography to aggregate signatures. I may write more about that someday, but not today.

As I wrote a couple months ago,

An elliptic curve E/Fq is said to be pairing-friendly if r divides qk − 1 for some small k. Here r is the size of the largest prime-order subgroup of the curve.

In the case of BLS12-381, r = z4 − z2 + 1 and r is a 255-bit prime number.

>>> r = z**4 - z**2 + 1
>>> isprime(r)
True
>>> 2**254 < r < 2**255
True

And now we can verify that that the embedding degree is 12, showing the BLS12-381 is a pairing-friendly curve.

>>> (p**12 - 1) % r
0

So what is being paired with what? And what is being embedded into what? The group G1 of order r, is a subgroup of BLS12-381. It is paired with another group G2 also of order r, and there is a bilinear mapping from (G1, G2) to the multiplicative group of the finite field with p12 elements. For more details, see the section on BLS12-381 in Ben Edginton’s book Upgrading Ethereum: A technical handbook on Ethereum’s move to proof of stake and beyond.

Related posts

Inverting matrices and bilinear functions

Posted on by John

The inverse of the matrix

M = \begin{bmatrix} a & b \\ c & d \end{bmatrix}

is the matrix

M^{-1} = \frac{1}{|M|} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix} = \frac{1}{ad - bc} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix}

assuming ad − bc ≠ 0.

Also, the inverse of the bilinear function (a.k.a. Möbius transformation)

f(z) = \frac{az + b}{cz + d}

is the function

f^{-1}(z) = \frac{dz - b}{-cz + a}

again assuming ad − bc ≠ 0.

The elementary takeaway is that here are two useful equations that are similar in appearance, so memorizing one makes it easy to memorize the other. We could stop there, but let’s dig a little deeper.

There is apparently an association between 2 × 2 matrices and Möbius transformations

\frac{az + b}{cz + d} \leftrightarrow \begin{bmatrix} a & b \\ c & d \end{bmatrix}

This association is so strong that we can use it to compute the inverse of a Möbius transformation by going to the associated matrix, inverting it, and going back to a Möbius transformation. In diagram form, we have the following

Now there are a few loose ends. First of all, we don’t really have a map between Möbius transformations and matrices per se; we have a map between a particular representation of a Möbius transformation and a 2 × 2 matrix. If we multiplied a, b, c, and d in a Möbius transformation by 10, for example, we’d still have the same transformation, just a different representation, but it would go to a different matrix.

What we really have is a map between Möbius transformations and equivalence classes of invertible matrices, where two matrices are equivalent if one is a non-zero multiple of the other. If we wanted to make the diagram above more rigorous, we’d replace ℂ2×2 with PL(2, ℂ), linear transformations on the complex projective plane. In sophisticated terms, our map between Möbius transformations and matrices is an isomorphism between automorphisms of the Riemann sphere and PL(2, ℂ).

Möbius transformations act a lot like linear transformations because they are linear transformations, but on the complex projective plane, not on the complex numbers. More on that here.

Generate random points inside a sphere

Posted on by John

For reasons I don’t understand, it’s more common, at least in my experience, to see references on how to generate points on a sphere than in a sphere. But the latter is easily obtained from the former.

The customary way to generate points on a sphere in n dimensions is to generate points from an n-dimensional normal (Gaussian) distribution and normalize by dividing by the length. And the way to generate samples from an n-dimensional normal is simply to generate a list of n one-dimensional normals [1].

The way to generate points in a sphere is to generate a point S on the sphere then randomly adjust its radius r. If you generate a uniform random point in [0, 1] and set r = u then you’ll over-sample points near the origin. Instead, the thing to do is to set

r = u1/n

That’s because volume is proportional to rn. Taking the nth root of u mades the volume inside a sphere of radius r uniformly distributed.

Demonstration

To illustrate this, I’ll use the same technique as the post on sampling from a tetrahedron: generate random points in a sphere and see whether the right proportion fall inside a (hyper)cube inside the sphere.

We’ll work in n = 4 dimensions and use a sphere of radius R = 5. Our test cube will be the cube with all vertex coordinates equal to 1 or 2. We’ll have to use a lot of sample points. (This is why naive Monte Carlo breaks down in high dimensions: it may be that few if any points land in the region of interest.)

Here’s the code to do the sampling.

import numpy as np
from scipy.stats import norm

def incube(v):
    retval = True
    for x in v:
        retval = retval and 1 <= x <= 2
    return retval

def sample(R, n):
    x = norm.rvs(size=n)
    x /= np.linalg.norm(x)
    u = np.random.random()
    x *= R*u**(1/n)
    return x

R = 5
n = 4
N = 100_000_000
V = R**n * np.pi**(n/2) / gamma(n/2 + 1)
print(V)

count = 0

for _ in range(N):
    v = sample(R, n)
    if incube(v):
        count += 1        

print(count/N, 1/V)

This prints

0.0003252 0.000324

The two results agree to four decimal places, which is about what we’d expect given 1/√N = 0.0001.

[1] If you’re porting the sampler above to an environment where you don’t have a function like norm.rvs for generating normal random samples, you can roll your own using the Box-Muller transformation. For example, here is a C++ implementation from scratch.

More on Carmichael

Posted on by John

This morning I took an old blog post and turned it into an X thread. I think the thread is easier to read. More expository and less rigorous.

The post and thread look at generalizations of the fact that every integer and its fifth power end in the same digit. The conclusion is that n and nk end in the same digit base b if b is square-free and k = φ(b) + 1. So, for example, 10 is square-free (i.e. not divisible by a non-trivial square) and φ(10) + 1 = 5.

Benjamin Clark replied suggesting looking at λ(b) + 1 in addition to φ(n) + 1.

Euler and Carmichael

To back up a bit, φ(n) is Euler’s totient function, the number of positive integers less than and relatively prime to n. Robert Carmichael’s totient function λ(n) is a closely related function, one that replaced Euler’s function in implementations of RSA encryption—more specifically in RSA decryption—something I wrote about earlier this week.

Euler’s generalization of Fermat’s little theorem says if a is relatively prime to m, then

aφ(n) = 1 (mod n).

Carmichael’s totient λ(n) is the smallest number such that

aλ(n) = 1 (mod n)

when a is relatively prime to n.

Sometimes φ(n) = λ(n), namely for n = 1, 2, 4, and every odd prime power. Otherwise λ(n) is a proper divisor of φ(n).

Carmichael’s conjecture

Carmichael conjectured that for every n there is at least one other integer m ≠ n such that φ(m) = φ(n). He proved that this is true at least for n less than 1037. The conjecture is now known to be true for n less than 101010.

Differential equation on a doughnut

Posted on by John

Here’s a differential equation from [1] that’s interesting to play with.

\begin{align*} x^\prime &= -my + nxz \\ y^\prime &= \phantom{-} mx + nyz \\ z^\prime &= (n/2)(1 + z^2 - x^2 - y^2) \\ \end{align*}

Even though it’s a nonlinear system, it has a closed-form solution, namely

\begin{align*} x(t) &= \frac{2a\cos(mt) - 2b \sin(mt)}{\Delta - 2c \sin(n t) + (2 - \Delta) \cos(n t)} \\ y(t) &= \frac{2a \sin(mt) + 2b \cos(mt)} {\Delta - 2 c \sin(n t) + (2 - \Delta) \cos(n t)} \\ z(t) &= \frac{2c \cos(nt) + (2 - \Delta) \sin(nt)}{\Delta - 2c \sin(nt) + (2 - \Delta) \cos(n t)} \\ \end{align*}

where (a, b, c) is the solution at t = 0 and Δ = 1 + a² + b² + c².

The solutions lie on the torus (doughnut). If m and n are coprime integers then the solutions form a closed loop. If the ratio m/n is not rational then the solutions are dense on the torus.

Here’s an example with parameters a = 1, b = 1, c = 3, m = 3, and n = 5.

And now with parameters a = 1, b = 1, c = 0.3, m = 4, and n = 5.

And finally with parameters a = 1, b = 1, c = 0.3, m = π, and n = 5.

Note that when m = 2 and n = 3 the trajectory traces out a trefoil knot.

Related posts

[1] Richard Parris. A Three-Dimensional System with Knotted Trajectories. The American Mathematical Monthly, Vol. 84, No. 6 (Jun. – Jul., 1977), pp. 468–469

RSA with multiple primes

Posted on by John

Typically RSA public keys are the product of two large primes, n = pq. But there’s no reason they couldn’t be the product of say three primes, n = pqr, or more primes, as long as φ(n), or λ(n), is calculated correctly.

Encryption is done the same way. Decryption could be done the same way, except there is the opportunity for it to be more efficient. The trick is to use the CRT (Chinese Remainder Theorem) in a way similar to Garner’s algorithm. This is why RSA with multiple primes is sometimes used for digital signatures.

The difficulty of factoring n using the GNFS algorithm doesn’t change depending on the number of factors n has, as long as all the factors are sufficiently large, far too large to find using trial division.

Daniel Bernstein’s post-quantum RSA paper was based on keys that are the product of a large number of 4096-bit primes. This way all the arithmetic is carried out modulo 4096-bit primes, not modulo terabyte primes.

A quiet change to RSA

Posted on by John

An RSA public key is a pair of numbers (e, n) where e is an exponent and n = pq where p and q are large prime numbers. The original RSA paper said choose a private key d and compute e. In practice now we choose e and compute d. Furthermore, e is now almost always 65537 for reasons given here. So the public key is essentially just n.

Euler’s totient function

The relationship between the exponent and the private decryption key in the original RSA paper was

ed = 1 mod φ(n).

It is easy to compute e given d, or d given e, when you know Euler’s totient function of n,

φ(n) = (p − 1)(q − 1).

The security of RSA encryption rests on the assumption that it is impractical to compute φ(n) unless you know p and q.

Carmichael’s totient function

Gradually over the course of several years, the private key d changed from being the solution to

ed = 1 mod φ(n)

to being the solution to

ed = 1 mod λ(n)

where Euler‘s totient function φ(n) was replaced with Carmichael‘s totient function λ(n).

The heart of the original RSA paper was Euler’s generalization of Fermat’s little theorem which says if a is relatively prime to m, then

aφ(n) = 1 (mod n)

Carmichael’s λ(n) is defined to be the smallest number that can replace φ(n) in the equation above. It follows that λ(n) divides φ(n).

Why the change?

Using Carmichael’s totient rather than Euler’s totient results in smaller private keys and thus faster decryption. When n = pq for odd primes p and q,

λ(n) = lcm( (p − 1), (q − 1) ) = (p − 1)(q − 1) / gcd( (p − 1), (q − 1) )

so λ(n) is smaller than φ(n) by a factor of gcd( (p − 1), (q − 1) ). At a minimum, this factor is at least 2 since p − 1 and q − 1 are even numbers.

However, an experiment suggests this was a trivial savings. When I generated ten RSA public keys the gcd was never more than 8.

from sympy import randprime, gcd

for _ in range(10):
    p = randprime(2**1023, 2**1024)
    q = randprime(2**1023, 2**1024)
    print(gcd(p-1, q-1))

I repeated the experiment with 100 samples. The median of the gcd’s was 2, the mean was 35.44, and the maximum was 2370. So the while gcd might be moderately large, but it is usually just 2 or 4.

Better efficiency

The efficiency gained from using Carmichael’s totient is minimal. More efficiency can be gained by using Garner’s algorithm.

Fermat primes and tangent numbers

Posted on by John

Fermat numbers

The nth Fermat number is defined by

F(n) = 2^{2^n} + 1

Pierre Fermat conjectured that the F(n) were prime for all n, and they are for n = 0, 1, 2, 3, and 4, but Leonard Euler factored F(5), showing that it is not prime.

Tangent numbers

The nth tangent number is defined by the Taylor series for tangent:

\tan(z) = \sum_{n=0}^\infty T(n) \frac{z^n}{n!}

Another way to put it is that the exponential generating function for T(n) is tan(z).

Fermat primes and tangent numbers

Here’s a remarkable connection between Fermat numbers and tangent numbers, discovered by Richard McIntosh as an undergraduate [1]:

F(n) is prime if and only if F(n) does not divide T(F(n) − 2).

That is, the nth Fermat number is prime if and only if it does not divide the (F(n) − 2)th tangent number.

We could duplicate Euler’s assessment that F(5) is not prime by showing that 4294967297 does not divide the 4294967295th tangent number. That doesn’t sound very practical, but it’s interesting.

Update: To see just how impractical the result in this post would be for testing whether a Fermat number is prime, I found an asymptotic estimate of tangent numbers on OEIS,  and estimated that the 4294967295th tangent number has about 80 billion digits.

[1] Richard McIntosh. A Necessary and Sufficient Condition for the Primality of Fermat Numbers. The American Mathematical Monthly, Vol. 90, No. 2 (Feb., 1983), pp. 98–99