Thomas writes stuff

IBM Ponder This challenge - June 2019

Unlike May’s challenge that was solved by throwing raw power at it, this month’s Ponder This challenge could be solved elegantly and efficiently with a very old and well-known algorithm: Euclid’s algorithm.

First let’s remove some complexity by considering only powers of two instead of smooth numbers. We’re then looking for distincts numbers and such that . If we assume that all the are even and the are odd1 we can rewrite this as

The problem is now to find a linear combination of and that is small enough. Everyone knows how to use Euclid’s algorithm to compute the greatest common divisor of two numbers but what the algorithm really does is find linear combinations of the two numbers that are increasingly small. When we give it integers we eventually reach zero but if we give it irrational numbers it can continue forever and find linear combinations arbitrarily close to zero. We can therefore use it to find and .

For the challenge, we have which is below what double precision floating-point numbers can handle2. Also the naive implementation of Euclid’s algorithm is not numerically stable so we need extra margin3:

from decimal import *

# Set the precision to 100 decimals, should be enough.
getcontext().prec = 100

eps = Decimal('1e-15')

def euclid(a, b):
  if a - b < eps:
    return [1, -1]
    c = int(a / b)
    [x, y] = euclid(b, a - b * c)
    return [y, x - c * y]

[b, a] = euclid(Decimal(2).sqrt(), 1)

We now need to decompose the coefficients and as sums of powers of two to obtain the solution:

def powers_of_2(n):
  ans = []
  p = 1
  while n > 0:
    if (n / p) % 2 == 1:
      n = n - p
    p = 2 * p
  return ans

set_a = [x * x for x in powers_of_2(abs(a))]
set_b = [x * x * 2 for x in powers_of_2(abs(b))]
print(set_a, set_b)

# Check that the difference is indeed small enough.
sum_a = sum([Decimal(n).sqrt() for n in set_a])
sum_b = sum([Decimal(n).sqrt() for n in set_b])
print(sum_a - sum_b)

We have now solved the challenge. But because we only used powers of two, the solution that we obtain is quite large. Let’s try to get a more compact solution by allowing all smooth numbers.

Just like we did previously with the powers of two, we can show that it is sufficient to find a very small linear combination of . This is a know problem and there are several algorithms to solve it. I happened to know the Lenstra–Lenstra–Lovász algorithm which is one of them4. Let’s use the PARI/GP implementation:

? smooth = [];
? roots = [];
? {
  for(i = 0, 1,
    for(j = 0, 1,
      for(k = 0, 1,
	for(l = 0, 1,
	  s = 2^i*3^j*5^k*7^l;
	  smooth = concat(smooth, s);
	  roots = concat(roots, sqrt(s));
? coefficients = lindep(roots, 16)
%4 = [1, 2, 0, 6, 1, 7, -1, -6, 3, 2, 1, -4, 4, 4, -5, 0]~
? roots * coefficients
%5 = 5.6579457352476430067472915744709114061 E-16

All the integer coefficients that we found are already smooth numbers so there is no need to decompose them as sums of smooth numbers, we just need to square them and multiply them by their associated smooth number to get our solution:

? for(i = 1, 16, print(smooth[i] * coefficients[i] * abs(coefficients[i])))

We now have a solution with only 14 smooth numbers all below 4000.

  1. We can show that if we have a solution where this is not the case, we can rearrange it to have even powers of two on one side and odd ones on the other side. 

  2. The difference between 1.0 and the smallest double precision floating-point number greater than 1.0 is 2.22044604925031308e-16

  3. We could also put some effort in making Euclid’s algorithm numerically stable but we’re lazy and it’s easier to just use a very high precision. 

  4. Finding integer relations is not the original purpose of the LLL algorithm but it does it too. And I couldn’t easily find a usable implementation of the PSLQ algorithm which is supposed to be better for this. 

You can send your comments to stuff at thomash dot fr.