Hofstadter-Conway $10,000 sequence
You are encouraged to solve this task according to the task description, using any language you may know.
The definition of the sequence is colloquially described as:
- Starting with the list [1,1],
- Take the last number in the list so far: 1, I'll call it x.
- Count forward x places from the beginning of the list to find the first number to add (1)
- Count backward x places from the end of the list to find the second number to add (1)
- Add the two indexed numbers from the list and the result becomes the next number in the list (1+1)
- This would then produce [1,1,2] where 2 is the third element of the sequence.
Note that indexing for the description above starts from alternately the left and right ends of the list and starts from an index of one.
A less wordy description of the sequence is:
a(1)=a(2)=1 a(n)=a(a(n-1))+a(n-a(n-1))
The sequence begins:
1, 1, 2, 2, 3, 4, 4, 4, 5, ...
Interesting features of the sequence are that:
- a(n)/n tends to 0.5 as n grows towards infinity.
- a(n)/n where n is a power of 2 is 0.5
- For n>4 the maximal value of a(n)/n between successive powers of 2 decreases.
The sequence is so named because John Conway offered a prize of $10,000 to the first person who could find the first position, p in the sequence where
|a(n)/n| < 0.55 for all n >= p.
It was later found that Hofstadter had also done prior work on the sequence.
The 'prize' was won quite quickly by Dr. Colin L. Mallows who proved the properties of the sequence and allowed him to find the value of n. (Which is much smaller than the 3,173,375,556. quoted in the NYT article)
The task is to:
- Create a routine to generate members of the Hofstadter-Conway $10,000 sequence.
- Use it to show the maxima between successive powers of two up to 2**20
- As a stretch goal: Compute the value of n that would have won the prize and confirm it is true for n up to 2**20
References:
- Conways Challenge Sequence, Mallows' own account.
- Mathworld Artcle.
C
<lang C>#include <stdio.h>
- include <stdlib.h>
int a_list[1<<20 + 1];
int doSqnc( int m) {
int max_df = 0; int p2_max = 2; int v, n; int k1 = 2; int lg2 = 1; double amax = 0; a_list[0] = -50000; a_list[1] = a_list[2] = 1; v = a_list[2];
for (n=3; n <= m; n++) {
v = a_list[n] = a_list[v] + a_list[n-v];
if ( amax < v*1.0/n) amax = v*1.0/n;
if ( 0 == (k1&n)) {
printf("Maximum between 2^%d and 2^%d was %f\n", lg2,lg2+1, amax);
amax = 0;
lg2++;
}
k1 = n;
}
return 1;
}</lang> Results
Maximum between 2^1 and 2^2 was 0.666667 Maximum between 2^2 and 2^3 was 0.666667 Maximum between 2^3 and 2^4 was 0.636364 Maximum between 2^4 and 2^5 was 0.608696 .... Maximum between 2^18 and 2^19 was 0.534645 Maximum between 2^19 and 2^20 was 0.533779
J
Solution (tacit): <lang j> hc10k =: , ] +/@:{~ (,&<: -.)@{: NB. Actual sequence a(n)
AnN =: % 1+i.@:# NB. a(n)/n MxAnN =: >./;.1~ 2 (=<.)@:^. 1+i.@# NB. Maxima of a(n)/n between successive powers of 2</lang>
Alternative solution (exponential growth):
The first, naive, formulation of hc10k grows by a single term every iteration; in this one, the series grows exponentially in the iterations.
<lang j> hc10kE =: 1 1 , expand@tail
expand =: 2+I.@;
tail =: copies&.>^:(<@>:`(<@,@2:))
copies =: >: |.@(#!.1 |.)~ 1 j. #;.1 #^:_1 ::1:~ ]~:{.</lang>
Example: <lang j> ] A=:1 1 hc10k @]^:[~ 2^20x 1 1 2 2 3 4 4 4 5 6 7 7 8 8 8 8 9 ...
AnN A
1 0.5 0.666667 0.5 0.6 0.666667 ...
MxAnN@AnN A
1 0.666667 0.666667 0.636364 ...
MxAnN@AnN@hc10kE 20
1 0.666667 0.666667 0.636364 ...</lang>
Oz
A direct implementation of the recursive definition with explicit memoization using a mutable map (dictionary): <lang oz>declare
local
Cache = {Dictionary.new}
Cache.1 := 1
Cache.2 := 1
in
fun {A N}
if {Not {Dictionary.member Cache N}} then
Cache.N := {A {A N-1}} + {A N-{A N-1}}
end
Cache.N
end
end
Float = Int.toFloat
for I in 0..19 do
Range = {List.number {Pow 2 I} {Pow 2 I+1} 1}
RelativeValues = {Map Range
fun {$ N}
{Float {A N}}
/ {Float N}
end}
Maximum = {FoldL RelativeValues Max 0.0}
in
{System.showInfo "Max. between 2^"#I#" and 2^"#I+1#": "#Maximum}
end</lang>
Output:
Max. between 2^0 and 2^1: 1.0 Max. between 2^1 and 2^2: 0.66667 Max. between 2^2 and 2^3: 0.66667 Max. between 2^3 and 2^4: 0.63636 Max. between 2^4 and 2^5: 0.6087 Max. between 2^5 and 2^6: 0.59091 Max. between 2^6 and 2^7: 0.57609 Max. between 2^7 and 2^8: 0.56742 Max. between 2^8 and 2^9: 0.55946 Max. between 2^9 and 2^10: 0.55494 Max. between 2^10 and 2^11: 0.5501 Max. between 2^11 and 2^12: 0.54746 Max. between 2^12 and 2^13: 0.54414 Max. between 2^13 and 2^14: 0.54244 Max. between 2^14 and 2^15: 0.54007 Max. between 2^15 and 2^16: 0.53878 Max. between 2^16 and 2^17: 0.53704 Max. between 2^17 and 2^18: 0.53602 Max. between 2^18 and 2^19: 0.53465 Max. between 2^19 and 2^20: 0.53378
PicoLisp
<lang PicoLisp>(de hofcon (N)
(cache '*HofCon (format (seed N))
(if (>= 2 N)
1
(+
(hofcon (hofcon (dec N)))
(hofcon (- N (hofcon (dec N)))) ) ) ) )
(setq *Scl 20)
(de sequence (M)
(let (Lim 4 Max 0 4k$ 0)
(for (N 3 (>= M N) (inc N))
(let V (*/ (hofcon N) 1.0 N)
(setq Max (max Max V))
(when (>= V 0.55)
(setq 4k$ N) )
(when (= N Lim)
(prinl
"Maximum between " (/ Lim 2)
" and " Lim
" was " (format Max `*Scl) )
(inc 'Lim Lim)
(zero Max) ) ) )
(prinl
"Win with " (inc 4k$)
" (the task requests 'n >= p')" ) ) )
(sequence (** 2 20))</lang> Output:
Maximum between 2 and 4 was 0.66666666666666666667 Maximum between 4 and 8 was 0.66666666666666666667 Maximum between 8 and 16 was 0.63636363636363636364 Maximum between 16 and 32 was 0.60869565217391304348 Maximum between 32 and 64 was 0.59090909090909090909 Maximum between 64 and 128 was 0.57608695652173913043 Maximum between 128 and 256 was 0.56741573033707865169 Maximum between 256 and 512 was 0.55945945945945945946 Maximum between 512 and 1024 was 0.55493741307371349096 Maximum between 1024 and 2048 was 0.55010087424344317418 Maximum between 2048 and 4096 was 0.54746289264756644805 Maximum between 4096 and 8192 was 0.54414474786396381303 Maximum between 8192 and 16384 was 0.54244270878036220067 Maximum between 16384 and 32768 was 0.54007109751158709445 Maximum between 32768 and 65536 was 0.53878402058425570614 Maximum between 65536 and 131072 was 0.53704365699986594575 Maximum between 131072 and 262144 was 0.53602006781156104419 Maximum between 262144 and 524288 was 0.53464543107811232092 Maximum between 524288 and 1048576 was 0.53377922996336783427 Win with 1490 (the task requests 'n >= p')
Python
Note how we allow for Pythons zero based indexing: <lang python>from __future__ import print_function, division
def hc10000(lst):
Modifies lst in-place, adding successive terms of the sequence
lst should start as [None, 1, 1] to allow for zero-based indexing
Doctest:
>>> lst = [None,1,1]
>>> hc10000(lst); lst
[None, 1, 1, 2]
>>> hc10000(lst); lst
[None, 1, 1, 2, 2]
>>> hc10000(lst); lst
[None, 1, 1, 2, 2, 3]
>>> hc10000(lst); lst
[None, 1, 1, 2, 2, 3, 4]
>>>
last = lst[-1]
lst.append(lst[last] + lst[-last])</lang>
Finding the maxima and the answer to the $10,000 question: <lang python>def maxandmallows(nmax=2**20):
power2, maxima, lst = 4, [0.5], [None,1,1]
for n in range(2,nmax):
last = lst[-1]
m = last/n
mallows = abs(m-0.5) -0.05
if mallows>0:
mm = n, mallows
if maxima[-1] < m:
maxima[-1] = m
if n==power2:
print("In the region %7i < n <= %7i: max a(n)/n = %6.4f" % (
power2/2, power2, maxima[-1]))
power2 *= 2
maxima.append(0.5)
hc10000(lst) # Quicker: lst.append(lst[last] + lst[-last])
return mm[0], maxima
mallowsmax, maxima = maxandmallows(2**20) print("\nYou too might have won $1000 with an answer of n =", mallowsmax)</lang>
Sample output:
In the region 2 < n <= 4: max a(n)/n = 0.6667 In the region 4 < n <= 8: max a(n)/n = 0.6667 In the region 8 < n <= 16: max a(n)/n = 0.6364 In the region 16 < n <= 32: max a(n)/n = 0.6087 In the region 32 < n <= 64: max a(n)/n = 0.5909 In the region 64 < n <= 128: max a(n)/n = 0.5761 In the region 128 < n <= 256: max a(n)/n = 0.5674 In the region 256 < n <= 512: max a(n)/n = 0.5595 In the region 512 < n <= 1024: max a(n)/n = 0.5549 In the region 1024 < n <= 2048: max a(n)/n = 0.5501 In the region 2048 < n <= 4096: max a(n)/n = 0.5475 In the region 4096 < n <= 8192: max a(n)/n = 0.5441 In the region 8192 < n <= 16384: max a(n)/n = 0.5424 In the region 16384 < n <= 32768: max a(n)/n = 0.5401 In the region 32768 < n <= 65536: max a(n)/n = 0.5388 In the region 65536 < n <= 131072: max a(n)/n = 0.5370 In the region 131072 < n <= 262144: max a(n)/n = 0.5360 In the region 262144 < n <= 524288: max a(n)/n = 0.5346 You too might have won $1000 with an answer of n = 1489
Tcl
The routine to return the nth member of the sequence. <lang tcl>package require Tcl 8.5
set hofcon10k {1 1} proc hofcon10k n {
global hofcon10k
if {$n < 1} {error "n must be at least 1"}
if {$n <= [llength $hofcon10k]} {
return [lindex $hofcon10k [expr {$n-1}]]
}
while {$n > [llength $hofcon10k]} {
set i [lindex $hofcon10k end] set a [lindex $hofcon10k [expr {$i-1}]] # Don't use end-based indexing here; faster to compute manually set b [lindex $hofcon10k [expr {[llength $hofcon10k]-$i}]] lappend hofcon10k [set c [expr {$a + $b}]]
} return $c
}</lang> The code to explore the sequence, looking for maxima in the ratio. <lang tcl>for {set p 1} {$p<20} {incr p} {
set end [expr {2**($p+1)}]
set maxI 0; set maxV 0
for {set i [expr {2**$p}]} {$i<=$end} {incr i} {
set v [expr {[hofcon10k $i] / double($i)}] if {$v > $maxV} {set maxV $v; set maxI $i}
}
puts "max in 2**$p..2**[expr {$p+1}] at $maxI : $maxV"
}</lang> Output:
max in 2**1..2**2 at 3 : 0.6666666666666666 max in 2**2..2**3 at 6 : 0.6666666666666666 max in 2**3..2**4 at 11 : 0.6363636363636364 max in 2**4..2**5 at 23 : 0.6086956521739131 max in 2**5..2**6 at 44 : 0.5909090909090909 max in 2**6..2**7 at 92 : 0.5760869565217391 max in 2**7..2**8 at 178 : 0.5674157303370787 max in 2**8..2**9 at 370 : 0.5594594594594594 max in 2**9..2**10 at 719 : 0.5549374130737135 max in 2**10..2**11 at 1487 : 0.5501008742434432 max in 2**11..2**12 at 2897 : 0.5474628926475664 max in 2**12..2**13 at 5969 : 0.5441447478639638 max in 2**13..2**14 at 11651 : 0.5424427087803622 max in 2**14..2**15 at 22223 : 0.5400710975115871 max in 2**15..2**16 at 45083 : 0.5387840205842557 max in 2**16..2**17 at 89516 : 0.5370436569998659 max in 2**17..2**18 at 181385 : 0.5360200678115611 max in 2**18..2**19 at 353683 : 0.5346454310781124 max in 2**19..2**20 at 722589 : 0.5337792299633678