Teaser 3200: Puzzling sum
From The Sunday Times, 21st January 2024 [link] [link]
Liam was given a set of addition sums as part of his homework. The teacher had cleverly chosen sums where the answer, together with the two numbers to be added, used all the digits from 1 to 9 just once (and there were no zeros). For instance, one of the sums was 193 + 275 = 468. I noticed that one of them was particularly interesting in that the two numbers to be added were perfect powers (squares, cubes etc).
In the interesting sum, what are the two numbers that were to be added?
[teaser3200]
Brian Gladman, 
GeoffR, and 1 other are discussing. 
S2T2 
Jim Randell 4:48 pm on 19 January 2024 Permalink |
This Python program uses the [[
SubstitutedExpression()]] solver from the enigma.py library to consider all possible sums that the teacher could have set, and then looks for those with summands that are perfect powers.It runs in 87ms.
Run: [ @replit ]
from enigma import ( SubstitutedExpression, decompose, irange, prime_factor, mgcd, unpack, sprintf as f ) # symbols for the 9 different digits syms = "ABCDEFGHI" # check a number is a perfect power (> 1) def is_ipower(n, gcd=unpack(mgcd)): if n == 1: return True return gcd(e for (p, e) in prime_factor(n)) > 1 # split the symbols into 2 addends and a result (x + y = z) for (a, b, c) in decompose(len(syms), 3, increasing=1, sep=0, min_v=1, max_v=4): # create the sum (x, y, z) = (syms[:a], syms[a:a + b], syms[a + b:]) # create an alphametic puzzle exprs = [ f("{x} + {y} = {z}") ] if len(x) == len(y): exprs.append(f("{x0} < {y0}", x0=x[0], y0=y[0])) p = SubstitutedExpression( exprs, digits=irange(1, 9), answer=f("({x}, {y}, {z})"), verbose='', ) # and solve it for (x, y, z) in p.answers(): if is_ipower(x) and is_ipower(y): # output solution print(f("{x} + {y} = {z}"))Solution: The interesting sum is: 243 + 576 = 819.
And:
There are 168 possible sums using the digits 1-9, but only this one has summands that are perfect powers. And all of them are of the form:
ABC + DEF = GHI.LikeLike
Jim Randell 8:10 pm on 19 January 2024 Permalink |
An alternative approach is just to look for two powers that give a sum with the required properties. (The code to generate powers is borrowed from Teaser 3119).
This Python program has an internal runtime of 1.2ms.
Run: [ @replit ]
from enigma import (Primes, nsplit, printf) # generate powers (x^y) in numerical order, x >= 0, y >= 2, def ipowers(): # powers less than 2^2 yield 0 yield 1 # powers from 2^2 upwards base = { 2: 2 } power = { 2: 4 } maxp = 2 primes = Primes(35, expandable=1).generate(maxp + 1) while True: # find the next power n = min(power.values()) yield n # what powers are involved ps = list(p for (p, v) in power.items() if v == n) # increase the powers for p in ps: base[p] += 1 power[p] = pow(base[p], p) # do we need to add in a new power? if maxp in ps: maxp = next(primes) base[maxp] = 2 power[maxp] = pow(2, maxp) # check for valid digits, must be all different and not include 0 valid = lambda ds: 0 not in ds and len(set(ds)) == len(ds) # collect powers ps = list() for p in ipowers(): ds = nsplit(p) if len(ds) > 4: break if not valid(ds): continue # look for a smaller power to pair with this one for (p1, ds1) in ps: t = p1 + p dst = nsplit(t) # check all 9 digits used ds9 = ds + ds1 + dst n = len(ds9) if n < 9: continue if n > 9: break if not valid(ds9): continue # output solution printf("{p1} + {p} = {t}") # record this power ps.append((p, ds))I will add the [[
ipowers()]] function to the enigma.py library (from version2024-01-19).LikeLike
Frits 1:53 am on 20 January 2024 Permalink |
@Jim, did you forget to check dst for zeroes?
LikeLike
Jim Randell 7:31 am on 20 January 2024 Permalink |
Thanks. Fixed now.
LikeLike
Brian Gladman 6:53 pm on 20 January 2024 Permalink |
@Jim. That is a very neat subroutine for generating integer powers in order.
LikeLike
Brian Gladman 9:21 pm on 21 January 2024 Permalink |
@Jim, I looked at alternatives to your ipowers() routine and found this quite a bit faster:
from heapq import heappush, heappop def ipowers(maths=True): if maths: yield 0 yield 1 heap, pv, m = [], None, 2 heappush(heap, (2 ** 2, 2, 2)) while True: m += 1 v = 2 ** m while heap[0][0] <= v: (v1, n1, m1) = heappop(heap) if v1 != v and v1 != pv: yield v1 pv = v1 heappush(heap, ((n1 + 1) ** m1, n1 + 1, m1)) heappush(heap, (2 ** m, 2, m))LikeLike
Jim Randell 10:37 am on 22 January 2024 Permalink |
@Brian: Yes, using a heap is faster (especially in CPython which I think has heaps implemented via a C module).
So, this is even faster:
from heapq import (heapify, heappush, heappop) def ipowers(): # powers less than 4 yield 0 yield 1 # powers from 4 upwards hi = 1 maxe = 2 pows = [(4, 2, 2)] heapify(pows) exps = Primes(35, expandable=1).generate(3) while True: # find the next power (p, b, e) = heappop(pows) if p > hi: yield p hi = p heappush(pows, (pow(b + 1, e), b + 1, e)) # do we need to add in a new exponent? if b == 2: maxe = next(exps) heappush(pows, (pow(2, maxe), 2, maxe))LikeLiked by 1 person
Brian Gladman 9:02 am on 23 January 2024 Permalink |
It is faster if you don’t count the cost of the extra import of Primes() which dwarfs the running cost for low power limits. Including import cost for powers below 1000 I measured 338 microseconds and 2.5 milliseconds on CPython. In my test, it overtakes the non-import version at a power limit of about 10^9.
LikeLike
Jim Randell 1:47 pm on 23 January 2024 Permalink |
I was more interested in algorithmic efficiency, so I was testing by generating the first 5 million perfect powers, and using primes for exponents was a clear winner on all the Python environments I tested.
I generally have the enigma module loaded, so there is essentially no overhead in using the primes generator, but the candidate exponents can be any superset of the primes (as the algorithm will discard duplicates), so you can use something like [[
irange(3, inf, step=2)]] if you don’t want to use primes, and still get a performance improvement.LikeLike
GeoffR 6:55 pm on 19 January 2024 Permalink |
from enigma import nsplit # A manual search for all three digit powers gave the following list d3_powers = [125, 128, 216, 243, 256, 512, 529, 576, 625, 729, 784, 961] for r in range(123, 987): # result of the sum g, h, i = nsplit(r) if len(set((g, h, i))) != 3:continue if r in d3_powers:continue # find the two summands, which must both be powers for n in d3_powers: for m in d3_powers: if n == m:continue a, b, c = nsplit(n) d, e, f = nsplit(m) if len(set((a, b, c, d, e, f))) == 6: if m + n == r: if m > n and len(set((a,b,c,d,e,f,g,h,i))) == 9: print(f"The sum is {m} + {n} = {r}.")LikeLike
GeoffR 8:27 pm on 19 January 2024 Permalink |
Seems faster to use the remaining digits, after the two squares, for the addition result.
from enigma import nsplit # A manual search for all three digit powers gave the following list d3_powers = [125, 128, 216, 243, 256, 512, 529, 576, 625, 729, 784, 961] # find the two summands, which must both be powers for n in d3_powers: for m in d3_powers: if n == m:continue a, b, c = nsplit(n) d, e, f = nsplit(m) if len(set((a, b, c, d, e, f))) == 6: # find result of additon from remaining digits g, h, i = set(range(1,10)).difference((set((a, b, c, d, e, f)))) if len(set((g, h, i))) != 3:continue r = 100*g + 10*h + i if r in d3_powers:continue if m + n == r: if m > n and len(set((a, b, c, d, e, f, g, h, i))) == 9: print(f"The sum is {m} + {n} = {r}.")LikeLike
GeoffR 2:50 pm on 22 January 2024 Permalink |
% A Solution in MiniZinc include "globals.mzn"; var 1..9:A; var 1..9:B; var 1..9:C; var 1..9:D; var 1..9:E; var 1..9:F; var 1..9:G; var 1..9:H; var 1..9:I; var 100..999:ABC == 100*A + 10*B + C; var 100..999:DEF == 100*D + 10*E + F; var 100..999:GHI == 100*G + 10*H + I; constraint all_different ([A,B,C,D,E,F,G,H,I]); constraint ABC + DEF == GHI; set of int: powers == {125,128,216,243,256,512,529,576,625,729,784,961}; constraint ABC > DEF /\ ABC in powers /\ DEF in powers; solve satisfy; output ["Sum is " ++ show(ABC) ++ " + " ++ show(DEF) ++ " = " ++ show(GHI)];LikeLike