Teaser 3050: Power struggle
From The Sunday Times, 7th March 2021 [link] [link]
Given any number, one can calculate how close it is to a perfect square or how close it is to a power of 2. For example, the number 7 is twice as far from its nearest perfect square as it is from its nearest power of 2. On the other hand, the number 40 is twice as far from its nearest power of 2 as it is from its nearest square.
I have quite easily found a larger number (odd and less than a million!) for which one of these distances is twice the other.
What is my number?
[teaser3050]
Frits, and 
Hugh Casement are discussing. 
S2T2 
Jim Randell 4:45 pm on 5 March 2021 Permalink |
A brute force approach in Python can search the solution space in 215ms. If we stop after the first candidate solution is found, then it only takes 72ms.
Run: [ @replit ]
from enigma import (isqrt, irange, printf) # distance to nearest square def dist_sq(n): x = isqrt(n) return min(n - x**2, (x + 1)**2 - n) # distance to nearest power of 2 def dist_p2(n): x = n.bit_length() return min(2**x - n, n - 2**(x - 1)) # find solutions for n in irange(41, 1000000, step=2): d1 = dist_sq(n) d2 = dist_p2(n) if d1 == 2 * d2 or d2 == 2 * d1: printf("{n} -> dist_sq = {d1}, dist_p2 = {d2}") break # we only need the first solutionSolution: Your number is 32775.
There aren’t many odd numbers with this property.
The first 5 are:
LikeLike
Jim Randell 12:18 pm on 6 March 2021 Permalink |
Here’s an alternative approach that considers numbers at increasing distances from the sets of “markers”.
Also we can see that powers of 2 (greater than 1) are all even, and we are looking for an odd number n > 40, so its distance d to a power of 2 must be odd, and its distance to a square must be 2d. This reduces the number of checks we need to make.
The following Python program runs in 50ms.
Run: [ @replit ]
from enigma import (first, lt, powers, irange, inf, tuples, intersect, printf) # markers less than 1 million M = 1000000 count = lt(M) squares = first(powers(0, inf, 2), count=count) pow2s = first((2**i for i in irange(0, inf)), count=count) # generate numbers at a distance <d> from the nearest marker in <ms> def dist(ms, d): for (ml, m, mr) in tuples(ms, 3): n = m - d if n - ml > d: yield n n = m + d if mr - n > d: yield n # generate odd numbers at twice the distance from # a square than from a power of 2 def solve(): # consider increasing odd distances to a power of 2 for d in irange(1, inf, step=2): d2 = 2 * d # distance d from pow2s, 2d from squares for n in intersect([dist(pow2s, d), dist(squares, d2)]): printf("[{n} -> dist_p2 = {d}, dist_sq = {d2}]") yield n # find the first odd solution, greater than 40, less than M for n in solve(): if 40 < n < M: printf("answer = {n}") breakLikeLike
Jim Randell 9:30 am on 14 March 2021 Permalink |
This program finds the 5 odd solutions given above, and larger ones (by default it outputs the first 10 it finds), by examining the two squares surrounding each odd power of 2:
Run: [ @replit ]
from enigma import (irange, inf, isqrt, div, first, arg, printf) # generate candidate solutions def generate(): # consider increasing odd powers of 2 for i in irange(3, inf, step=2): p = 2**i # and the square below n = isqrt(p) # look for candidate solutions n2 = n * n n21 = n2 + 2 * n + 1 p2 = 2 * p # and check viability for x in [div(p2 + n2, 3), p2 - n2]: if x and x % 2 == 1 and x - n2 < n21 - x: printf("[{x} -> dist_p2 = {dp2}, dist_sq = {dsq}]", dp2=abs(p - x), dsq=x - n2) yield x for x in [div(p2 + n21, 3), p2 - n21]: if x and x % 2 == 1 and n21 - x < x - n2: printf("[{x} -> dist_p2 = {dp2}, dist_sq = {dsq}]", dp2=abs(p - x), dsq=n21 - x) yield x # find the first N solutions N = arg(10, 0, int) first(generate(), N)LikeLike
Hugh Casement 2:07 pm on 6 March 2021 Permalink |
I assume that “number” means integer: the setters of these puzzles have never been able to tell the difference.
The logarithm of n to base 2, q = ln(n) / ln(2). We round to the nearest integer and take 2^q.
Having written and run a Basic program in about one minute, I decided to look for even numbers.
There are dozens, starting with 54, 114, 278, 546, 982, 2002, 4182, 8370, …
The program was ten lines: I never expected Basic to be more compact than Python.
LikeLike
Frits 1:06 pm on 7 March 2021 Permalink |
A variation of Jim’s second program (more efficient).
from enigma import first, powers, irange, inf, tuples, intersect, printf, \ isqrt # markers less than 1 million M = 1000000 squares = first(powers(0, inf, 2), count=(lambda x: x < M)) pow2s = first((2 ** i for i in irange(0, inf)), count=(lambda x: x < M)) # generate numbers at a distance <d> from the nearest marker in <ms> def dist(ms, d): for (ml, m, mr) in tuples(ms, 3): n = m - d if n - ml > d: yield n n = m + d if mr - n > d: yield n # generate odd numbers at twice the distance from # a square than from a power of 2 def solve(): # consider increasing odd distances to a power of 2 for d in irange(1, inf, step=2): d2 = 2 * d for n in dist(pow2s, d): # is there a square at distance 2d from n? if n - d2 in squares or n + d2 in squares: # is it the nearest square? x = isqrt(n) x2 = x * x x2n = x2 + 2 * x + 1 nearest = x2 if n - x2 < x2n - n else x2n if nearest in {n - d2, n + d2}: printf("[{n} -> dist_p2 = {d}, dist_sq = {d2}]") yield n # find the first odd solution, greater than 40, less than M for n in solve(): if 40 < n < M: printf("answer = {n}") breakLikeLike
Jim Randell 2:53 pm on 7 March 2021 Permalink |
@Frits: I thought that a “hybrid” approach of my two programs would probably be faster. Although my second program has an internal runtime of 2.74ms, so I didn’t think it was worth the extra code.
Interestingly when I timed the code you posted it ran slightly slower than my version (3.03ms). It might be faster the other way around (generating distances from squares and then calculating distances from powers of two), as I would assume that [[
int.bit_length()]] is probably faster than [[math.isqrt()]] (certainly the code in enigma.py for [[isqrt()]] starts by calling [[int.bit_length()]], and then does more stuff. (Although if [[math.isqrt()]] is available it just uses that, which should be implemented by native code).LikeLike
Frits 5:01 pm on 7 March 2021 Permalink |
@Jim, I have tested with the Brian Gladman’s Timer function (which uses perf_counter_ns ) doing 200 runs (jra3 is my variation on your jra2 and jra4 a variation on jra2 the other way around (looking at distances from squares and then checking distances from powers of two)).
I would have been surprised in jra4 was faster than jra3 as the distances from squares generates many elements (like 1900).
LikeLike
Jim Randell 12:06 pm on 8 March 2021 Permalink |
@Frits: I suppose it goes to show that different timing strategies produce different results.
For a piece of Python code that is only going to be executed once I think that a 3ms internal runtime is perfectly fine. This run time is dwarfed by the amount of additional time the Python interpreter takes getting ready to run it, so in my preferred measure of execution time both the programs give an external runtime of 50ms.
For code that is going to be executed thousands of times it may be worth shaving a few milliseconds off it, but code that is only executed once I’m not so sure.
I wanted my programs to show two different approaches. My first program defines functions that calculate the two distances and then compares them. But there are nearly 500,000 candidate numbers, so this could involve a lot of calculation.
The second program makes a set of markers, and then looks at numbers that are fixed distances from the markers, so it is essentially a searching strategy. In this case the distances involved in the solution aren’t very large, so the searching strategy wins.
LikeLike
Frits 3:05 pm on 8 March 2021 Permalink |
@Jim,
Talking about isqrt:
have you considered using math.sqrt in enigma.py’s isqrt() if math.isqrt is not available (like in an up-to-date PyPy) for numbers under a certain limit? See Brian Gladman’s number_theory.py.
LikeLike
Jim Randell 6:46 pm on 8 March 2021 Permalink |
My original version of
isqrt()used a float square root to provide an initial guess, and then applied Newton’s method until the result settled down to an integer. (See this comment on Enigma 1759).But I switched that to a binary technique I found on Wikipedia [ @wikipedia ] which performed better for numbers that overflowed floats.
But looking at it again it turns out that that was because I was choosing a poor initial guess in the overflow case. It turns out you can do a much better guess using the fact that
int.bit_length()behaves roughly likelog2().So I think I’ll switch back to using Newton’s method for cases where
math.isqrt()is not available.BTW: For a very fast implementation of
isqrt()check outgmpy2.isqrt().LikeLike
Hugh Casement 1:20 pm on 7 March 2021 Permalink |
Isn’t this simpler?
defdbl a - z for n = 41.0 to 999999.0 step 2.0 p = int(log(n)/log(2.0) + 0.5) pp = 2.0^p dp = abs(pp - n) s = int(sqr(n) + 0.5) ss = s*s ds = abs(ss - n) if dp = 2.0*ds or ds = 2.0*dp then print n pp ss next nTurbo Basic has a log2 function, which would have made the third line a bit more compact.
Of course some expressions could be combined, but at a slight loss of transparency.
The .0 are not strictly necessary but make sure the numbers are interpreted as floating-point (single-length integers go only to 32767).
LikeLike
Jim Randell 1:55 pm on 7 March 2021 Permalink |
@Hugh: (I put your code into
[code] ... [/code]tags so it should look better now).It is simpler, but it doesn’t calculate accurate values:
e.g. The nearest power of 2 to 46 is 32, which gives a distance of 14:
>>> dist_p2(46)
14
>>> f = lambda n: abs(2 ** int(log2(n) + 0.5) – n)
>>> f(46)
18
You have to go much higher before the calculation of the nearest square breaks down (much more than a million), probably due to exceeding the precision of floating point numbers. (e.g. It fails at 10^37).
LikeLike
Frits 2:00 pm on 7 March 2021 Permalink |
Hi Hugh,
Your program suffers from the same problem Erling Torkildsen had at PuzzlingInPython (dp should be 107 iso 149 for n = 363).
from math import log2 for n in range(41, 1000000, 2): #p = int(log2(n) + 0.5) p = round(log2(n)) pp = 2 ** p dp = abs(pp - n) #s = int(n ** 0.5 + 0.5) s = round(n ** 0.5) ss = s * s ds = abs(ss - n) if n == 363: print("n =", n, "p =", p, "pp =", pp, "dp =", dp, "ds =", ds, "s =", s, "ss =", ss) if dp == 2 * ds or ds == 2 * dp: print(n, pp, ss) break # n = 363 p = 9 pp = 512 dp = 149 ds = 2 s = 19 ss = 361LikeLike
Frits 4:52 pm on 9 March 2021 Permalink |
Looping over the power of 2 entries.
# suppose p2 (2 ** p) is the power of 2 is the nearest power # let s^2 be the nearest square smaller equal to p2 # so that s^2 <= p2 <= (s + 1)^2 # example: 484 <= 512 <= 529 with s = 22 and p2 = 2 ** 9 # +--- 2s-1 ---+----- 2s + 1 ----+--- 2s + 3 --+ # | | | | # ---- (s - 1)^2 ---- s^2 --- p2 ---- (s + 1)^2 --- (s + 2)^2 --- # # even power exponents are also squares so ignore them for p in range(5, 21, 2): p2 = 2 ** p s = int(sqrt(p2)) s2 = s ** 2 # can x be between (s - 1)^2 and s^2 ? # then ds is maximal s - 1 and d2 is maximal (s - 1) / 2 minx = min(p2 - ((s - 1) // 2), s2) # make sure minx is odd and ... minx = minx + (minx % 2 == 0) # ... ds is not a multiple of 4 and minx = minx + 2 * (minx % 4 == 1) # if s2 is even then first entries have an odd ds if s2 % 2 == 0: minx += s - 2 # can x be between (s + 1)^2 and (s + 2)^2 ? # then ds is maximal s + 1 and d2 is maximal (s + 1) / 2 maxx = max(p2 + ((s + 1) // 2), s2 + 2 * s + 1) # if s2 is odd then last entries have an odd ds if s2 % 2 : maxx -= s + 1 # consider odd x's so that ds is not a multiple of 4 for x in range(minx, maxx + 1, 4): # determine nearest square y = int(sqrt(x)) y2 = y ** 2 ds = min(x - y2, y2 + 2 * y + 1 - x) # determine nearest power of 2 y2min1 = 2 ** (x.bit_length() - 1) d2 = min(2 * y2min1 - x, x - y2min1) # distance from x to power of 2 is always odd # so ds = 2 * d2 (and not d2 = 2 * ds) if ds == 2 * d2: print(f"number is {x}; distance to power of 2: {d2}; to square: {ds}") exit() # we only need the firstLikeLike