Problem 167 on Project Euler works with Ulam numbers, also referred to as 1-additive sequences in several mathematical papers. Reproduction of the problem statement:
For two positive integers a and b, the Ulam sequence \(U(a,b)\) is defined by \(U(a,b)_1 = a\), \(U(a,b)_2 = b\) and for \(k > 2\), \(U(a,b)_k\) is the smallest integer greater than \(U(a,b)_{k-1}\) which can be written in exactly one way as the sum of two distinct previous members of \(U(a,b)\).
For example, the sequence \(U(1,2)\) begins with
$$1, 2, 3 = 1 + 2, 4 = 1 + 3, 6 = 2 + 4, 8 = 2 + 6, 11 = 3 + 8;$$5 does not belong to it because \(5 = 1 + 4 = 2 + 3\) has two representations as the sum of two previous members, likewise \(7 = 1 + 6 = 3 + 4\).
Find \(\sum_{n=2}^10 U(2,2n+1)_k\), where \(k = 10^11\).
The first thing to note is that all of the sequences in the summation have precisely two even numbers: \(2\) and \(4n + 4\). This is proven here, so the property can be relied upon to hold true well into larger values.
Additionally, it is proven here that any 1-additive sequence with finitely even terms is eventually regular. A sequence is regular if successive differences repeat.
These properties lead to a two-phase solution. Compute enough values to determine a repeating period, using the period size and cumulative difference to compute the value for \(k=10^11\).
Because there are only two even values in the sequence, one of these two values must be a component in the sum of every odd number. Whenever a new odd value is computed for a sequence, it introduces two candidates into the pool of potential values, which are stored into a priority queue. The next value in the sequence can then be determined by pulling values out of the queue until a non-repeating entry is found.
from Queue import PriorityQueue
class U(object):
def __init__(self, v):
self.computed = []
self.candidates = PriorityQueue()
self.candidates.put(2)
self.candidates.put(v)
self.candidates.put(2*v + 2)
def compute(self, k):
while len(self.computed) < k:
candidate = self.candidates.get()
valid = True
while not self.candidates.empty() and self.candidates.queue[0] == candidate:
valid = False
self.candidates.get()
if valid:
self.computed.append(candidate)
if candidate % 2 == 1:
self.candidates.put(2 + candidate)
self.candidates.put(2*self.computed[1] + 2 + candidate)
return self.computed[k - 1]With a fast way of directly computing values, enough can be generated to determine a period. This is accomplished with a lazily-evaluated sequence of differences, checking increasingly large period sizes until one is found that repeats.
class Period(object):
def __init__(self, u, period_start):
self.period_start = period_start
self.minimum_period = 20
self.diffs = []
self.previous = u.compute(period_start)
self.u = u
self.period = self.compute_period()
def compute_period(self):
period_length = self.minimum_period
while not self.verify_period(period_length):
period_length += 1
return self.diffs[0:period_length]
def verify_period(self, period_length):
i = period_length
while i <= 2*period_length:
if self.lookup_diff(i) != self.lookup_diff(i % period_length):
return False
i += 1
return True
def lookup_diff(self, i):
while i >= len(self.diffs):
s = self.u.compute(self.period_start + 1 + len(self.diffs))
self.diffs.append(s - self.previous)
self.previous = s
return self.diffs[i]
def compute(self, k):
if k <= self.period_start:
return self.u.compute(k)
fundamental_difference = sum(self.period)
base = self.u.compute(self.period_start) + fundamental_difference*((k - self.period_start)/len(self.period))
return base + sum(self.period[0:(k - self.period_start)%len(self.period)])Now it's possible to compute the final answer, which takes about two minutes, most of which is spent computing the period in the last couple iterations. The intermediate values are printed so they can be copied to unit tests used during algorithm refinement.
result = 0
for n in range(2, 11):
v = 2*n + 1
u = U(v)
period = Period(u, 100)
x = period.compute(10**11)
print str(v) + ": " + str(x)
print len(period.period)
result += x
print result