# A bevy of ones

Take any positive integer d that is not a multiple of 2 or 5. Then there is some integer k such that d × k has only 1’s in its decimal representation. For example, take d = 13. We have

13 × 8457 = 111111.

Or if we take d = 27,

27 × 4115226337448559670781893 = 111111111111111111111111111.

Let’s change our perspective and start with the string of 1’s. If d is not a multiple of 2 or 5, then there is some number made up of only 1’s that is divisible by d. And in fact, the number of 1’s is no more than d.

This theorem generalizes to any integer base b > 1. If d is relatively prime to b, then there is a base b number with d or fewer “digits” which is divisible by d [1].

The following Python code allows us to find the number k such that d × k has only 1’s in its base b representation, provided k is relatively prime to b. It returns the number of 1’s we need to string together to find a multiple of k. If k shares a factor with b, the code returns 0 because no string of 1’s will ever be divisible by k.

    from math import gcd

def all_ones(n, b = 10):
return sum(b**i for i in range(n))

def find_multiple(k, b = 10):
if gcd(k, b) > 1:
return 0
for n in range(1, k+1):
if all_ones(n, b) % k == 0:
return n


The two Python functions above default to base 10 if a base isn’t provided.

We could find a multiple of 5 whose hexadecimal representation is all 1’s by calling

    print(find_multiple(5, 16))

and this tells us that 11111sixteen is a multiple of 5, and in fact

5 × 369dsixteen = 11111sixteen.

[1] Elmer K. Hayashi. Factoring Integers Whose Digits Are all Ones. Mathematics Magazine, Vol. 49, No. 1 (Jan., 1976), pp. 19-22

# Binomial coefficients mod primes

Imagine seeing the following calculation:

The correct result is

and so the first calculation is off by 25 orders of magnitude.

But there’s a variation on the calculation above that is correct! A theorem by Édouard Lucas from 1872 that says for p prime and for any nonnegative integers m and n,

So while the initial calculation was grossly wrong as stated, it is perfectly correct mod 19. If you divide 487343696971437395556698010 by 19 you’ll get a remainder of 10.

A stronger versions of Lucas’ theorem [1] says that if p is at least 5, then you can replace mod p with mod p³. This is a stronger conclusion because it says not only is the difference between the left and right side of the congruence divisible by p, it’s also divisible by p² and p³.

In our example, not only is the remainder 10 when 487343696971437395556698010 is divided by 19, the remainder is also 10 when dividing by 19² = 361 and 19³ = 6859.

## More on binomial coefficients

[1] V. Brun, J. O. Stubban, J. E. Fjeldstad, L. Tambs, K. E. Aubert, W. Ljunggren, E. Jacobsthal. On the divisibility of the difference between two binomial coefficients, Den 11te Skandinaviske Matematikerkongress, Trondheim, 1949, 42–54.

# Bit flipping to primes

Someone asked an interesting question on MathOverflow: given an odd number, can you always flip a bit in its binary representation to make it prime?

It turns out the answer is no, but apparently it is very often the case an odd number is just a bit flip away from being prime. I find that surprising.

Someone pointed out that $$2131099$$ is not a bit flip away from a prime, and that this may be the smallest example [1]. The counterexample $$2131099$$ is itself prime, so you could ask whether an odd number is either a prime or a bit flip away from a prime. Is this always the case? If not, is it often the case?

The MathOverflow question was stated in terms of Hamming distance, counting the number of bits in which two bit sequences differ. It asked whether odd numbers are always Hamming distance 1 away from a prime. My restatement of the question asks whether the Hamming distance is always at most 1, or how often it is no more than 1.

You could ask more generally about the Hamming distance to the nearest prime. Is it bounded, if not by 1, then by another finite number? If so, what is the smallest such bound? What is the probability that its value is 1? Etc.

This ties into a couple of other things I’ve blogged about. A few weeks ago I wrote about new work on the problem of finding the proportion of odd numbers that can be written as the sum of a power of 2 and a prime. That’s a little different problem since bit flipping is taking the XOR (exclusive or) and not always the same as addition. It also leaves out the possibility of flipping a bit beyond the most significant bit of the number, i.e. adding to a number n a power of 2 greater than n.

Another related post is on the Rowhammer attack on public key cryptography. By flipping a bit in the product of two primes, you can produce a number which is much easier to factor.

These two posts suggest a variation on the original problem where we disallow flipping bits higher than the most significant bit of n. So giving a k-bit number n, how often can we flip one of its k bits and produce a prime?

[1] Note that the bit flipped may be higher than the most significant bit of the number, unless ruled out as in the paragraph above. Several people have asked “What about 85?” It is true that flipping any of the seven lowest bits of 85 will not yield a prime. But flipping a zero bit in a more significant position will give a prime. For example, 1024 + 85 is prime. Bur for $$2131099$$ it is not possible to add any larger power of 2 to the number and produce a prime.

# Square roots of Gaussian integers

In previous post I showed how to compute the square root of a complex number. I gave as an example that computed the square root of 5 + 12i to be 3 + 2i.

(Of course complex numbers have two square roots, but for convenience I’ll speak of the output of the algorithm as the square root. The other root is just its negative.)

I chose zxiy in my example so that x² + y² would be a perfect square because this simplified the exposition.That is, I designed the example so that the first step would yield an integer. But I didn’t expect that the next two steps in the algorithm would also yield integers. Does that always happen or did I get lucky?

It does not always happen. My example was based on the Pythagorean triple (5, 12, 13). But (50, 120, 130) is also a Pythagorean triple, and the square root of z = 50 + 130i does not have integer real and imaginary parts.

At this point it will help to introduce a little terminology. A Gaussian integer is a complex number with integer real and imaginary parts. We could summarize the discussion above by saying the square root of 5 + 12i is a Gaussian integer but the square root of 50 + 120i is not.

While (5, 12, 13) and (50, 120, 130) are both are Pythagorean triples, only the first is a primitive Pythagorean triple, one in which the first two components are relatively prime. So maybe we should look at primitive Pythagorean triples.

There is a theorem going all the way back to Euclid that primitive Pythagorean triples have the form

(m² – n², 2mn, m² + n²)

where m and n are relatively prime integers and one of m and n is even.

Using the algorithm in the previous post, we can show that if x and y are the first components of a triple of the form above, then the square root of the Gaussian integer x + iy is also a Gaussian integer. Specifically we have

ℓ = m² + n²
u = √((m² + n² + m² – n²)/2) = m
v = sign(y) √((m² + n² – m² + n²)/2) = sign(y) n

This means u and v are integers and so the square root of x + iy is the Gaussian integer u + iv.

There is one wrinkle in the discussion above. Our statement of Euclid’s theorem left out the assumption that we’ve ordered the sides of our triangle so that the odd side comes first. Primitive Pythagorean triples could also have the form

(2mnm² – n², m² + n²)

and in that case our proof would break down.

Surprisingly, there’s an asymmetry between the real and imaginary parts of the number whose square root we want to compute. It matters that x is odd and y is even. For example, √(5 + 12i) is a Gaussian integer but √(12 + 5i) is not!

So a more careful statement of the theorem above is that if x and y are the first two components of a primitive Pythagorean triple, the square root of xiy is a Gaussian integer.

# How probable is a probable prime?

A probable prime is a number that passes a test that all primes pass and that most composite numbers fail. Specifically, a Fermat probable prime is a number that passes Fermat’s primality test. Fermat’s test is the most commonly used, so that’s nearly always what anyone means by probable prime unless they’re more specific.

A number n is a Fermat probable prime for base b if

bn-1 = 1 (mod n).

This test isn’t conclusive, but it can be implemented efficiently and it weeds out most composite numbers. To read more on probable primes, see this post.

If a number is a probable prime, how probable is it that it really is prime? This post will briefly summarize some results from a paper [1] that makes this precise. From that paper:

… let P(x) denote the probability that an integer n is composite given that

1. n is chosen at random with 1 < nx, n odd,
2. b is chosen at random with 1 < b < n − 1, and
3. n is a probable prime to the base b.

The authors give upper bounds on P(x) for x equal to various powers of 2. In particular they report

P(2250) ≤ 5.876 × 10−6

and

P(2260) ≤ 3.412 × 10−6

and so a the chances that a 256-bit probable prime is composite are in the neighborhood of 4 in a million. In practice, one would test with multiple b‘s. The tests for various b‘s aren’t entirely independent, but running the tests for multiple bases does mean that fewer composite numbers slip through. There are a few pesky numbers, the Carmichael numbers, that are Fermat probable primes for nearly all bases (see footnote [2] here for more details), but these are rare.

I looked through the paper for results for larger powers of 2 to get results that would be applicable to, for example, RSA keys. The largest result explicit in the paper is

P(2330) ≤ 8.713 × 10−8

though I would like to know P(21024) and P(21536) since RSA keys are the products of (probable) primes this large. Presumably the results in [1] could be used to compute P(x) for larger values of x, but I haven’t read the paper closely enough to how much personal effort or computational effort that would require. I imagine it would be difficult or else the authors would have reported results for probable primes of the size frequently used in applications.

## Related posts

[1] Jared Ducker Lichtman and Carl Pomerance. Improved error bounds for the Fermat primality test on random inputs. Mathematics of Computation. Volume 87, Number 314, November 2018, pages 2871–2890.

# Relatively prime determinants

Suppose you fill two n×n matrices with random integers. What is the probability that the determinants of the two matrices are relatively prime? By “random integers” we mean that the integers are chosen from a finite interval, and we take the limit as the size of the interval grows to encompass all integers.

Let Δ(n) be the probability that two random integer matrices of size n have relatively prime determinants. The function Δ(n) is a strictly decreasing function of n.

The value of Δ(1) is known exactly. It is the probability that two random integers are relatively prime, which is well known to be 6/π². I’ve probably blogged about this before.

The limit of Δ(n) as n goes to infinity is known as the Hafner-Sarnak-McCurley constant [1], which has been computed to be 0.3532363719…

Since Δ(n) is a decreasing function, the limit is also a lower bound for all n.

## Python simulation

Here is some Python code to experiment with the math discussed above. We’ll first do a simulation to show that we get close to 6/π² for the proportion of relatively prime pairs of integers. Then we look at random 2×2 determinants.

    from sympy import gcd
from numpy.random import randint
from numpy import pi

def coprime(a, b):
return gcd(a, b) == 1

def random_int(N):
return randint(-N, N)

def random_det(N):
a, b, c, d = randint(-N, N, 4)
return a*d - b*c

count = 0
N = 10000000 # draw integers from [-N, N)
num_reps = 1000000

for _ in range(num_reps):
count += coprime(random_int(N), random_int(N))
print("Simulation: ", count/num_reps)
print("Theory: ", 6*pi**-2)


This code printed

    Simulation:  0.607757
Theory:  0.6079271018540267


when I ran it, so our simulation agreed with theory to three figures, the most you could expect from 106 repetitions.

The analogous code for 2×2 matrices introduces a function random_det.

    def random_det(N):
a, b, c, d = randint(-N, N, 4, dtype=int64)
return a*d - b*c


I specified the dtype because the default is to use (32 bit) int as the type, which lead to Python complaining “RuntimeWarning: overflow encountered in long_scalars”.

I replaced random_int with random_det and reran the code above. This produced 0.452042. The exact value isn’t known in closed form, but we can see that it is between the bounds Δ(1) = 0.6079 and Δ(∞) = 0.3532.

## Theory

In [1] the authors show that

This expression is only known to have a closed form when n = 1.

## Related posts

[1] Hafner, J. L.; Sarnak, P. & McCurley, K. (1993), “Relatively Prime Values of Polynomials”, in Knopp, M. & Seingorn, M. (eds.), A Tribute to Emil Grosswald: Number Theory and Related Analysis, Providence, RI: Amer. Math. Soc., ISBN 0-8218-5155-1.

# Prime plus power of 2

A new article [1] looks at the problem of determining the proportion of odd numbers that can be written as a power of 2 and a prime. A. de Polignac conjectured in 1849 that all odd numbers have this form.

A counterexample is 127, and so apparently the conjecture was that every odd number is either prime or a power of 2 plus a prime [2]. Alternatively, you could say that a prime is a prime plus a power of 2 where the exponent is -∞.

The smallest counterexample to Polignac’s conjecture is 905. It is clearly not prime, nor are any of the values you get by subtracting powers of 2: 393, 649, 777, 841, 873, 889, 897, 901, and 903.

The proportion of numbers of the form 2n plus a prime is known to be between 0.09368 and 0.49095. In [1] the authors give evidence that the proportion is around 0.437.

You could generalize Polignac’s conjecture by asking how many powers of 2 you need to add to a prime in order to represent any odd number. Clearly every odd number x is the sum of some number of powers of 2 since you can write numbers in binary. But is there a finite upper bound k on the number of powers of 2 needed, independent of x? I imagine someone has already asked this question.

Polignac conjectured that k = 1. The example x = 905 above shows that k = 1 won’t work. Would k = 2 work? For example, 905 = 2 + 16 + 887 and 887 is prime. Apparently k = 2 is sufficient for numbers less than 1,000,000,000.

[1] Gianna M. Del Corso, Ilaria Del Corso, Roberto Dvornicich, and Francesco Romani. On computing the density of integers of the form 2n + p. Mathematics of Computation. https://doi.org/10.1090/mcom/3537 Article electronically published on April 28, 2020

# Powers that don’t change the last digit

If you raise any number to the fifth power, the last digit doesn’t change.

Here’s a little Python code to verify this claim.

    >>> [n**5 for n in range(10)]
[0, 1, 32, 243, 1024, 3125, 7776, 16807, 32768, 59049]


In case you’re not familiar with Python, or familiar with Python but not familiar with list comprehensions, the code is very similar to set builder notation in math:

{n5 | n = 0, 1, 2, … 9}

To emphasize the last digits, we could take the remainder by 10:

    >>> [n**5 % 10 for n in range(10)]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]


The reason this works is related to Fermat’s little theorem.

We could generalize this to the last digits in other bases. For example, in base 21, raising numbers to the 13th power doesn’t change the last digit. You could verify this by running

    [x**13 % 21 for x in range(21)]

Why do we take 5th powers in base 10 and 13th powers in base 21? Because φ(10) + 1 = 5 and φ(21) + 1 = 13 where φ(m) is Euler’s “totient” function applied to the base. The function φ(m) is defined as the number of positive integers less than m that are relatively prime to m.

If we’re in base 30, φ(30) = 8, and so we’d expect 9th powers to leave the last digit unchanged. And we can evaluate

    [n**9 % 30 for n in range(30)]

to see that this is the case.

Now let’s try base 16. We have φ(16) = 8, so let’s raise everything to the 9th power.

    >>> [n**9 % 16 for n in range(16)]
[0, 1, 0, 3, 0, 5, 0, 7, 0, 9, 0, 11, 0, 13, 0, 15]


That didn’t work out the way the other examples did. Why’s that? The bases 10, 21, and 30 are the product of distinct primes, but 16 is a prime power.

However, there’s something else interesting going on in base 16. Instead of taking the remainder by 16, i.e. looking at the last hexadecimal digit, let’s look at all the digits.

    >>> [hex(n**9) for n in range(16)]
['0x0', '0x1', '0x200', '0x4ce3', '0x40000', ..., '0x8f3671c8f']


In base 16, the last digit of n9 is not the same as the last digit of n. But the last non-zero digit is the same as the last digit of n.

# Various kinds of pseudoprimes

Every condition that all primes satisfy has a corresponding primality test and corresponding class of pseudoprimes.

The most famous example is Fermat’s little theorem. That theorem says that for all primes p, a certain congruence holds. The corresponding notion of pseudoprime is a composite number that also satisfies Fermat’s congruence.

To make a useful primality test, a criterion should be efficient to compute, and the set of pseudoprimes should be small.

## Fermat pseudoprimes

A Fermat pseudoprime base b is a composite number n such that

bn-1 = 1 (mod n).

The smallest example is 341 because

2340 = 1 mod 341

and 341 = 11*31 is composite. More on Fermat pseudoprimes here.

## Wilson pseudoprimes

Fermat’s little theorem leads to a practical primality test, but not all theorems about primes lead to such a practical test. For example, Wilson’s theorem says that if p is prime, then

(p – 1)! = -1 mod p.

A Wilson pseudoprime would be a composite number n such that

(n – 1)! = -1 mod n.

The good news is that the set of Wilson pseudoprimes is small. In fact, it’s empty!

The full statement of Wilson’s theorem is that

(p – 1)! = -1 mod p

if and only if p is prime and so there are no Wilson pseudoprimes.

The bad news is that computing (p – 1)! is more work than testing whether p is prime.

## Euler pseudoprimes

Let p be an odd prime and let b be a positive integer less than p. Then a theorem of Euler says that

where the symbol on the right is the Legendre symbol. It looks like a fraction in parentheses, but that’s not what it means. The symbol is defined to be 1 if b is a square mod p and -1 otherwise. It’s an unfortunate but well-established bit of notation.

We could turn Euler’s theorem around and convert it into a primality test by saying that if a number p does not the congruence above, it cannot be prime.

There’s one complication with this: the Legendre symbol is only defined if p is prime. However, there is a generalization of the Legendre symbol, the Jacobi symbol, which is defined as long as the top argument is relatively prime to the bottom argument. Euler pseudoprimes are also known as Euler-Jacobi pseudoprimes.

The Jacobi symbol uses the same notation as the Legendre symbol, but there is no ambiguity because if the bottom argument is prime, the Jacobi symbol equal the Legendre symbol. More on Legendre and Jacobi symbols here.

An odd composite number n is a Euler pseudoprime base b if n and b are relatively prime and

where the symbol on the right is now the Jacobi symbol.

There are efficient algorithms for computing Jacobi symbols, and so the criterion above is a practical way to prove a number is not prime. But there are Euler pseudoprimes that slip past this test. A small example is n = 9 and b = 10. Euler pseudoprimes are less common than Fermat pseudoprimes.

## Other examples

There are many other examples: Lucas pseudoprimes, elliptic pseudoprimes, etc. Any theorem about prime numbers can be recast as a primality test, and the composite numbers that slip through labeled pseduoprimes relative to that test. The question is whether such a test is useful, i.e. efficient to compute and without too many false positives.

# Finding large pseudoprimes

Fermat’s little theorem says that if p is a prime number, then for any integer b,

bp-1 = 1 (mod p).

This gives a necessary but not sufficient test for a number to be prime. A number that satisfies the equation above but is not prime is called a pseudoprime [1] for base b. For example, 341 is a pseudoprime base 2 because

2340 = 1 mod 341

and clearly 341 = 11*31.

How might you go about hunting for pseudoprimes? This page gives a way to find pseduoprimes base 2, quoting the following theorem.

If both numbers q and 2q – 1 are primes and n = q(2q-1) then 2n-1 = 1 (mod n) if and only if q is of the form 12k + 1.

Here’s Python code that will find pseudoprimes of a given size in bits.

    from secrets import randbits
from sympy import isprime

def find_pseudoprime(numbits):
n = numbits // 2
while True:
q = randbits(n)
if isprime(q) and q%12 == 1 and isprime(2*q-1):
return q*(2*q-1)

print( find_pseudoprime(256) )


This returned

47362882341088596725068562696893704769436677460225591859092704246296157080253

in under a second.

Clearly the return value of find_pseudoprime is composite since it is computed as the product of two integers. You could confirm it returns a pseudoprime by computing

    pow(2, p-1, p)


on the return value to verify that it equals 1.

***

[1] Sometimes such a number is called a Fermat pseudoprime to distinguish it from other kinds of pseudo primes, i.e. numbers that satisfy other necessary conditions for being prime without actually being prime. For example, Perrin pseudoprimes. See the next post for more examples.