Archive for June 25th, 2009

Problem 14 – 3n+1

June 25, 2009

The “3n+1 problem” is one I mention to my calculus class every time we get to sequences. It’s just so fun.

I’ve coded up a few solutions. My first is a basic one where I have a function that takes a number and generates the sequence starting at that number. Then I just loop up to 1000000 looking at the lengths of these sequences.

```def seq(num):
""" Returns the 3n+1 sequence for num """
ret = [num]
while(not ret[-1] == 1):
ret.append(ret[-1]%2 == 0 and ret[-1]//2 or 3*ret[-1]+1)
return ret

def solve(bound):
""" Find the number less than bound with largest 3n+1 sequence

Returns the number, the length of its sequence, and the computation time
"""
t = time.time()
retn,retl = 0,0
for n in xrange(1,bound):
s = seq(n)
if(len(s) > retl):
retn,retl = n,len(s)
t = time.time() - t
return ((retn,retl), t)

if __name__ == "__main__":
print solve(int(sys.argv[1]))
```

This can be moderately improved by not passing around the entire sequence. Instead, just code up the helper function (“seq”, above) to return the length, and use it accordingly in the main loop (“solve”, above). I’ll not bother copying that code here.

The only real improvement I’ve made on this is to store values as we find them, so we can look them up again later. Here’s the way I coded that up:

```def lenseq(n, mem):
""" Recursively calculate the length of n using a memory dictionary

This function will modify mem, adding new values to it
"""
if(n in mem):
return mem[n]
suc = n%2 and 3*n+1 or n//2
lenseq(suc, mem)
mem[n] = mem[suc] + 1
return mem[n]

def solve(bound):
""" Find the number less than bound with longest 3n+1 sequence

Return ((the number, the length of its sequence), computation time)
"""
t = time.time()
retn,retl = 0,0
mem = {1:1}
for n in xrange(1,bound):
l = lenseq(n, mem)
if(l > retl):
retn,retl = n,l
t = time.time() - t
return ((retn,retl),t)

if __name__ == "__main__":
print solve(int(sys.argv[1]))
```

And I did a little performance comparison, because it was easy. In the chart below, “Basic” is the first set of code above, “BasicPassLength” is the modification that doesn’t pass sequences around, just their lengths, and “StoreLengths” is the second set of code above.

The scale is log in both axes. The tickmark, $n$, on the $x$-axis refers to an upper bound of $10^n$ (the last tickmark, then, being the stated problem). The $y$-axis is less precise, but it’s along the lines of the log of the time taken, in seconds. The “Basic” code took about 2.5 minutes on the stated problem (bound of 1000000), “BasicPassLength” about 1.5 minutes, and the “StoreLengths” took about 20 seconds.

Just out of curiosity, I checked how many values ended up getting stored. There were 2168611 values stored, ranging up to 56991483520 (which, in case you are curious, has length 174). I haven’t thought up an obvious way to trim down which values get stored (maybe don’t do those above the bound?), or other ways to improve on this solution.