The Goldilocks Prime

In Ham15 Mike Hamburg introduced a class of primes of the form $p={\phi }^2-\phi -1$ and named them Goldilocks primes. In those prime fields $\phi$ satisfies the Golden ratio relation ${\phi }^2=\phi +1$. When $\phi$ is a power of two this allows for very efficient implementation. Unfortunately, the $-1$ at the end means it will have very few power-of-two roots of unity. We'd like to have many so we can do big number theoretic transforms.

We will instead consider primes of the form $p={\phi }^2-\phi +1$. Following the plonky2 project we keep the name Goldilocks prime. Now $\phi$ satisfies a modified Golden ration relation ${\phi }^2=\phi -1$. More importantly, the $+1$ at the end guarantees we will have a $\phi$-th root of unity.

The values of $\phi$ and $p$ where $p$ is prime are the sequences A055494 and A002383 respectively:

$$\begin{array}{rl}\phi & =2,3,4,6,7,9,13,15,16,18,21,22,25,28,\dots ,\\ p& =3,7,13,31,43,73,157,211,241,307,421,\dots \end{array}$$

As mentioned, in the prime field ${\mathbb{F}}_p$ we have ${\phi }^2=\phi -1$. From this further useful relations on the powers of $\phi$ follow:

$$\begin{array}{rlrlrlrl}{\phi }^3& =-1& {\phi }^4& =-\phi & {\phi }^5& =1-\phi & {\phi }^6& =1\end{array}$$

And the pattern repeats. This shows hat $\phi$ is a primitive 6th root of unity. These relations leads to a particularly efficient modular reduction algorithm. First we write the number $n$ in base $\phi$:

$$n=n_0+n_1\cdot \phi +n_2\cdot {\phi }^2+n_3\cdot {\phi }^3+n_4\cdot {\phi }^4+\cdots$$

We get rid of all the terms $n_{>3}$ by folding them into $n_0$, $n_1$, $n_2$ using the relations on ${\phi }^3$ and ${\phi }^6$.

$$\begin{array}{rl}{n}_0^{\prime }& =n_0-n_3+n_6-n_9+\cdots \\ n_1^{\prime }& =n_1-n_4+n_7-n_{10}+\cdots \\ n_2^{\prime }& =n_2-n_5+n_8-n_{11}+\cdots \end{array}$$

To reduce this from three to two terms we use the relation on ${\phi }^2$.

$$\begin{array}{r}{n}_0^{\prime }=n_0-n_2\\ n_1^{\prime }=n_1+n_2\end{array}$$

The reduced number is $n_0+n_1\cdot \phi$. All that remains is carry propagation and some conditional subtractions of the modulus. Note that the reduction requires no multiplication operations.

To multiply two numbers $a,b\in {\mathbb{F}}_p$ we can write them in base $\phi$ and make use of Karatsuba multiplication:

$$\begin{array}{rl}a\cdot b& =(a_0+a_1\cdot \phi )\cdot (b_0+b_1\cdot \phi )\\ & =a_0\cdot b_0+(a_0\cdot b_1+a_1\cdot b_0)\cdot \phi +a_1\cdot b_1\cdot {\phi }^2\\ & =a_0\cdot b_0-a_1\cdot b_1+(a_0\cdot b_1+a_1\cdot b_0+a_1\cdot b_1)\cdot \phi \\ & =a_0\cdot b_0-a_1\cdot b_1+((a_0+a_1)\cdot (b_0+b_1)-a_0\cdot b_0)\cdot \phi \end{array}$$

This requires three full multiplications in base $\phi$ and results in three terms. We reduce this to two terms using the final step of the modular reduction procedure.

Question. What do Barrett and Montgomery reduction look like?

A lucky binary base

To implement base $\phi$ operations efficiently we'd like it to be a power of two, i.e. $\phi =2^k$. There are not many such primes, the only solutions with $k<4000$ are $k=1$, $4$, and $32$. The first to are too small to be practically useful which leaves us with $32$. $\phi =2^{32}$ makes an especially useful prime $p=2^{64}-2^{32}+1$ or written in full

$$p=18446744069414584321=0x\mathtt{ffff}\mathtt{ffff}00000001$$

This is a great size. It is small enough to fit in a 64-bit number, but large enough to do a lot of useful math in. For example, if we have three 32-bit numbers $a,b,c\in [0,2^{32}-1]$ we can compute $a\cdot b+c$ without overflowing:

$$a\cdot b+c\le (2^{32}-1{)}^2+(2^{32}-1)=2^{64}-2^{32}<p$$

When computing on 32/64-bit computers we have some useful interactions with the binary rings:

$$\begin{array}{rlrlrlrlr}{\left[2^{n\cdot 32}\right]}_p& ={\left[{\phi }^n\right]}_p& {\left[2^{64}\right]}_p& ={\left[-1\right]}_{2^{32}}& {\left[-p\right]}_{2^{64}}& ={\left[-1\right]}_{2^{32}}& {\left[\phi \right]}_{2^{32}}& =0& \end{array}$$

The first identity gives us a very efficient way to implement bit-shifts. This in turn leads to an efficient inversion algorithm from the half-extended binary GCD:

def inverse(a):
    (u, a) = (MODULUS, a)
    (t0, t1) = (0, 1)
    twos = trailing_zeros(v)
    v >>= twos
    twos += 96 # 2^96 = -1
    while u != v:
        u -= v
        t0 += t1
        count = trailing_zeros(u)
        u >>= count
        t1 <<= count
        twos += count
        if u < v:
            (u, v) = (v, u)
            (t0, t1) = (t1, t0)
            twos += 96 # 2^96 = -1
    return t0 << ((191 * twos) % 192) # 2^-1 = 2^191, 2^192 = 1

Multiplicative group

The multiplicative group has order $p-1$ and factors into a power of two and, interestingly, all five known Fermat primes.

$$p-1=2^{32}\cdot 3\cdot 5\cdot 17\cdot 257\cdot 65537$$

The factor $2^{32}$ means we can do efficient NTTs of power-of-two sizes using the Cooley-Tukey algorithm. The Fermat primes are not just curious, they allow us to use Rader's FFT algorithm to turn the Fermat prime sized NTTs into power-of-two NTTs. With efficient NTTs for all prime factors we can use the Good-Thomans algorithm to efficiently combine these without requirng twiddle factors. Overall this gives decent NTT performance for any divisor of $p-1$.

But what if we want a different size? The usual method is to pad the data until it is a power of two size. This is wasteful as can mean almost doubling the size. Instead we could go beyond powers of two and include the other small primes $3$, $5$ things improve quite a bit

For example if we do a $2^{32}\cdot 3$ mixed radix, then in half of the cases it is 25% smaller than if we only use $2^{32}$. If we also add $5$ then in a quarter of the cases it is an extra 16% smaller. These two seem worthwhile. Adding $17$ still gives a visible improvement, but is probably not worth the implementation complexity. And $257$ and $65537$ do not contribute anything noticeable.

NTTs require lot's of multiplications by roots of unity. Goldilocks again helps us make this efficient. We already found that $\phi =2^{32}$ is a $6$th root of unity, and multiplying by it is simply a left shift of 32 bits. But if $2^{32}$ is a $6$th root, then $2$ itself must be a $6\cdot 32=192$th root. Great! This means multiplying by $192$th roots is just cheap bit shift. And since $192=2^6\cdot 3$ has many divisors we can now create a whole set of $n$-th roots ${\omega }_n$ that can be implemented using bit shifts:

$$\begin{array}{rlrlrlrlrl}{\omega }_2& =2^{96}& {\omega }_6& =2^{32}& {\omega }_{16}& =2^{12}& {\omega }_{48}& =2^4& {\omega }_{192}& =2\\ {\omega }_3& =2^{64}& {\omega }_8& =2^{24}& {\omega }_{24}& =2^8& {\omega }_{64}& =2^3& & \\ {\omega }_4& =2^{48}& {\omega }_{12}& =2^{16}& {\omega }_{32}& =2^6& {\omega }_{96}& =2^2& & \end{array}$$

We can go even further and construct a $384$th root using two bit shifts, but for that we need a trick. Schönhage observed that given an $8$th root ${\omega }_8$, the value ${\omega }_8+{\omega }_8^7$ is a square root of $2$. Let's quickly check this

$$({\omega }_8+{\omega }_8^7{)}^2={\omega }_8^2+{\omega }_8^{7\cdot 2}+2\cdot {\omega }_8^{1+7}={\omega }_8^2-{\omega }_8^2+2=2$$

In Goldilocks $2^{24}$ is an $8$th root and $2^{96}=-1$ so $n=2^{24}-2^{72}$ is a square root of two. Furthermore since $n^2=2$ and $2$ is a $192$th root, it follows that $n$ is a $384$th root ${\omega }_{384}$. If we don't want that factor of $3$ we can compute

$${\omega }_{128}={\omega }_{384}^3={\omega }_{384}^2\cdot {\omega }_{384}=2\cdot (2^{24}-2^{72})=2^{25}-2^{73}$$

This also shows that all even powers of ${\omega }_{384}$ have one bit set and all odd powers two.

$${\omega }_{384}^i=\{\begin{array}{ll}{2}^k& i=2\cdot k\\ 2^{k+24}-2^{k+73}& i=2\cdot k+1\end{array}$$

This $\sqrt{2}$ trick is useful, let's understand it better and see if we can take it further.

Schönhage's trick

To intuitively understand why this works, forget for the moment that we are working in a prime field and imagine the 8th roots of unity in the complex plane.

We can construct the sum ${\omega }_8+{\omega }_8^7$ by adding vectors. Pythagoras theorem shows it is $\sqrt{2}$. This procedure generalizes this to $n$-roots with ${\omega }_n+{\omega }_n^{n-1}$, where the angles between roots are $\frac{2\cdot \pi }{n}$. A bit of trigonometry shows the solutions

$${\omega }_n+{\omega }_n^{n-1}=2\cdot \cos \frac{2\cdot \pi }{n}$$

For certain values of $n$ there are simple exact values for $\cos$, see a table of exact trigonometric values. We indeed get $\sqrt{2}$ and a few more. Most look quite uninteresting, but these two I like

$$\begin{array}{rlrl}{\omega }_{12}+{\omega }_{12}^{11}& =\sqrt{3}& {\omega }_{10}+{\omega }_{10}^9& =\frac{1+\sqrt{5}}{2}\end{array}$$

The last one is also known as the Golden ratio. It turns out these work in the Goldilocks field as well. Since ${\omega }_{12}=2^{16}$ we have $\sqrt{3}=2^{16}-2^{80}$. For ${\omega }_{10}$ I don't know a multiplication free construction.

Question. It is possible to geometrically construct $\sqrt[4]{2}$ using the geometric mean theorem. Is there a construction that leads to an efficient (multiplication free) implementation?

Adam P. Goucher explains that it is not possible to construct $\sqrt[4]{2}$ as a rational function of roots of unity.

Generators

So far we've cleverly constructed roots up to ${\omega }_{384}$, but for our NTTs we need all the roots up to the $2^{32}$th. The standard way to do this is to find a generator of the multiplicative group. A generator is a field element $g$ such that all non-zero numbers can be written as powers $g^i$ for some $i$. There are many such numbers and the smallest generator in Goldilocks in $7$. By Fermat's little theorem we have $g^{p-1}=1$ and thus $g$ is a primitive $(p-1)$-th root of unity. Since this is the highest possible root of unity, we can derive all others from it using

$${\omega }_n=g^{\frac{p-1}{n}}$$

It is convention to use the smallest generator, but if we compute ${\omega }_{192}$ with $g=7$ we get $35184372080640$ instead of $2$. This happens because $n$-th primitive roots are not unique and there are many choices that work. While this is still a power of two, ${\omega }_{192}=2^{77}$, it makes things more complicated than they need to be and we get a more complicated formula for ${\omega }_{384}^i$.

A priory there is no reason to pick one choice of roots over the other. But we do need our choices to be consistent, meaning we want

$${\omega }_n={\omega }_{k\cdot n}^k$$

to hold where applicable. Using a generator full fills that requirement. We'd also like a choice where we have an easy formula for multiplication by ${\omega }_{384}^i$ in terms of bit shifts. The smallest generator that satisfies this is $554$. But this one has ${\omega }_{384}=-({\omega }_8+{\omega }_8^7)$, it takes the negative square root of two. The smallest generator that picks the positive root is $2717$.

Question. Are there larger values of $k$ for which $2^{2\cdot k}-2^k+1$ is prime or is $32$ the largest?

Further references

• Gou21. Adam P. Goucher (2021). "An efficient prime for number-theoretic transforms".
• Po22. Thomas Pornin (2022).

https://github.com/mir-protocol/plonky2/issues/1#issuecomment-809323437