Introduction

This is a (very) shortened version of my Bachelor thesis in Computer Science.

Miller-Rabin primality test

The Miller-Rabin primality test checks for a number's primality by ensuring that there an non-trivial square roots of $11$ modulo this number. For a positive odd integer $nn$ and number of rounds $rr$ the test works like this:

• Factor $n-1n\; -\; 1$ as $2^{k}d2^kd$.
• For each round repeat this steps:
  1. Select a random value $aa$ in the interval $\{2,\dots ,n-2\}\backslash \{2,\; \backslash ldots,\; n\; -\; 2\backslash \}$ called base.
  2. Calculate $a^d(\mathrm{m}\mathrm{o}\mathrm{d}n)a^d\; \backslash pmod\; n$, if this is equal to $11$ modulo $nn$ go back to step 1, otherwise proceed to step 3.
  3. If $a^{2^{i}d}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}n)a^\{2^id\}\; \backslash equiv\; -1\; \backslash pmod\; n$ for any go back to step 1 (if it is not the last round), otherwise exit and declare $nn$ composite.
• If no round deemed $nn$ composite the test deems it prime.

Arnault's method

Arnault's method constructs an odd integer $nn$ product of an odd number of primes $p_1,p_2,\dots ,p_{h}p\_1,p\_2,\backslash ldots,p\_h$ such that $nn$ is a strong pseudoprimes to a given set of prime bases $\mathcal{A}=\{a_1,a_2,\dots ,a_t\}\backslash mathcal\{A\}\; =\; \backslash \{a\_1,a\_2,\backslash ldots,a\_t\backslash \}$. The method construct $p_{1}p\_1$ via a system of equations and then each $p_{i}p\_i$ as $k_i(p_1-1)+1k\_i(p\_1\; -\; 1)\; +\; 1$ for $k_{i}k\_i$ randomly chosen to be coprime to each base (except for $k_1=1k\_1\; =\; 1$) for $1\le i\le h1\; \backslash le\; i\; \backslash le\; h$. For $nn$ to be a strong pseudoprimes to bases $a_1,a_2,\dots ,a_{t}a\_1,a\_2,\backslash ldots,a\_t$ it is sufficient that each of its prime factors $pp$ satisfies $\left(\frac{a_i}{p}\right)=-1\backslash left(\backslash frac\{a\_i\}\{p\}\backslash right)\; =\; -1$ (so $a_{i}a\_i$ is a quadratic nonresidue modulo $pp$) for $1\le i\le t1\; \backslash le\; i\; \backslash le\; t$, where $\left(\frac{x}{y}\right)\backslash left(\backslash frac\{x\}\{y\}\backslash right)$ is the Legendre symbol of $xx$ modulo $yy$.

Arnault's method starts by generating sets $S_{a}S\_a$ for each base $aa$ such that:

Since $p_{1}p\_1$ has to satisfy the conditions of each set and so do $p_2,\dots ,p_{h}p\_2,\backslash ldots,p\_h$ stricter conditions can be derived:

From each of these sets random residues $r_1,r_2,\dots ,r_{t}r\_1,r\_2,\backslash ldots,r\_t$ are chosen to solve a system of equations:

The CRT tells us a solution is not guaranteed since moduli are not coprime so we may have to try several times (giving up after a few iterations to not fall into the Sunk cost fallacy).

To restrict the field of search a few more conditions can be added, resulting in polynomials for which $p_{1}p\_1$ has to be a root:

After getting a solution $(r,m)(r,\; m)$ we look for values of $qq$ that make $p_1=r+qmp\_1\; =\; r\; +\; qm$ prime (by simply bruteforcing) and calculate $p_2,\dots ,p_{h}p\_2,\backslash ldots,p\_h$ as described above, if all are prime their product $nn$ will be a strong pseudoprimes to the bases $\mathcal{A}=\{a_1,a_2,\dots ,a_t\}\backslash mathcal\{A\}\; =\; \backslash \{a\_1,a\_2,\backslash ldots,a\_t\backslash \}$.

Optimizations

Some sections of Arnault's method slow it down significantly, but with a few clever tricks some optimizations can be found.

Finding $p_{1}p\_1$

To get $p_{1}p\_1$ we must solve the system of equations we saw above, which might take several iterations of the random selection of residues, but a faster strategy is in front of us: instead of picking all residues at once we can select them one at a time, repicking only when there is no solution and continuing like so until the end ^[1].

Constructing $S_{a}S\_a$

While implementing Arnault's method I noticed a pattern I desperately had to confirm: for $a>2a\; >\; 2$ each set $S_{a}S\_a$ contains exactly $a-1a\; -\; 1$ elements. After some testing which seemed to confirm my hypothesis I wanted to prove it, which eventually resulted in a two part proof, but first some lemmas:

1. Dirichlet's theorem an arithmetic progressions states that there are infinitely many primes of the form $r+kmr\; +\; km$ for coprime integers $rr$ and $mm$.
2. The law of quadratic reciprocity states that $\left(\frac{p}{q}\right)=(-1{)}^{\frac{p-1}{2}\frac{q-1}{2}}\left(\frac{q}{p}\right)\backslash left(\backslash frac\{p\}\{q\}\backslash right)\; =\; (-1)^\{\backslash frac\{p\; -\; 1\}\{2\}\backslash frac\{q\; -\; 1\}\{2\}\}\backslash left(\backslash frac\{q\}\{p\}\backslash right)$.

The proof is split in two cases, depending on the value of $aa$ modulo $44$:

1. If $a\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}4)a\; \backslash equiv\; 1\; \backslash pmod\; 4$ we have $\left(\frac{a}{p}\right)=\left(\frac{p}{a}\right)\backslash left(\backslash frac\{a\}\{p\}\backslash right)\; =\; \backslash left(\backslash frac\{p\}\{a\}\backslash right)$ and since $aa$ is fixed for the set $S_{a}S\_a$ it is easier to work with the right-hand side. In $\{1,\dots ,a-1\}\backslash \{1,\; \backslash ldots,\; a\; -\; 1\backslash \}$ there are exactly $\frac{a-1}{2}\backslash frac\{a\; -\; 1\}\{2\}$ quadratic residues and quadratic nonresidues, we define sets $XX$ and $YY$ with $XX$ containing the odd residues and $YY$ the even ones. Since $S_{a}S\_a$ contains primes reduced modulo $4a4a$ we do not cannot allow even values, meaning elements from the sets $\{x+a\mid x\in X\},\{x+3a\mid x\in X\},\{y\mid y\in Y\},\{y+2a\mid y\in Y\}\backslash \{x\; +\; a\; \backslash mid\; x\; \backslash in\; X\backslash \},\; \backslash \{x\; +\; 3a\; \backslash mid\; x\; \backslash in\; X\backslash \},\; \backslash \{y\; \backslash mid\; y\; \backslash in\; Y\backslash \},\; \backslash \{y\; +\; 2a\; \backslash mid\; y\; \backslash in\; Y\backslash \}$, leaving to fill $S_{a}S\_a$: $X,X^{\mathrm{\prime }}=\{x+2a\mid x\in X\},Y^{\mathrm{\prime }}=\{y+a\mid y\in Y\},Y^{\mathrm{\prime }\mathrm{\prime }}=\{y+3a\mid y\in Y\}X,\; X\text{'}\; =\; \backslash \{x\; +\; 2a\; \backslash mid\; x\; \backslash in\; X\backslash \},\; Y\text{'}\; =\; \backslash \{y\; +\; a\; \backslash mid\; y\; \backslash in\; Y\backslash \},\; Y\text{'}\text{'}\; =\; \backslash \{y\; +\; 3a\; \backslash mid\; y\; \backslash in\; Y\backslash \}$, which are trivially disjoint. Therefore $\mid S_a\mid =\mid X\cup X^{\mathrm{\prime }}\cup Y^{\mathrm{\prime }}\cup Y^{\mathrm{\prime }\mathrm{\prime }}\mid =2\mid X\mid +2\mid Y\mid =2(\mid X\mid +\mid Y\mid )=2\frac{a-1}{2}=a-1\backslash lvert\; S\_a\; \backslash rvert\; =\; \backslash lvert\; X\; \backslash cup\; X\text{'}\; \backslash cup\; Y\text{'}\; \backslash cup\; Y\text{'}\text{'}\; \backslash rvert\; =\; 2\backslash lvert\; X\; \backslash rvert\; +\; 2\backslash lvert\; Y\; \backslash rvert\; =\; 2(\backslash lvert\; X\; \backslash rvert\; +\; \backslash lvert\; Y\; \backslash rvert)\; =\; 2\backslash frac\{a\; -\; 1\}\{2\}\; =\; a\; -\; 1$.

2. For $a\equiv 3(\mathrm{m}\mathrm{o}\mathrm{d}4)a\; \backslash equiv\; 3\; \backslash pmod\; 4$ we do something similar by splitting the set $\{1,\dots ,4a-1\}\backslash \{1,\; \backslash ldots,\; 4a\; -\; 1\backslash \}$ into two disjoint sets: , with $\mid X\mid =\mid Y\mid =a\backslash lvert\; X\; \backslash rvert\; =\; \backslash lvert\; Y\; \backslash rvert\; =\; a$ trivially. Notice that both sets contain all residues modulo $aa$ exactly once and: elements of $XX$ are congruent to $11$ modulo $44$, implying $\left(\frac{a}{x}\right)=\left(\frac{x}{a}\right),\mathrm{\forall }x\in X\backslash left(\backslash frac\{a\}\{x\}\backslash right)\; =\; \backslash left(\backslash frac\{x\}\{a\}\backslash right),\; \backslash forall\; x\; \backslash in\; X$, meaning in $XX$ there are $\frac{a-1}{2}\backslash frac\{a\; -\; 1\}\{2\}$ quadratic nonresidues; meanwhile elements of $YY$ are congruent to $33$ modulo $44$, implying $\left(\frac{a}{y}\right)=-\left(\frac{y}{a}\right),\mathrm{\forall }y\in y\backslash left(\backslash frac\{a\}\{y\}\backslash right)\; =\; -\backslash left(\backslash frac\{y\}\{a\}\backslash right),\; \backslash forall\; y\; \backslash in\; y$, meaning we need the quadratic residues modulo $aa$ which are once again $\frac{a-1}{2}\backslash frac\{a\; -\; 1\}\{2\}$. Combining $XX$ and $YY$ to get $S_{a}S\_a$ we finally get $\mid S_a\mid =\mid X\mid +\mid Y\mid =a-1\backslash lvert\; S\_a\; \backslash rvert\; =\; \backslash lvert\; X\; \backslash rvert\; +\; \backslash lvert\; Y\; \backslash rvert\; =\; a\; -\; 1$.

The steps used in the proof can be made into an algorithm which constructs each set in $\mathcal{O}(a)\backslash mathcal\{O\}(a)$ time:

import gmpy2

def find_Sa(a):
    assert gmpy2.is_prime(a)

    if a % 4 == 1:
        X = set(p for p in range(3, a, 2) if gmpy2.legendre(p, a) == -1)
        Y1 = set(p for p in range(2 + a, 2*a, 2) if gmpy2.legendre(p, a) == -1)
        X1 = set(x + 2*a for x in X)
        Y2 = set(y + 2*a for y in Y1)

        Sa = X | Y1 | X1 | Y2
    else:
        remove = 1 + 4*(-1 * pow(4, -1, a) % a)

        X = set((1 + 4*t) for t in range(a)) - set((remove,))

        X_qnr = set(filter(lambda x: gmpy2.legendre(x, a) == -1, X))
        Y_qr = set((x + 2*a) % (4*a) for x in (X - X_qnr))

        Sa = X_qnr | Y_qr

    return Sa

The naive strategy tests the condition $\left(\frac{a}{p}\right)=-1\backslash left(\backslash frac\{a\}\{p\}\backslash right)\; =\; -1$ for primes in increasing fashion until some arbitrary threshold, which may not find all $a-1a\; -\; 1$ of them; stopping instead once $a-1a\; -\; 1$ have been found is what the figure is showcasing, which is $\mathrm{\Omega }(a\log a)\backslash Omega(a\backslash log\; a)$ (empirically).

Conclusions

Fooling a primality test has huge consequences, for example a discrete logarithm can be solved with a much faster method if the modulus used is a composite with known factorization, jeopardizing protocols such as the Diffie-Hellman key exchange which is ubiquitous in public-key cryptography. My implementation of Arnault's method can be found on GitHub.

---

1. I had first seen this idea on AMPS18, which was the main inspiration for my thesis. ↩︎