*Article* **On the Number of Witnesses in the Miller–Rabin Primality Test**

## **Shamil Talgatovich Ishmukhametov \*, Bulat Gazinurovich Mubarakov and Ramilya Gakilevna Rubtsova**

Institute of Computational Mathematics and Information Technology, Kazan Federal University, Kremlevskya St. 35, Kazan 420008, Russia; mubbulat@mail.ru (B.G.M.); Ramilya.Rubtsova@kpfu.ru (R.G.R.)

**\*** Correspondence: Shamil.Ishmukhametov@kpfu.ru

Received: 24 February 2020; Accepted: 28 March 2020; Published: 1 June 2020

**Abstract:** In this paper, we investigate the popular Miller–Rabin primality test and study its effectiveness. The ability of the test to determine prime integers is based on the difference of the number of primality witnesses for composite and prime integers. Let *W*(*n*) denote the set of all primality witnesses for odd *n*. By Rabin's theorem, if *n* is prime, then each positive integer *a* < *n* is a primality witness for *n*. For composite *n*, the power of *W*(*n*) is less than or equal to *ϕ*(*n*)/4 where *ϕ*(*n*) is Euler's Totient function. We derive new exact formulas for the power of *W*(*n*) depending on the number of factors of tested integers. In addition, we study the average probability of errors in the Miller–Rabin test and show that it decreases when the length of tested integers increases. This allows us to reduce estimations for the probability of the Miller–Rabin test errors and increase its efficiency.

**Keywords:** prime numbers; primality test; Miller–Rabin primality test; strong pseudoprimes; primality witnesses

## **1. Introduction**

The MillerRabin primality test is an algorithm that checks whether a given number is prime or composite. Its original version, due to Gary L. Miller, was deterministic and relied on the unproved extended Riemann Hypothesis [1]. Michael O. Rabin modified it to obtain a probabilistic algorithm [2].

**Definition 1.** *Let <sup>m</sup> be a positive integer represented as <sup>m</sup>* = <sup>2</sup>*<sup>s</sup>* · *<sup>u</sup> where <sup>u</sup> is odd. We introduce two auxiliary functions bin*(*m*) = *s and odd*(*m*) = *u.*

**Definition 2.** *Let n be an odd natural, n* > 9*. An integer a*, 1 ≤ *a* < *n*, *is called a primality witness for n if it is co-prime to n and one of the following conditions holds:*

$$\begin{aligned} 1. \,\, a^{\text{odd}(n-1)} &\equiv 1 \text{ mod } n, \\ 2. \,\, a^{\text{odd}(n-1)2^i} &\equiv -1 \text{ mod } n \,\,\, for \,\, some \,\, i, \,\, 0 \le i < \,\, b \,\, in(n-1), \end{aligned} \tag{1}$$

(We replaced original Rabin's definition of the compositeness witnesses by the opposite relation). For generality, we count 1 and *n* − 1 as primality witnesses and call them trivial witnesses since they satisfy (1) for any *n*.

Let *W*(*n*) denote the set of all primality witnesses for *n*. The Rabin theorem [2] asserts that if number *n* is prime then each non-zero integer, *a* < *n* is a primality witness for *n*, and therefore, the number of all witnesses |*W*(*n*)| = *n* − 1. For composite *n*, it satisfies inequality |*W*(*n*)| ≤ *ϕ*(*n*)/4 where *ϕ*(*n*) is Euler's totient function. Since Rabin did not consider 1 as a witness, then he stated the strict inequality |*W*(*n*)| < *ϕ*(*n*)/4.

Later, Gary Miller [1] developed a primality test that takes any integer *a*, 1 < *a* < *n*, checks if *a* is not a factor of *n* (otherwise, *n* is trivially composite), and whether *a* is a primality witness for *n*, that is, lies in the set *W*(*n*). If the answer is positive, then *n* is probable prime with probability exceeding 3/4. If we need in a more exact result, we should repeat this procedure several times taking different numbers *a* < *n*.

The researchers refer to this algorithm as to the Miller and Rabin primality test. We abbreviate it to *MR test*.

**Definition 3.** *Parameters a which are used in Miller's algorithm are called* bases*. They are chosen randomly from interval* [1; *n* − 1]*. If, for a given odd integer, n relation (1) holds at a base a, we say, n passes the MR test at base a. Otherwise, we call a a compositeness witness for n and deduce that n is certainly composite.*

The probability of error after *k* successful iterations becomes less than 1/4*k*. The only type of error in the Rabin' procedure is defining a composite integer as prime.

More details on the Miller–Rabin test can be found in Chapter 3 of text-book [3] by Crandall and Pomerance. We abbreviate Miller–Rabin test as MR test.

**Definition 4.** *Composite integers qualifying by MR test as probable prime at a base a are called strong pseudoprimes relative to base a. Composite integers being probably prime relative to all a from a set A of bases are called strong probable prime relative to set of bases A.*

Investigation of pseudoprime integers has a long history in the Computational Number Theory. We outline main advantages in this direction in the next section.

#### **2. Some History Remarks**

Fist attempts to find fast primality algorithms were based on Fermat's Little Theorem asserting that for prime *n* and for any positive integer *a*, the following relation holds

$$a^n \equiv a \bmod n \tag{2}$$

Indeed, many composite integers do not satisfy (2) and can be discarded after the first check. Composite *n* that satisfy (2) are called Fermat pseudoprimes relative to base *a*.

It is important to note that all strong pseudoprimes relative to a base *a* are also Fermat pseudoprimes relative to *a*.

We can decrease the number of false decisions by Fermat's test by checking the relation (2) with several different *a*. However, this does not allow us to completely avoid false conclusions since so-called Carmichael numbers exist.

Integer *n* is called a Carmichael number if it satisfies (2) for all *a*. Carmichael numbers appear relatively rarely and the least Carmichael number is 561 = 3 · 7 · 11. It is known that Carmichael numbers are exactly those integers which satisfy Korselt's criterion:

Korselt Criterion (1899). A positive compositeinteger *n* is a Carmichael number if and only if *n* is square-free, and for all prime divisors *p* of *n*, it is true that *p* − 1|*n* − 1.

One of the interesting problems is to find for a given odd integer *n* the least witness. In 1994 Alford, Granville and Pomerance proved [4] that such witnesses exceed (log *n*)1/(3 log log log *<sup>n</sup>*) for infinitely many n. We also show that there are finite sets of odd composites which do not have a reliable witness, namely a common witness for all of the numbers in the set.

MR test discards a Carmichael number *n*, if the base was chosen from [1; *n* − 1]\*W*(*n*).

Let us fix a base *a* and let *na* be a least composite integer that the MR Test accepts at the base *a*. Then, any odd *n* < *na* for which *a* is a primality witness, is definitely prime. This means that when we know *na*, we can definitely check any *n* < *na* for primality using only one round of the MR procedure. The corresponding integer *na* is small. But if we take a set *A* of several different bases *a* and find a

least composite *nA* for which all *a* ∈ *A* are primality witness, this *nA* can be very large. Candidates for bases *a* can be any positive integers that are not squares. However, historically, candidates for special bases are chosen from the set of primes.

Let *Pk* denote the set of the first *k* primes *Pk* = {2, 3, 5, 7, ... , *pk*}, and let *ψ<sup>k</sup>* be a least strong pseudoprime relative to *Pk* for a *k* ≥ 1. Function *ψ<sup>k</sup>* is well defined and is exponentially computable. Its computation began already 40 years ago.

First four values of *ψ<sup>k</sup>* have been found by C. Pomerance, J. Selfridge, and S.Waggstaff [5] in 1980.

A systematic calculation of *ψ<sup>k</sup>* for larger *k* has been initiated by J. Jaeschke [6] who elaborated basic algorithms helpful for searching for strong pseudoprimes of different forms. In 1993 Jaeschke calculated *ψ<sup>k</sup>* for 5 ≤ *k* ≤ 8 and proposed upper bounds for *ψ<sup>k</sup>* at 9 ≤ *k* ≤ 11.

F. Arnault in papers [7,8] described another algorithm to search for Carmichael numbers and strong pseudoprimes integers.

Jaeschke' hypothesis have been improved in 2001 by Z. Zang [9] who constructed a lesser 19-digits decimal integer *Q*<sup>11</sup> = 3825123056546413051 bounding above *ψ*11. Z.Zang conjectures that values *ψ<sup>k</sup>* for 9 ≤ *k* ≤ 11 are equal to each other and coincide with *Q*11.

In 2012 J. Jiang and Y. Deng [10] confirmed Zang's Hypothesis by showing that *Q*<sup>11</sup> = *ψ*<sup>9</sup> = *ψ*<sup>10</sup> = *ψ*11.

The last record is reached by J. Sorenson and J. Webster [11] in 2016 . They found *ψ*<sup>12</sup> and *ψ*13, where *<sup>ψ</sup>*<sup>13</sup> = <sup>3317044064679887385961981</sup> ≈ 3.3 · <sup>10</sup>24. So at the moment we can successfully determine prime integers less than 3.3 · 1024 by only 13 rounds of the MR test. But this bound is much less than integers used in Cryptography. For example, DSS algorithm uses prime integers of length 256 bits (≈80 decimal digits).

Another branch of investigations in connected with the problem of distribution of Fermat pseudoprimes and strong pseudoprimes. Let *F*(*n*) denote set

$$F(n) = \{ a \bmod n \; : \; a^{n-1} \equiv 1 \bmod n \} \dots$$

Clearly, *F*(*n*) ⊇ *W*(*n*).

In 1985 P. Erdos and C. Pomerance [12] studied an asymptotic behavior of average function

$$A(x) = \frac{1}{x} \sum\_{n \le x} \prime |F(n)|$$

where sum is counted over odd integers. They showed using complex number-theoretical calculations that *A*(*x*) is a growing function bounded below by *x*15/23.

Our average function *Avg*(*x*) looks close to *A*(*x*) but we show that for almost all composite *n W*(*n*) consists of only two elements 1 and *n* − 1 and function *Avg*(*x*) tends to zero with *x* tending to infinity.

Average number of errors in the MR test was also studied in 1993 by I. Damgard, P. Landrock and C Pomerance. In paper [13] they studied an average probability of the false decision by the MR test in the following procedure:

Fix *k* > 0 and *t* > 0 and choose randomly *k*-bit odd integer *n*. Check it with *t* rounds of MR test with randomly chosen bases from [1; *n* − 1]. If *n* was discarded during the procedure (that is, found *a* ∈ *W*(*n*)), take another *n*. Continue until *n* was found passed *t* rounds. Let *pk*,*<sup>t</sup>* be the probability that the procedure returns a composite integer.

The authors found explicit upper bounds for various *k* and *t*. In particular they proved that *pk*,1 ≤ *<sup>k</sup>*242<sup>−</sup> <sup>√</sup>*<sup>k</sup>* for *<sup>k</sup>* <sup>≥</sup> 2. Their results show that the probability of false decisions of the MR test depends on the length of tested numbers and it decreases if the length of the numbers increases.

#### **3. Counting Number of Witnesses**

In this section we deduce exact formulas for the number of primality witnesses for different types of composite integers.

We begin our investigation with a little proposition improving Rabin's estimate.

**Theorem 1.** *If a* ∈ *W*(*n*)*, then n* − *a* ∈ *W*(*n*)*.*

**Proof.** Let *<sup>k</sup>* <sup>=</sup> *ordn*(*a*). If *<sup>k</sup>* is odd, then *<sup>a</sup>odd*(*n*−1) mod *<sup>n</sup>* <sup>=</sup> 1, and (*<sup>n</sup>* <sup>−</sup> *<sup>a</sup>*)*odd*(*n*−1) ≡ −<sup>1</sup> mod *<sup>n</sup>*, therefore, *n* − *a* is also a witness.

If *<sup>k</sup>* is even, then *<sup>a</sup>k*/2 ≡ −<sup>1</sup> mod *<sup>n</sup>*. If *<sup>k</sup>*/2 is even, then (*<sup>n</sup>* − *<sup>a</sup>*)*k*/2 ≡ *<sup>a</sup>k*/2 ≡ −<sup>1</sup> mod *<sup>n</sup>*, and (*n* − *a*) is a witness.

Finally, if *<sup>k</sup>*/2 is odd, then (*<sup>n</sup>* <sup>−</sup> *<sup>a</sup>*)*k*/2 ≡ −*ak*/2 <sup>≡</sup> <sup>1</sup> mod *<sup>n</sup>*. Since *<sup>k</sup>*/2 <sup>|</sup> *odd*(*<sup>n</sup>* <sup>−</sup> <sup>1</sup>), then *<sup>a</sup>odd*(*n*−1) <sup>≡</sup> 1 mod *n*, and (*n* − *a*) again is a witness.

This completes the proof.

**Corollary 1.** *(The Improved Rabin Theorem). Let n be a natural, and A be an arbitrary set of bases less than n, co-prime to n, such that for any a* ∈ *A, n* − *a is not in A. If all bases a* ∈ *A are primality witnesses of n, then n is probable prime with probability of error less than or equal to* 1/16*k.*

Indeed, when we found a primality witness *a* for integer *n*, we get two primality witnesses for *n*, namely, *<sup>a</sup>* and *<sup>n</sup>* − *<sup>a</sup>*. So, this reduces the probability of error by a factor of 42 = 16.

Let *Nw*(*n*) = |*W*(*n*)| be the power of number of primality witnesses *W*(*n*). As mentioned earlier, for prime *n Nw*(*n*) = *n* − 1, and for composite *n Nw* ≤ *ϕ*(*n*)/4.

Below we estimate function *Nw*(*n*) more exactly. First we formulate a theorem restricting possible witnesses for a composite *n*.

**Theorem 2.** *Let n* = *u* · *v for co-prime factors u and v (possibly, composite), and a* ∈ *W*(*n*)*. Then,*

$$\begin{array}{l} 1.\,\,ord\_{\mathbb{U}}(a) \mid \,\,\,GCD(\varphi(u),\,(u-\varrho(u))v-1),\\ 2.\,\,ord\_{\mathbb{U}}(a) \mid \,\,\,GCD(\varphi(v),\,(v-\varrho(v))u-1),\\ 3.\,\,bin(ord\_{\mathbb{U}}(a)) = bin(ord\_{\mathbb{V}}(b)). \end{array} \tag{3}$$

**Proof.** 1. Since *<sup>a</sup>* is a primality witness for *<sup>n</sup>* then *<sup>a</sup>n*−<sup>1</sup> ≡ <sup>1</sup> mod *<sup>n</sup>* and *<sup>a</sup>n*−<sup>1</sup> ≡ <sup>1</sup> mod *<sup>u</sup>*. Besides, *n* − 1 = *uv* − 1 = *ϕ*(*u*)*v* + (*u* − *ϕ*(*u*))*v* − 1, so

$$1 \equiv a^{\nu - 1} \equiv a^{\varrho(\mu)\upsilon + (\mu - \varrho(\mu))\upsilon - 1} \equiv a^{(\mu - \varrho(\mu))\upsilon - 1} \mod \mu\_\prime$$

since *<sup>a</sup>ϕ*(*u*) <sup>≡</sup> 1 mod *<sup>u</sup>* by Euler's Theorem.

2. By symmetry.

3. If *ordu*(*a*) is odd, then *<sup>a</sup>odd*(*n*−1) <sup>≡</sup> <sup>1</sup> mod *<sup>n</sup>* (otherwise, *<sup>a</sup>* satisfies the second clause of the MRT, and *ordu*(*a*) should be even). Then *<sup>a</sup>odd*(*n*−1) <sup>≡</sup> 1 mod *<sup>v</sup>* and *ordv*(*a*) is odd.

If *bin*(*ordu*(*a*)) = *i* for 0 < *i* < *bin*(*n* − 1), then *a* is a witness by second clause of the MRT, so *aodd*(*n*−1)2*i*−<sup>1</sup> ≡ −<sup>1</sup> mod *<sup>n</sup>*, *<sup>a</sup>odd*(*n*−1)2*i*−<sup>1</sup> ≡ −<sup>1</sup> mod *<sup>v</sup>*, and *<sup>a</sup>odd*(*n*−1)2*<sup>i</sup>* ≡ 1 mod *v*, so *ordv*(*a*) = *odd*(*<sup>n</sup>* − <sup>1</sup>)2*<sup>i</sup>* and *bin*(*ordv*(*a*)) is equal to *<sup>i</sup>*.

The theorem is proved.

**Example 1.** *Let n* = 15 · 19 = 285, *and a* ∈ *W*(*n*)*. By Theorem 2:*

1. *ordu*(*a*) | *GCD*(*ϕ*(*u*), (*u* − *ϕ*(*u*))*v* − 1) = *GCD*(8, 132) = 4, 2. *ordv*(*a*) | *GCD*(*ϕ*(*v*), (*v* − *ϕ*(*v*))*u* − 1) = *GCD*(18, 14) = 2, 3. *bin*(*ordu*(*a*)) = *bin*(*ordv*(*b*)).

*So, possible a satisfies* (*ordu*(*a*), *ordv*(*a*)) = (1, 1)*, or,* (*ordu*(*a*), *ordv*(*a*)) = (2, 2)*, so n* = 285 *has only trivial witnesses* 1 *and n* − 1*.*

**Theorem 3.** *Let n* = *<sup>p</sup><sup>k</sup> be a degree of prime p, then Nw*(*n*) = *<sup>p</sup>* − <sup>1</sup>*.*

**Proof.** Let *<sup>a</sup>* be a witness for *<sup>n</sup>* = *<sup>p</sup>k*, then *orda*(*n*) | *GCD*(*ϕ*(*n*), *<sup>n</sup>* − <sup>1</sup>) = *GCD*(*pk*−1(*<sup>p</sup>* − <sup>1</sup>), *<sup>p</sup><sup>k</sup>* − <sup>1</sup>) = *p* − 1.

Besides, any *a* satisfying *ap*−<sup>1</sup> mod *n* = 1 is a witness of *n*. Indeed, let *ap*−<sup>1</sup> mod *n* = 1. Then, *<sup>m</sup>* = *ordn*(*a*) is a factor of *<sup>n</sup>* − <sup>1</sup> = *<sup>p</sup><sup>k</sup>* − 1. Let *<sup>n</sup>* − <sup>1</sup> = <sup>2</sup>*<sup>s</sup>* · *<sup>t</sup>* for odd *<sup>t</sup>*, therefore, *<sup>m</sup>* = <sup>2</sup>*s*<sup>1</sup> · *<sup>t</sup>*1, where *s*<sup>1</sup> ≤ *s* and *t*<sup>1</sup> is a factor of *t*.

If *s*<sup>1</sup> = 0, then *at*<sup>1</sup> mod *n* = 1, *a<sup>t</sup>* mod *n* = 1 and *a* is a witness by the first clause of the MRT. Otherwise, let 0 <sup>≤</sup> *<sup>r</sup>* <sup>≤</sup> *<sup>s</sup>*<sup>1</sup> be such that *<sup>a</sup>t*12*<sup>r</sup>* ≡ −<sup>1</sup> mod *<sup>n</sup>*. Then *<sup>a</sup>t*2*<sup>r</sup>* ≡ −1 mod *n* and *a* is a witness by the second clause of the MRT. This completes the proof.

We call integer *n semiprime* if it is a product of two distinct primes *n* = *pq*, *p* < *q*. Semiprimes are close to primes, and we prove below that they have a maximal number of primality witnesses among composite numbers.

**Theorem 4.** *Number of witnesses of semiprime n* = *pq is equal to*

$$N\_w(pq) = (odd(d))^2 \cdot (4^{bin(d)} + 2) / 3,\tag{4}$$

*where d* = *GCD*(*p* − 1, *q* − 1)*.*

We begin with example of application of this formula.

**Example 2.** *Let <sup>n</sup>* = <sup>11</sup> · <sup>31</sup> = <sup>341</sup>*. Then <sup>d</sup>* = *GCD*(*<sup>p</sup>* − 1, *<sup>q</sup>* − <sup>1</sup>) = <sup>10</sup> = <sup>5</sup> · <sup>2</sup>1*, odd*(*d*) = <sup>5</sup>*, <sup>s</sup>* = *bin*(*d*) = 1*. By the theorem,*

$$N\_w(31) = 5^2 \cdot (4+2)/3 = 50.5$$

**Proof.** Let *d* = *GCD*(*p* − 1, *q* − 1). Applying Theorem 2 to *n* = *pq* we obtain

$$1. \operatorname{ord}\_p(a) |d, \operatorname{ord}\_q(a) |d, \\\ 2. \operatorname{bin}(\operatorname{ord}\_\mathbb{N}(a)) = \operatorname{bin}(\operatorname{ord}\_\mathbb{U}(b)).$$

We distribute all *n*-witnesses *a* into *s* + 1 classes *Wi*, 0 ≤ *i* ≤ *s*, where class *Wi* consists of *a* with *bin*(*ordp*(*a*)) = *bin*(*ordq*(*a*)) = *i*.

Class *W*<sup>0</sup> contains such *a* that both *ordp*(*a*) and *ordq*(*a*) are odd. Let *a* ∈ *W*0, and (*i*, *j*) = (*ordp*(*a*), *ordq*(*a*)). Numbers *i* and *j* are factors of *u* = *odd*(*d*) by the choice of *a*. Conversely, each integer *a* < *n* satisfying *ordp*(*a*) | *u*, *ordq*(*a*) | *u*, is a witness of *n* and lies in *W*0.

Let fix a pair (*i*, *j*), *i*|*d*, *j*|*d*. By Euler's theorem, in *Zp* there are exactly *ϕ*(*i*) elements of multiplicative order *i*, and in *Zq* there are *ϕ*(*j*) elements of multiplicative order *j*, so, there exist exactly *ϕ*(*i*) · *ϕ*(*j*) pairs (*x*, *y*), 0 < *x* < *p*, 0 < *y* < *q*, such that (*ordp*(*x*), *ordq*(*y*)) = (*i*, *j*). But for each such pair (*x*, *y*) there exists a unique *a* < *n* with (*a* mod *p*, *a* mod *q*)=(*x*, *y*), so there is a injective correspondence between witnesses *a* of *n* with odd orders *ordp*(*a*), *ordq*(*a*), and pairs (*x*, *y*) with *x*|*u*, *y*|*u*. Therefore, the power of *W*<sup>0</sup> is equal to

$$|W\_0| = \sum\_{\mathbf{x}|u,\ \mathbf{y}|u} \varphi(\mathbf{x}) \cdot \varphi(\mathbf{y}) = \left(\sum\_{\mathbf{x}|u} \varphi(\mathbf{x})\right)\left(\sum\_{\mathbf{y}|u} \varphi(\mathbf{y})\right) = u^2,$$

since by a known theorem of Euler for any natural *<sup>m</sup>* <sup>∑</sup>*v*|*<sup>m</sup> <sup>ϕ</sup>*(*v*) = *<sup>m</sup>*.

The next class *W*<sup>1</sup> has the same power *u*<sup>2</sup> since is consists of witnesses *a* with *bin*(*ordp*(*a*)) = *bin*(*ordq*(*a*)) = 1, and

$$|W\_1| = \sum\_{\substack{\mathbf{x} \mid d,\ \mathbf{y} \mid d}} \wp(2\mathbf{x}) \cdot \wp(2\mathbf{y}) = \mu^2,$$

since *ϕ*(2*z*) = *ϕ*(*z*) for odd *z*.

The power of class *Wi* is equal to

$$\sum\_{\substack{\mathbf{x}|d\_r \ y|d}} \varrho(\mathbf{2}^i \mathbf{x}) \cdot \varrho(\mathbf{2}^i y) = 4^{i-1} \boldsymbol{\mu}^2.$$

Therefore, the number of all witnesses *Nw*(*n*) = *<sup>u</sup>*2(<sup>1</sup> + <sup>1</sup> + <sup>4</sup> + ... + <sup>4</sup>*s*−1) = *<sup>u</sup>*<sup>2</sup> · (4*<sup>s</sup>* + <sup>2</sup>)/3. This completes the proof.

**Corollary 2.** *(Rabin's theorem for semiprimes). The number of witnesses of n* = *pq*, *p* ≤ *q*, *is less or equal to ϕ*(*n*)/4*.*

**Proof.** If *p* = *q*, then *Nw*(*n*) = *p* − 1 by Theorem 3, and *ϕ*(*n*)/4 = *p*(*p* − 1)/4, so *Nw*(*n*) < *ϕ*(*n*)/4 at *p* ≥ 5.

Let *p* < *q*. Ratio *Nw*(*n*)/*n* reaches its maximum when *GCD*(*p* − 1; *q* − 1) = *p* − 1, *q* = 2*p* − 1, and *bin*(*p* − 1) = 1. Indeed, *odd*(*n*) is diminishing in two times when *bin*(*p* − 1) is added by 1, and the whole expression in (4) becomes less. Then, *max odd*(*d*)=(*p* − 1)/2, so

$$\max \mathcal{N}\_w(pq) = N\_w(p(2p-1)) = \frac{(p-1)^2}{2} = \frac{\varrho(n)}{4}.$$

**Example 3.** *Let n* = <sup>7</sup> · <sup>13</sup> = <sup>91</sup>*. Nw*(91) = <sup>3</sup><sup>2</sup> · <sup>2</sup> = <sup>18</sup> = *<sup>ϕ</sup>*(91)/4*.*

Now we study function *Nw*(*n*) at products of *k* distinct primes. The general result for such products is formulated below:

**Theorem 5.** *Let n* = *p*<sup>1</sup> · *p*<sup>2</sup> · ... *pk be the product of k distinct primes. Then*

$$N\_{\overline{w}}(n) = u\_1 \cdot u\_2 \cdot \dots \cdot u\_k \cdot \left(1 + \frac{2^{ks} - 1}{2^k - 1}\right), \text{ where}$$

$$s = \min\{\bin(d\_1), \varinjlim(d\_2), \dots, \varinjlim(d\_k)\}, \ d\_i = GCD\left(p\_i - 1; \varinjlim\_{j \neq i} p\_j - 1\right),$$

*ui* = *odd*(*di*)*.*

Let us begin with an example *n* = 7 · 13 · 31 = 2821. The corresponding restrictions are listed below:

$$\begin{array}{l} 1. or d\_p(a) \mid d\_1 = \text{GCD}(p-1; qr - 1) = 6, \; \mu\_1 = 3, \\ 2. or d\_q(a) \mid d\_2 = \text{GCD}(q-1; pr - 1) = 12, \; \mu\_2 = 3, \\ 3. or d\_r(a) \mid d\_3 = \text{GCD}(r-1; pq - 1) = 30, \; \mu\_3 = 15, \\ 4. \; \text{bin}(\text{ord}\_p(a)) = \text{bin}(\text{ord}\_q(b)) = \text{bin}(\text{ord}\_r(b)). \end{array}$$

Since *s* = *min*{*bin*(*d*1), *bin*(*d*2), *bin*(*d*3)} = *min*{1, 2, 1} = 1, we obtain

$$N\_w(2821) = 3 \cdot 3 \cdot 15 \left(1 + \frac{2^3 - 1}{2^3 - 1} \right) = 2700$$

(compare with *ϕ*(*n*)/4 = 6 · 12 · 30/4 = 540).

**Proof.** Let *ui* = *odd*(*di*) and *k*-tuple (*x*1, *x*2, ... , *xk*) contains components *xi* | *ui*, 1 ≤ *i* ≤ *k*. There are *ϕ*(*x*1) · ... · *ϕ*(*xk*) witnesses of *n* with *ordpi* (*a*) = *xi* for 1 ≤ *i* ≤ *k*. So,

$$|\mathcal{W}\_0| = \sum\_{(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_k), \mathbf{x}\_i \mid \boldsymbol{\mu}\_i} \varphi(\mathbf{x}\_1) \cdot \dots \cdot \varphi(\mathbf{x}\_k) =$$

$$= \left(\sum\_{\mathbf{x} \mid \boldsymbol{\mu}\_1} \varphi(\mathbf{x})\right) \cdot \left(\sum\_{\mathbf{x} \mid \boldsymbol{\mu}\_2} \varphi(\mathbf{x})\right) \cdot \dots \left(\sum\_{\mathbf{x} \mid \boldsymbol{\mu}\_k} \varphi(\mathbf{x})\right) = \boldsymbol{\mu}\_1 \cdot \boldsymbol{\mu}\_2 \cdot \dots \cdot \boldsymbol{\mu}\_k.$$

As in the previous theorem, the power of class *W*<sup>1</sup> is equal to power of *W*<sup>0</sup> = *u*<sup>1</sup> · *u*<sup>2</sup> · ... · *uk*, while the power of the each further class *Wi*+<sup>1</sup> is equal to the power of the previous one multiplied by *ϕ*(2*k*) = 2*k*−<sup>1</sup> since each additive *ϕ*(2*<sup>i</sup> <sup>x</sup>*1) · ... · *<sup>ϕ</sup>*(2*<sup>i</sup> xk*) in the previous class corresponds to additive *<sup>ϕ</sup>*(2*i*+1*x*1) · ... · *<sup>ϕ</sup>*(2*i*+<sup>1</sup>*xk*) and their ratio *ri* is

$$r\_i = \frac{\varrho(2^{i+1}\mathbf{x}\_1) \cdot \dots \cdot \varrho(2^{i+1}\mathbf{x}\_k)}{\varrho(2^i\mathbf{x}\_1) \cdot \dots \cdot \varrho(2^i\mathbf{x}\_k)} = 2^k.$$

The proof is complete.

#### **4. Frequency Function**

In this part we introduce a notion of *frequency function* that characterizes the probability to find at one attempt a primality witness for a given integer *n*.

Let define frequency function *Fr*(*n*) as follows

$$Fr(n) = \frac{N\_w(n)}{\varphi(n)}.$$

According to Rabin's theorem, *Fr*(*n*) = 1 for prime *n*, and *Fr*(*n*) ≤ 1/4 for composite *n*. We study distribution of values *Fr*(*n*) for semiprime integers *n* = *pq*, *p* < *q*.

1. We begin our research with case *q* − 1 = *k*(*p* − 1) for *k* ≥ 2. Numbers of this type appear frequently among strong pseudoprimes. Let rewrite *p* and *q* in form *p* = 2*su* + 1, *q* = 2*sku* + 1, where *u* is odd, *s* ≥ 1, and consider different *s*:

\*\*Case 1.\*\*  $s = 1$ ,  $u = 
od d(d) = (p - 1)/2$ ,  $N\_w(pq) = 2u^2 = (p - 1)^2/2$ , 
$$Fr(n) = 
\frac{(p - 1)^2/2}{(p - 1)(q - 1)} = 
\frac{2u^2}{2u 
\cdot 2ku} = 
\frac{1}{2k}$$
.

*Function Fr*(*n*) *reaches its maximum* 1/4 *at k* = 2*:* (*p*, *q*)=(2*u* + 1, 4*u* + 1)*. Since, both p and q are prime then u* ≡ 0 mod 3*, so* (*p*, *q*)=(6*t* + 1, 12*t* + 1)*, t* ≥ 1*. Such pairs form a sequence*

(7, 13), (19, 37),(31, 61), (37, 73),....

**Case 2.** *<sup>s</sup>* = 2, *<sup>u</sup>* = *odd*(*d*)=(*<sup>p</sup>* − <sup>1</sup>)/4, *Nw*(*pq*) = <sup>6</sup>*u*2, *and*

$$Fr(n) = \frac{6\mu^2}{(p-1)(q-1)} = \frac{6\mu^2}{4\mu \cdot 4ku} = \frac{3}{8k}.$$

*Maximum of Fr*(*n*) *is now* 3/16 = 0.1875 *at k* = 2*.*

**Case 3.** *s* ≥ 1*, At arbitrary s we have*

$$Fr(n) = \frac{(1 + (4^s - 1)/3)\mu^2}{(p - 1)(q - 1)} = \frac{(1 + (4^s - 1)/3)\mu^2}{2^s \mu \cdot 2^s k \mu} = \frac{1}{3ku^2 \cdot 2^{2s - 1}} + \frac{1}{3k}.$$

*Thus, function Fr*(*n*) *at semiprimes n* = *pq, q* − 1 = *k*(*p* − 1), *is located in the interval*

$$\frac{1}{3k} < Fr(n) \le \frac{1}{2k'} \text{ } k \ge 2.\tag{5}$$

2. Now, we turn to a common case *n* = *pq*:

$$p = 1 + k\_1 \mu\_\prime \, q = 1 + k\_2 \mu\_\prime \, \text{GCD}(k\_1, k\_2) = 1, \,\,\mu = \text{t2}^s, \,\,\text{t} \,\, \text{odd}.$$

For such *n*

$$N\_w(n) = t^2(4^s + 2) / 3,\ \phi(n) = k\_1 k\_2 t^2 4^s,\ Fr(n) = \frac{4^s + 2}{3k\_1 k\_2 \cdot 4^s}.$$

So,

$$\frac{1}{3k\_1k\_2} < Fr(n) \le \frac{1}{2k\_1k\_2}$$

Conclusion. Function *Fr*(*n*) at semiprimes *n* = *pq* depends mostly on values *k*<sup>1</sup> and *k*<sup>2</sup> in representation *p* = *k*1*u* + 1, *q* = *k*2*u* + 1. *Fr*(*n*) takes maximal values close to 1/4 only at small *k*<sup>1</sup> and *k*2. This completely corresponds to experimental data. Among values *ψ<sup>k</sup>* the most expected are pseudoprimes of form *u* = (*u* + 1)(2*u* + 1) with minimal values *k*<sup>1</sup> = 1 and *k*<sup>2</sup> = 2.

An important question connecting with efficiency of MRT is the average frequency of witnesses for composite numbers. As earlier, we study this problem for semiprime integers.

Let fix any prime *p* and a board *B*. We count average frequency of integers *pq*, *q* > *p*, *pq* ≤ *B*. For convenience, we assume that *B* = *p*(*p* + (*p* − 1)*k*) for a positive *k* ∈ **Z**.

For simplicity we explain all deductions at example *p* = 11. Every prime *q* has *d* = *GCD*(*p* − 1, *q* − 1) equal either 2, or 10.

Let *d* = 10. Corresponding *q* lie in the set {21, 31, 41, 51, 61, 71, 81, 91, 101, ... , 10*k* + 11}, where 10*k* + 11 = *B*/*p*. Each third integer in the sequence is a multiple of 3, some others are multiples of 7, 11 etc. Since *q* should be prime we need to remove them from the sequence. The rest consists of integers

$$Q\_B = \{31, 41, 61, 71, 101, 113 \dots\}. \tag{6}$$

We assume that primes *q* ∈ *QB* are distributed uniformly in the interval [1, *B*/*p*]. Then the average frequency can be estimated as

$$Avg(Fr(n)) \approx \frac{1}{k} \left(\frac{1}{4} + \frac{1}{6} + \dots \, \, \, \, \, \frac{1}{2k}\right) = \frac{1}{2k} \left(1 + \frac{1}{2} + \frac{1}{3} + \dots \, \, \, \, \frac{1}{k}\right).$$

(we remind that *Fr*(*p*(*i*(*p* − 1) + *p*) = 1/2(*i* + 1)).

The expression in the last brackets is a partial sum of the Harmonic Series. Its value is

$$\sum\_{i=1}^{k} \frac{1}{i} < \sum\_{i=1}^{k+1} \frac{1}{i} = \ln k + \gamma + \varepsilon\_{\nu\nu}$$

where *<sup>γ</sup>* = 0.5772... is the Euler—Mascheroni constant and lim*k*→<sup>∞</sup> *<sup>ε</sup><sup>n</sup>* = 0. Constant *<sup>γ</sup>* and additive *<sup>ε</sup><sup>n</sup>* can be ignored so

$$Avg(Fr(n)) < \frac{\ln k}{2k}$$

Since (*<sup>p</sup>* − <sup>1</sup>)*<sup>k</sup>* + <sup>1</sup> = *<sup>B</sup>*/*p*, then *<sup>k</sup>* > *<sup>B</sup>*/*p*<sup>2</sup> − 1 and ln *<sup>k</sup>* < ln *<sup>B</sup>*, so

$$Avg(Fr(n)) < \frac{\ln B}{2(B - p^2)} \cdot p^2 \tag{7}$$

Let us move now to primes *q* of type *d* = *GCD*(*p* − 1, *q* − 1) = 2. They lie in the sequence

$$q \in \{13, 15, 17, 19, 23, 25, 27, 29, \dots, 2m+1\}$$

where 2*m* + 1 = *B*/*p*, *q* = 2*i* + 1, *GCD*(*i*, 5) = 1. When we remove composite integers, the rest contains at least half members.

Integers *n* = *pq* with *GCD*(*p* − 1, *q* − 1) = 2 have only trivial witnesses 1 and *n* − 1 so their frequency function takes values

$$Fr(n) = \frac{2}{(p-1)(q-1)}$$

.

Assuming that such *n* are distributed uniformly in the interval [*p*2; *B*] we estimate the average frequency by expression

$$Avg(Fr) \approx \left(2\sum\_{p \le k \le m} \frac{1}{(p-1)(2k+1)}\right) / \left(\frac{2m+1-p}{2}\right) <$$

$$\frac{4}{(2m+1-p)(p-1)} \cdot \frac{1}{2} \cdot \sum\_{i=(p+1)/2}^{m} \frac{1}{i} < \frac{2}{(2m+1-p)(p-1)} \cdot \ln m$$

Substituting in the last expression 2*m* + 1 = *B*/*p* we get

$$Avg(Fr) < \frac{2p \ln B}{(B - p^2)(p - 1)}\tag{8}$$

Expressions (7) and (8) give upper bounds for two types of integers *n* = *pq*. In the second case the estimation is lesser so average estimation for the united class of all *n* = *pq* ≤ *B*, *p* < *q*, can be set by the upper bound of (7). This assertion does not depend on a special *p* = 11 so we can state the following theorem.

**Theorem 6.** *Let p be a prime and B satisfy B* > *p*2*. Then the average frequency of witnesses in the class of semiprimes n* = *pq* ≤ *B*, *q* > *p*, *has an upper bound*

$$Avg(Fr(n)) < \frac{p^2 \ln B}{2(B - p^2)}$$

Note than limit of the average function is 0 as *B* → ∞. This explains the phenomenon that the number of false conclusions in the Miller–Rabin test decreases when length of tested integers increases.

#### **5. Numbers with Maximal Frequency of Witnesses**

In this section we study composite *n* with maximal frequency *Fr*(*n*) = 1/4. Let *n* = *p*<sup>1</sup> *p*<sup>2</sup> ... *pk* be the product of *k* different primes.

We begin with case *k* = 2. As we see from the previous section, integers *n* = *pq* have maximal frequency only in case when *q* = 2*p* − 1. Such pairs appear comparatively often, and their quantity is diminishing together with their size.

Table <sup>1</sup> contains number of semiprimes with maximal frequency in intervals [(*<sup>i</sup>* − <sup>1</sup>) · 105; *<sup>i</sup>* · <sup>10</sup>5; ], 1 ≤ *<sup>i</sup>* < 10.

**Table 1.** Distribution of semiprimes with maximal frequency below 106.


Case *k* = 3 is more interesting. In order function *Fr*(*pqr*) reached its maximum = 0.25, we need satisfaction of four requirements:

$$\begin{array}{l} 1.\,\,GCD(p-1;qr-1) = p-1, \\ 2.\,\,GCD(q-1;pr-1) = q-1, \\ 3.\,\,GCD(r-1;pq-1) = r-1. \\ 4.\,\,b
in(p-1) = b
in(q-1) = 1.\end{array} \tag{9}$$

Such triples exist, and an example of it was already given in Rabin's paper [2] *n* = 487 · 1531 · 2683 = 2000436751. Rabin himself estimated *Fr*(*n*) as 0.2493, but the difference is due to the fact that he did not include 1 in the list of witnesses.

Such triples appear much more seldom and have a form

$$n = (2k\_1 + 1)\mu \cdot (2k\_2 + 1)\mu \cdot (2k\_3 + 1)\mu \text{ for } \mu \in \mathcal{N}.$$

We arranged the search of such triples at a computer and found 160 such integers not exceeding <sup>2</sup> · <sup>10</sup>14. The least triple we found is

$$n = 19 \cdot 199 \cdot 271 = 1024651... $$

The largest found triple has a form *n* = (*u* + 1)(3*u* + 1)(5*u* + 1) at *u* = 24102:

$$n = 24103 \cdot 72307 \cdot 120511 = 21002 \text{ 84533 02331}.$$

Let us study the form *u*, 3*u*, 5*u* and find restrictions on *u* in order to *n* = (*u* + 1)(3*u* + 1)(5*u* + 1) satisfies first 3 conditions of (9). The first requirement is satisfied automatically. The second and third requirement are listed below:

$$(3u+1)-1 \mid (u+1)(5u+1)-1 \to u \equiv 0 \bmod 3.$$

$$(5u+1)-1 \mid (u+1)(3u+1)-1 \to 3u+4 \equiv 0 \bmod 5.$$

so *u* = 6 + 15*t* for *t* ≥ 1. If we add requirements *p* ≡ *q* ≡ *r* ≡ 3 mod 4 we obtain

$$115t + 7 \equiv 3 \bmod 4 \to t \equiv 1 \bmod 4,\ u = 6 + 15(1 + 4t\_1) = 21 + 60t\_1.$$

Let now consider products of *k* primes where *k* ≥ 4. The maximum of frequency of such products is 1/2*k*−1, since it is reached when for any *<sup>i</sup>* <sup>≤</sup> *<sup>k</sup>* (*pi* <sup>−</sup> <sup>1</sup>)/2 is odd, and (*pi* <sup>−</sup> <sup>1</sup>) <sup>|</sup> <sup>∏</sup>(*pj*=*<sup>i</sup>* <sup>−</sup> <sup>1</sup>). Then,

$$Fr(p) = 2 \cdot \prod\_{i=1}^{k} \frac{p\_i - 1}{2} = \frac{\wp(n)}{2^{k-1}}.$$

A quick search of tuples *n* = *pqrt* below 10<sup>12</sup> gave 70 examples of them. The least 4-tuple was

$$n = 19 \cdot 31 \cdot 127 \cdot 547 = 40917241\_{\prime}$$

while the largest was

$$n = 19 \cdot 127 \cdot 14071 \cdot 29347 = 99 \, 64281 \, 70081\ldots$$

Some computational results on distribution of strong semiprime integers can be found in [14].

## **6. Conclusions**

In this section we will summarize the main results of the paper.


Since such integers have maximal values of *F*(*n*) among all composites, this opens a way in future investigations to find exact upper bounds for average values of frequency function among all *k*-bit odd integers for any *k*.

4. Finally, we described possible forms of composites with maximal values of frequency function for products of *k* distinct primes at *k* ≥ 2 and using computer calculations found their examples and their quantity at initial intervals of set of all naturals.

**Author Contributions:** S.T.I. gave impetus to the research and proved Theorems 1 and 2. B.G.M. proved other theorems and propositions, and R.G.R. developed software for testing results. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by RFBR grant 18-47-16005. This investigation was supported by the grant of Scientific and Educational Mathematical Center of the Volga Federal District, agreement No. 075-02-2020-1478.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
