Teaser 2994: Consecutive sums
From The Sunday Times, 9th February 2020 [link] [link]
Amelia noticed that 15 is equal to 1+2+3+4+5 or 4+5+6 or 7+8, so there are three possible ways that it can be expressed as the sum of consecutive whole numbers. She then told Ben that she had found a three-digit number which can be expressed as the sum of consecutive whole numbers in just two different ways. “That’s interesting”, said Ben. “I’ve done the same, but my number is one more than yours”.
What is Ben’s number?
[teaser2994]
Frits, and 
S2T2 
Jim Randell 5:36 pm on 7 February 2020 Permalink |
The following Python program runs in 221ms.
Run: [ @replit ]
from enigma import (defaultdict, subsets, irange, tri, printf) # collect consecutive sums for 3-digit numbers d = defaultdict(list) for (i, j) in subsets(irange(1, 500), size=2): n = tri(j) - tri(i - 1) if 99 < n < 1000: d[n].append((i, j)) # find consecutive numbers expressible in 2 different ways for (n, vs) in d.items(): if len(vs) == 2: vs1 = d.get(n + 1, None) if vs1 and len(vs1) == 2: printf("{n} -> {vs}; {n1} -> {vs1}", n1=n + 1)Solution: Ben’s number is 289.
289 can be expressed as:
Amelia’s number is 288, which can be expressed as:
LikeLike
Jim Randell 6:36 pm on 7 February 2020 Permalink |
Here is a neater program that uses the technique described in the comments of Enigma 1. It runs in 83ms.
Run: [ @replit ]
from enigma import (divisors, irange, cache, printf) S = cache(lambda n: sum(d % 2 for d in divisors(n))) # find 2 consecutive 3-digit numbers with S(n) == 3 for a in irange(100, 998): b = a + 1 if S(a) == 3 and S(b) == 3: printf("{a}, {b}")Analytically:
In order for the numbers to have exactly 3 odd divisors they must be of the form:
where p is an odd prime.
And we are looking for two consecutive numbers, so one will be even and one will be odd (k = 0).
We can look for 3-digit numbers manually, or programatically:
from enigma import (primes, tuples, printf) # collect 3-digit numbers with S(n) = 3 ns = list() # consider p^2 term for p in primes.between(3, 31): n = p * p while n < 1000: if n > 99: ns.append(n) n *= 2 ns.sort() # find consecutive pairs for (a, b) in tuples(ns, 2): if b == a + 1: printf("{a}, {b}")We find the following 3-digit numbers with exactly 3 odd divisors:
And there are only two consecutive numbers in this list:
LikeLike
Frits 10:52 pm on 3 August 2020 Permalink |
Over the first 500 million numbers (> 99) only 2 consecutive numbers exist.
from enigma import Primes, printf, irange # collect 3-digit numbers with S(n) = 3 ns = list() test = 500000000 # consider p^2 term for p in Primes(test).irange(3, test): n = p * p #print(n) while n < test**2: if n > 99: ns.append(n) n *= 2 ns.sort() # find consecutive pairs for i in irange(0, len(ns) - 2): if ns[i+1] == ns[i] + 1: printf("{ns[i]}, {ns[i+1]}")Output:
288, 289
1681, 1682
LikeLike
Jim Randell 10:24 am on 4 August 2020 Permalink |
@Frits: I found some notes (and a program) that I wrote at the time on looking for larger solutions:
For consecutive numbers one of the following equations must hold:
where p and q are odd primes.
So we can consider solutions to the following forms of Pell’s Equation:
and look for solutions where x is an odd prime and y is an odd prime multiplied by a power of 2.
I have the following code (in a file called pells.py) that uses SymPy to solve Pell’s equations:
try: # newer versions of SymPy from sympy.solvers.diophantine.diophantine import diop_DN except ImportError: # older versions of SymPy from sympy.solvers.diophantine import diop_DN # interleave a sequence of iterators def interleave(ss): while ss: xs = list() for (i, s) in enumerate(ss): try: yield next(s) except StopIteration: xs.insert(0, i) for i in xs: del ss[i] # initial values from a sequence of iterators def initial(ss): for s in ss: try: yield next(s) except StopIteration: pass def solutions(D, N, x, y): # initial solution yield (D, N, x, y) # to go further we need the fundamental solution to (D, 1) if N == 1: (u, v) = (x, y) else: [s] = diop_DN(D, 1) (u, v) = map(int, s) while True: t = (x, y) (x, y) = (u * x + D * v * y, v * x + u * y) if (x, y) == t: break assert x * x - D * y * y == N yield (D, N, x, y) # generate solutions to Pell's equations # # x^2 - D.y^2 = N # # each equation may have multiple infinite families of solutions # # return values are (D, N, x, y), satisfying the above equation # # fn=interleave will generate solutions from each family interleaved # # fn=initial will generate just the first solution from each family # # Python 3 style definition: def pells(*DNs, fn=interleave): def pells(*DNs, **kw): fn = kw.get('fn', interleave) ss = list() for (D, N) in DNs: # consider fundamental solutions for s in diop_DN(D, N): (x, y) = map(int, s) ss.append(solutions(D, N, x, y)) return fn(ss)Which I then use in the following program:
Run: [ @replit ]
from pells import pells from enigma import (irange, is_prime_mr, drop_factors, first, printf) is_prime = lambda n: is_prime_mr(n, r=10) # look for possible solutions (from the first 60) for (D, N, x, y) in first(pells((2, 1), (2, -1)), 60): # extract p, q p = x if not (p > 2 and is_prime(p)): continue (k, q) = drop_factors(y, 2) if not (q > 2 and is_prime(q)): continue # output solution (a, b) = (x * x, 2 * y * y) printf("{a} - {b} = {N} [p={p} q={q}]")I used the program to check up to 6,000 solutions to the Pell’s equations (by which time the numbers under consideration are huge), and found the following pairs:
LikeLike
Frits 8:51 am on 7 August 2020 Permalink |
@Jim,
On PuzzlingInPython [{broken link removed}] the Pell equation has been reduced to:
(x_{k+1}, y_{k+1}) = (3x_k+4y_k,2x_k+3y_k)In an overnight run with this formula I checked up to 7000 solutions (only using one core) finding no more than you did.
LikeLike
Frits 8:10 pm on 3 August 2020 Permalink |
Where did you find formula n = (2^k)(p^2) ?
I only found a reference at OEIS A072502 ?
from enigma import irange, printf, tau # The requested numbers appear in OEIS as A072502 # Gary Detlefs suggested the tau equation S = lambda n: tau(2*n) == (tau(n) + 3) for a in irange(100, 998): b = a + 1 if S(a) and S(b) : printf("{a}, {b}")LikeLike
Jim Randell 9:58 am on 4 August 2020 Permalink |
@Frits:
If a number is a power of 2 (= 2^n) then it will have only 1 odd divisor: 1.
If a number has exactly 1 odd prime factor (= (2^n)p), then it will have exactly 2 odd divisors: 1, p.
If a number has exactly 2 different odd prime factors (= (2^n)(p.q)), then it will have exactly 4 odd divisors: 1, p, q, p.q.
But if the 2 odd prime factors are the same (so the number is (2^n)(p^2)), then it will have exactly 3 odd divisors: 1, p, p².
LikeLike
Jim Randell 7:39 am on 14 June 2025 Permalink |
I’ve written some code for solving the following general quadratic diophantine equation in 2 variables, without using SymPy:
The code is available from the enigma_plus repository as pells.py.
The program given above can then be revised to use the new solver:
from enigma import (is_prime_mr, drop_factors, merge, ifirst, printf) import pells is_prime = lambda n: is_prime_mr(n, r=10) # look for possible solutions (from the first 60 candidates) sols = (pells.diop_quad(1, -2, N) for N in [+1, -1]) for (x, y) in ifirst(merge(sols, uniq=1), 60): # extract p, q p = x if not (p > 2 and is_prime(p)): continue (k, q) = drop_factors(y, 2) if not (q > 2 and is_prime(q)): continue # output solution (a, b) = (x * x, 2 * y * y) printf("{a} - {b} = {N} [p={p} q={q}]", N=a - b)LikeLike