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 wellknown 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 odd^{1} 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 floatingpoint numbers can handle^{2}. Also the naive implementation of Euclid’s algorithm is not numerically stable so we need extra margin^{3}:
We now need to decompose the coefficients and as sums of powers of two to obtain the solution:
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 them^{4}. 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 E16
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.

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. ↩

The difference between
1.0
and the smallest double precision floatingpoint number greater than1.0
is2.22044604925031308e16
. ↩ 
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. ↩

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. ↩