Hofstadter-Conway $10,000 sequence

From Rosetta Code
Revision as of 16:24, 8 December 2010 by 87.194.249.72 (talk) (Added C# version.)
Task
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.
a(n) / n for n in 1..256
a(n) / n for n in 1..256


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:

  1. Create a routine to generate members of the Hofstadter-Conway $10,000 sequence.
  2. Use it to show the maxima of a(n)/n between successive powers of two up to 2**20
  3. 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:

Ada

<lang Ada> with Ada.Text_IO; use Ada.Text_IO;

procedure Test_Hofstadter_Conway is

  type Sequence is array (Positive range <>) of Positive;
  function Hofstadter_Conway (N : Positive) return Sequence is
  begin
     return A : Sequence (1..N) do
        A (1) := 1;
        if N > 1 then
           A (2) := 1;
           for K in 3..N loop
              A (K) := A ( A (K - 1)) + A (K - A (K - 1));
           end loop;
        end if;
     end return;
  end Hofstadter_Conway;
  
  S : Sequence := Hofstadter_Conway (2**20);
  Mallows : Positive;

begin

  for N in 1..19 loop
     declare
        Max   : Float := 0.0;
        This  : Float;
        Where : Positive;
     begin
        for K in 2**N..2**(N+1) loop
           This := Float (S (K)) / Float (K);
           if Max < This then
              Max := This;
              Where := K;
           end if;
        end loop;
        if Max > 0.55 then
           Mallows := Where;
        end if;
        Put (Integer'Image (N) & " .." & Integer'Image (N + 1) & ":");
        Put_Line (Float'Image (Max) & " at" & Integer'Image (Where));
     end;
  end loop;
  Put_Line ("Mallows number" & Integer'Image (Mallows));

end Test_Hofstadter_Conway; </lang> Sample output:

 1 .. 2: 6.66667E-01 at 3
 2 .. 3: 6.66667E-01 at 6
 3 .. 4: 6.36364E-01 at 11
 4 .. 5: 6.08696E-01 at 23
 5 .. 6: 5.90909E-01 at 44
 6 .. 7: 5.76087E-01 at 92
 7 .. 8: 5.67416E-01 at 178
 8 .. 9: 5.59459E-01 at 370
 9 .. 10: 5.54937E-01 at 719
 10 .. 11: 5.50101E-01 at 1487
 11 .. 12: 5.47463E-01 at 2897
 12 .. 13: 5.44145E-01 at 5969
 13 .. 14: 5.42443E-01 at 11651
 14 .. 15: 5.40071E-01 at 22223
 15 .. 16: 5.38784E-01 at 45083
 16 .. 17: 5.37044E-01 at 89516
 17 .. 18: 5.36020E-01 at 181385
 18 .. 19: 5.34645E-01 at 353683
 19 .. 20: 5.33779E-01 at 722589
Mallows number 1487

ALGOL 68

Translation of: C

- note: This specimen retains the original C coding style.

Works with: ALGOL 68 version Revision 1 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

<lang algol68>PROC do sqnc = (INT max)INT: (

   [max]INT a list;
   INT v;
   INT k1 := 2;
   INT lg2 := 1;
   INT nmax; 
   LONG REAL amax := 0;
   INT out mallows number;
   v := a list[1] := a list[2] := 1;

   FOR n FROM 3 TO max DO
       v := a list[n] := a list[v] + a list[n-v];
       IF amax < v/n THEN 
         amax := v/n; 
         nmax := n ;
         IF v/n >= 0.55 THEN out mallows number := n FI
       FI;
       IF 2r0 = (BIN k1 AND BIN n) THEN
           printf(($"Maximum between 2^"g(0)" and 2^"g(0)" is about "g(-8,6)" at "g(0)l$, lg2,lg2+1, amax, nmax));
           amax := 0;
           lg2 +:= 1
       FI;
       k1 := n
   OD;
   out mallows number

);

INT mallows number = do sqnc(2**20);

printf(($"You too might have won $1000 with an answer of n = "g(0)$, mallows number))</lang> Output:

Maximum between 2^1 and 2^2 is about 0.666667 at 3
Maximum between 2^2 and 2^3 is about 0.666667 at 6
Maximum between 2^3 and 2^4 is about 0.636364 at 11
Maximum between 2^4 and 2^5 is about 0.608696 at 23
Maximum between 2^5 and 2^6 is about 0.590909 at 44
Maximum between 2^6 and 2^7 is about 0.576087 at 92
Maximum between 2^7 and 2^8 is about 0.567416 at 178
Maximum between 2^8 and 2^9 is about 0.559459 at 370
Maximum between 2^9 and 2^10 is about 0.554937 at 719
Maximum between 2^10 and 2^11 is about 0.550101 at 1487
Maximum between 2^11 and 2^12 is about 0.547463 at 2897
Maximum between 2^12 and 2^13 is about 0.544145 at 5969
Maximum between 2^13 and 2^14 is about 0.542443 at 11651
Maximum between 2^14 and 2^15 is about 0.540071 at 22223
Maximum between 2^15 and 2^16 is about 0.538784 at 45083
Maximum between 2^16 and 2^17 is about 0.537044 at 89516
Maximum between 2^17 and 2^18 is about 0.536020 at 181385
Maximum between 2^18 and 2^19 is about 0.534645 at 353683
Maximum between 2^19 and 2^20 is about 0.533779 at 722589
You too might have won $1000 with an answer of n = 1487

AutoHotkey

<lang autohotkey>Progress, b2 w150 zh0 fs9, CreateLists ... CreateLists(2 ** (Max:=20))

Progress,, Find Maxima ... Loop, % Max - 1

   msg .= "Maximum between 2^" A_Index " and 2^" A_Index + 1
       .  " is " GetMax(2 ** A_Index, 2 ** (A_Index + 1), n)
       .  " for n = " n "`n"

Progress,, Find Mallows Number ... Loop, % 2 ** Max

   If (n_%A_Index% > 0.55)
       MallowsNumber := A_Index

msg .= "Mallows Number = " MallowsNumber

Progress, Off MsgBox, %msg%

---------------------------------------------------------------------------

GetMax(a, b, ByRef Item) { ; return max value of a(n)/n between a and b

---------------------------------------------------------------------------
   Loop {
       IfGreater, a, %b%, Break
       If (Maximum < n_%a%)
           Maximum := n_%a%, Item := a
       a++
   }
   Return, Maximum

}

---------------------------------------------------------------------------

CreateLists(Lenght) { ; Hofstadter-Conway sequences (using lookups)

---------------------------------------------------------------------------
   ; create the sequence  a_%A_Index%  [ a(n)   ]
   ;   and  the sequence  n_%A_Index%  [ a(n)/n ]
   ;-----------------------------------------------------------------------
   global
   a_1 := a_2 := n_1 := 1, n_2 := 1 / 2
   Loop, %Lenght% {
       IfLess, A_Index, 3, Continue
       n1 := A_Index - 1
       an1 := a_%n1%
       nan1 := A_Index - an1
       a_%A_Index% := a_%an1% + a_%nan1%
       n_%A_Index% := a_%A_Index% / A_Index
   }

}</lang> Message box shows:

Maximum between 2^1 and 2^2 is 0.666667 for n = 3
Maximum between 2^2 and 2^3 is 0.666667 for n = 6
Maximum between 2^3 and 2^4 is 0.636364 for n = 11
Maximum between 2^4 and 2^5 is 0.608696 for n = 23
Maximum between 2^5 and 2^6 is 0.590909 for n = 44
Maximum between 2^6 and 2^7 is 0.576087 for n = 92
Maximum between 2^7 and 2^8 is 0.567416 for n = 178
Maximum between 2^8 and 2^9 is 0.559459 for n = 370
Maximum between 2^9 and 2^10 is 0.554937 for n = 719
Maximum between 2^10 and 2^11 is 0.550101 for n = 1487
Maximum between 2^11 and 2^12 is 0.547463 for n = 2897
Maximum between 2^12 and 2^13 is 0.544145 for n = 5969
Maximum between 2^13 and 2^14 is 0.542443 for n = 11651
Maximum between 2^14 and 2^15 is 0.540071 for n = 22223
Maximum between 2^15 and 2^16 is 0.538784 for n = 45083
Maximum between 2^16 and 2^17 is 0.537044 for n = 89516
Maximum between 2^17 and 2^18 is 0.536020 for n = 181385
Maximum between 2^18 and 2^19 is 0.534645 for n = 353683
Maximum between 2^19 and 2^20 is 0.533779 for n = 722589
Mallows Number = 1489

C

<lang C>#include <stdio.h>

  1. 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

C++

<lang cpp>

  1. include <deque>
  2. include <iostream>

int hcseq(int n) {

 static std::deque<int> seq(2, 1);
 while (seq.size() < n)
 {
   int x = seq.back();
   seq.push_back(seq[x-1] + seq[seq.size()-x]);
 }
 return seq[n-1];

}

int main() {

 int pow2 = 1;
 for (int i = 0; i < 20; ++i)
 {
   int pow2next = 2*pow2;
   double max = 0;
   for (int n = pow2; n < pow2next; ++n)
   {
     double anon = hcseq(n)/double(n);
     if (anon > max)
       max = anon;
   }
   std::cout << "maximum of a(n)/n between 2^" << i
             << " (" << pow2 << ") and 2^" << i+1
             << " (" << pow2next << ") is " << max << "\n";
   pow2 = pow2next;
 }

} </lang> Output:

maximum of a(n)/n between 2^0 (1) and 2^1 (2) is 1
maximum of a(n)/n between 2^1 (2) and 2^2 (4) is 0.666667
maximum of a(n)/n between 2^2 (4) and 2^3 (8) is 0.666667
maximum of a(n)/n between 2^3 (8) and 2^4 (16) is 0.636364
maximum of a(n)/n between 2^4 (16) and 2^5 (32) is 0.608696
maximum of a(n)/n between 2^5 (32) and 2^6 (64) is 0.590909
maximum of a(n)/n between 2^6 (64) and 2^7 (128) is 0.576087
maximum of a(n)/n between 2^7 (128) and 2^8 (256) is 0.567416
maximum of a(n)/n between 2^8 (256) and 2^9 (512) is 0.559459
maximum of a(n)/n between 2^9 (512) and 2^10 (1024) is 0.554937
maximum of a(n)/n between 2^10 (1024) and 2^11 (2048) is 0.550101
maximum of a(n)/n between 2^11 (2048) and 2^12 (4096) is 0.547463
maximum of a(n)/n between 2^12 (4096) and 2^13 (8192) is 0.544145
maximum of a(n)/n between 2^13 (8192) and 2^14 (16384) is 0.542443
maximum of a(n)/n between 2^14 (16384) and 2^15 (32768) is 0.540071
maximum of a(n)/n between 2^15 (32768) and 2^16 (65536) is 0.538784
maximum of a(n)/n between 2^16 (65536) and 2^17 (131072) is 0.537044
maximum of a(n)/n between 2^17 (131072) and 2^18 (262144) is 0.53602
maximum of a(n)/n between 2^18 (262144) and 2^19 (524288) is 0.534645
maximum of a(n)/n between 2^19 (524288) and 2^20 (1048576) is 0.533779

C#

Works with: C# version 3.0

<lang csharp> using System; using System.Linq;

namespace HofstadterConway {

   class Program
   {
       static int[] GenHofstadterConway(int max)
       {
           int[] result = new int[max];
           result[0]=result[1]=1;
           for (int ix = 2; ix < max; ix++)
               result[ix] = result[result[ix - 1] - 1] + result[ix - result[ix - 1]];
           return result;
       }
       static void Main(string[] args)
       {
           double[] adiv = new double[1 << 20];
           {
               int[] a = GenHofstadterConway(1 << 20);
               for (int i = 0; i < 1 << 20; i++)
                   adiv[i] = a[i] / (double)(i + 1);
           }
           for (int p = 2; p <= 20; p++)
           {
               var max = Enumerable.Range(
                    (1 << (p - 1)) - 1,
                    (1 << p) - (1 << (p - 1))
                    )
                    .Select(ix => new { I = ix + 1, A = adiv[ix] })
                    .OrderByDescending(x => x.A)
                    .First();
               Console.WriteLine("Maximum from 2^{0} to 2^{1} is {2} at {3}",
                   p - 1, p, max.A, max.I);
           }
           Console.WriteLine("The winning number is {0}.",
               Enumerable.Range(0, 1 << 20)
                   .Last(i => (adiv[i] > 0.55)) + 1
               );
       }
   }

} </lang> Output:-

Maximum from 2^1 to 2^2 is 0.666666666666667 at 3
Maximum from 2^2 to 2^3 is 0.666666666666667 at 6
Maximum from 2^3 to 2^4 is 0.636363636363636 at 11
Maximum from 2^4 to 2^5 is 0.608695652173913 at 23
Maximum from 2^5 to 2^6 is 0.590909090909091 at 44
Maximum from 2^6 to 2^7 is 0.576086956521739 at 92
Maximum from 2^7 to 2^8 is 0.567415730337079 at 178
Maximum from 2^8 to 2^9 is 0.559459459459459 at 370
Maximum from 2^9 to 2^10 is 0.554937413073713 at 719
Maximum from 2^10 to 2^11 is 0.550100874243443 at 1487
Maximum from 2^11 to 2^12 is 0.547462892647566 at 2897
Maximum from 2^12 to 2^13 is 0.544144747863964 at 5969
Maximum from 2^13 to 2^14 is 0.542442708780362 at 11651
Maximum from 2^14 to 2^15 is 0.540071097511587 at 22223
Maximum from 2^15 to 2^16 is 0.538784020584256 at 45083
Maximum from 2^16 to 2^17 is 0.537043656999866 at 89516
Maximum from 2^17 to 2^18 is 0.536020067811561 at 181385
Maximum from 2^18 to 2^19 is 0.534645431078112 at 353683
Maximum from 2^19 to 2^20 is 0.533779229963368 at 722589
The winning number is 1489.

Clojure

<lang Clojure>(use 'clojure.contrib.math) ; for expt

A literal transcription of the definition, with memoize doing the heavy lifting

(def conway

    (memoize
     (fn [x]
       (if (< x 3)
         1
         (+ (conway (conway (dec x)))
            (conway (- x (conway (dec x)))))))))

(def N (drop 1 (range))) ; natural numbers

This is enough to compute them all. The rest is grouping and max-finding

(def all-conways (map conway N))

All the powers of two

(def pow2 (map #(expt 2 %) N))

Find the lowest power of two > n

(defn lowest-pow-higher [n]

 (some #(and (> % n) %) pow2))
Split the natural numbers into groups at each power of two

(def groups (partition-by lowest-pow-higher N))

The conway numbers of each number in the group

(def C (map #(map conway %) groups))

Each conway number divided by its index

(def ratios (map #(map / %1 %2) C groups))

The largest value in each group of ratios

(def maxima (map #(apply max %) ratios))

(take 4 maxima) ; no rounding errors: still using ratios

yields (1 2/3 2/3 7/11)

(def first-20 (take 20 maxima)) (apply >= first-20) ; yields true - the sequence is decreasing (map double first-20) ; to get them in decimal form</lang>

Common Lisp

<lang lisp>(defparameter *hof-con*

 (make-array '(2) :initial-contents '(1 1) :adjustable t

:element-type 'integer :fill-pointer 2))

(defparameter *hof-con-ratios*

 (make-array '(2) :initial-contents '(1.0 0.5) :adjustable t

:element-type 'single-float :fill-pointer 2))

(defun hof-con (n)

 (let ((l (length *hof-con*)))
   (if (<= n l) (aref *hof-con* (1- n))

(extend-hof-con-sequence l n))))

(defun extend-hof-con-sequence (l n)

 (loop for i from l below n do
      (let* ((x (aref *hof-con* (1- i)))

(hc (+ (aref *hof-con* (1- x)) (aref *hof-con* (- i x))))) (vector-push-extend hc *hof-con*) (vector-push-extend (/ hc (+ i 1.0)) *hof-con-ratios*)))

 (aref *hof-con* (1- n)))

(defun max-in-array-range (arr id1 id2)

 (let ((m 0) (id 0))
   (loop for i from (1- id1) to (1- id2) do

(let ((n (aref arr i))) (if (> n m) (setq m n id i))))

   (values m (1+ id))))

(defun maxima (po2)

 (hof-con (expt 2 po2))
 (loop for i from 1 below po2 do
      (let ((id1 (expt 2 i)) (id2 (expt 2 (1+ i))))

(multiple-value-bind (m id) (max-in-array-range *hof-con-ratios* id1 id2) (format t "Local maximum in [~A .. ~A]: ~A at n = ~A~%" id1 id2 m id)))))

(defun mallows (po2)

 (let ((n (expt 2 po2)))
   (hof-con n)
   (do ((i (1- n) (1- i)))

((> (aref *hof-con-ratios* i) 0.55) (+ i 2)))))</lang> Sample session:

ROSETTA> (maxima 20)
Local maximum in [2 .. 4]: 0.6666667 at n = 3
Local maximum in [4 .. 8]: 0.6666667 at n = 6
Local maximum in [8 .. 16]: 0.6363636 at n = 11
Local maximum in [16 .. 32]: 0.6086956 at n = 23
Local maximum in [32 .. 64]: 0.59090906 at n = 44
Local maximum in [64 .. 128]: 0.57608694 at n = 92
Local maximum in [128 .. 256]: 0.5674157 at n = 178
Local maximum in [256 .. 512]: 0.55945945 at n = 370
Local maximum in [512 .. 1024]: 0.5549374 at n = 719
Local maximum in [1024 .. 2048]: 0.55010086 at n = 1487
Local maximum in [2048 .. 4096]: 0.5474629 at n = 2897
Local maximum in [4096 .. 8192]: 0.54414475 at n = 5969
Local maximum in [8192 .. 16384]: 0.5424427 at n = 11651
Local maximum in [16384 .. 32768]: 0.54007107 at n = 22223
Local maximum in [32768 .. 65536]: 0.538784 at n = 45083
Local maximum in [65536 .. 131072]: 0.53704363 at n = 89516
Local maximum in [131072 .. 262144]: 0.53602004 at n = 181385
Local maximum in [262144 .. 524288]: 0.53464544 at n = 353683
Local maximum in [524288 .. 1048576]: 0.5337792 at n = 722589
NIL
ROSETTA> (mallows 20)
1490

D

D v.2 <lang D>import std.stdio: writefln;

void hofstadterConwaySequence(int m) {

   auto alist = new int[m + 1];
   alist[0 .. 2] = 1;
   double amax = 0.0;
   auto v = alist[2];
   int k1 = 2;
   int lg2 = 1;
   foreach (n; 2 .. m + 1) {
       v = alist[n] = alist[v] + alist[n - v];
       if (amax < v * 1.0 / n)
           amax = v * 1.0 / n;
       if ((k1 & n) == 0) {
           writefln("Maximum between 2^%d and 2^%d was %f", lg2, lg2 + 1, amax);
           amax = 0;
           lg2++;
       }
       k1 = n;
   }

}

void main() {

   hofstadterConwaySequence(1 << 20);

}</lang>

Output:

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^5 and 2^6 was 0.590909
Maximum between 2^6 and 2^7 was 0.576087
Maximum between 2^7 and 2^8 was 0.567416
Maximum between 2^8 and 2^9 was 0.559459
Maximum between 2^9 and 2^10 was 0.554937
Maximum between 2^10 and 2^11 was 0.550101
Maximum between 2^11 and 2^12 was 0.547463
Maximum between 2^12 and 2^13 was 0.544145
Maximum between 2^13 and 2^14 was 0.542443
Maximum between 2^14 and 2^15 was 0.540071
Maximum between 2^15 and 2^16 was 0.538784
Maximum between 2^16 and 2^17 was 0.537044
Maximum between 2^17 and 2^18 was 0.536020
Maximum between 2^18 and 2^19 was 0.534645
Maximum between 2^19 and 2^20 was 0.533779

Fortran

<lang Fortran> program conway

 implicit none
 integer :: a(2**20)  ! The sequence a(n)
 real    :: b(2**20)  ! The sequence a(n)/n
 real    :: v         ! Max value in the range [2*i, 2**(i+1)]
 integer :: nl(1)     ! The location of v in the array b(n)
 integer :: i, N, first, second, last, m
 ! Populate a(n) and b(n)
 a(1:2) = [1, 1]
 b(1:2) = [1.0e0, 0.5e0]
 N = 2
 do i=1,2**20
    last = a(N)
    first = a(last)
    second = a(N-last+1)
    N = N+1
    a(N:N) = first + second
    b(N:N) = a(N:N)/real(N)
 end do
 ! Calculate the max values in the logarithmic ranges
 m = 0
 do i=1,19
    v = maxval(b(2**i:2**(i+1)))
    nl = maxloc(b(2**i:2**(i+1)))
    write(*,'(2(a,i0),a,f8.6,a,i0)') &
         'Max. between 2**', i,      &
         ' and 2**', (i+1),          &
         ' is ', v,                  &
         ' at n = ', 2**i+nl(1)-1
    if (m == 0 .and. v < 0.55e0) then
       m = i-1
    end if
 end do
 ! Calculate Mallows number
 do i=2**(m+1), 2**m,-1
    if (b(i) > 0.55e0) then
       exit
    end if
 end do
 write(*,'(a,i0)') 'Mallows number = ',i

end program conway </lang>

Output:

Max. between 2**1 and 2**2 is 0.666667 at n = 3
Max. between 2**2 and 2**3 is 0.666667 at n = 6
Max. between 2**3 and 2**4 is 0.636364 at n = 11
Max. between 2**4 and 2**5 is 0.608696 at n = 23
Max. between 2**5 and 2**6 is 0.590909 at n = 44
Max. between 2**6 and 2**7 is 0.576087 at n = 92
Max. between 2**7 and 2**8 is 0.567416 at n = 178
Max. between 2**8 and 2**9 is 0.559459 at n = 370
Max. between 2**9 and 2**10 is 0.554937 at n = 719
Max. between 2**10 and 2**11 is 0.550101 at n = 1487
Max. between 2**11 and 2**12 is 0.547463 at n = 2897
Max. between 2**12 and 2**13 is 0.544145 at n = 5969
Max. between 2**13 and 2**14 is 0.542443 at n = 11651
Max. between 2**14 and 2**15 is 0.540071 at n = 22223
Max. between 2**15 and 2**16 is 0.538784 at n = 45083
Max. between 2**16 and 2**17 is 0.537044 at n = 89516
Max. between 2**17 and 2**18 is 0.536020 at n = 181385
Max. between 2**18 and 2**19 is 0.534645 at n = 353683
Max. between 2**19 and 2**20 is 0.533779 at n = 722589
Mallows number = 1489

F#

<lang fsharp>let a = ResizeArray[0; 1; 1] while a.Count <= (1 <<< 20) do

 a.[a.[a.Count - 1]] + a.[a.Count - a.[a.Count - 1]] |> a.Add

for p = 1 to 19 do

 Seq.max [|for i in 1 <<< p .. 1 <<< p+1 -> float a.[i] / float i|]
 |> printf "Maximum in %6d..%7d is %g\n" (1 <<< p) (1 <<< p+1)</lang>

Outputs:

Maximum in      2..      4 is 0.666667
Maximum in      4..      8 is 0.666667
Maximum in      8..     16 is 0.636364
Maximum in     16..     32 is 0.608696
Maximum in     32..     64 is 0.590909
Maximum in     64..    128 is 0.576087
Maximum in    128..    256 is 0.567416
Maximum in    256..    512 is 0.559459
Maximum in    512..   1024 is 0.554937
Maximum in   1024..   2048 is 0.550101
Maximum in   2048..   4096 is 0.547463
Maximum in   4096..   8192 is 0.544145
Maximum in   8192..  16384 is 0.542443
Maximum in  16384..  32768 is 0.540071
Maximum in  32768..  65536 is 0.538784
Maximum in  65536.. 131072 is 0.537044
Maximum in 131072.. 262144 is 0.53602
Maximum in 262144.. 524288 is 0.534645
Maximum in 524288..1048576 is 0.533779

Haskell

<lang haskell>import Data.List import Data.Ord import Data.Array import Text.Printf

hc :: Int -> Array Int Int hc n = arr

 where arr = array (1, n) $ (1, 1) : (2, 1) : map (f (arr!)) [3 .. n]
       f a i = (i, a (a $ i - 1) + a (i - a (i - 1)))

printMaxima :: (Int, (Int, Double)) -> IO () printMaxima (n, (pos, m)) =

   printf "Max between 2^%-2d and 2^%-2d is %1.5f at n = %6d\n"
                            n    (n + 1)        m        pos

main = do

   mapM_ printMaxima maxima
   printf "Mallows's number is %d\n" mallows
 where
   hca = hc $ 2^20
   hc' n  = fromIntegral (hca!n) / fromIntegral n
   maxima = zip [0..] $ map max powers
   max seq = maximumBy (comparing snd) $ zip seq (map hc' seq)
   powers = map (\n -> [2^n .. 2^(n + 1) - 1]) [0 .. 19]
   mallows = last.takeWhile ((< 0.55) . hc') $ [2^20, 2^20 - 1 .. 1]</lang>

Icon and Unicon

Icon

<lang icon>procedure main(args)

   m       := integer(!args) | 20
   nextNum := create put(A := [], 1 | 1 | |A[A[*A]]+A[-A[*A]])[*A]
   p2      := 2 ^ (p := 1)
   maxv    := 0
   every n := 1 to (2^m) do {
       if maxv <:= (x := @nextNum / real(n)) then maxm := n
       if x >= 0.55 then mallow := n+1   # Want *next* n!
       if n = p2 then {
           write("Max between 2^",p-1," and 2^",p," is ",maxv," at n = ",maxm)
           p2 := 2 ^ (p +:= 1)
           maxv := 0
           }
       }
   write("Mallow's number is ",\mallow | "NOT found!")

end</lang>

Output:

->hc
Max between 2^0 and 2^1 is 1.0 at n = 1
Max between 2^1 and 2^2 is 0.6666666667 at n = 3
Max between 2^2 and 2^3 is 0.6666666667 at n = 6
Max between 2^3 and 2^4 is 0.6363636364 at n = 11
Max between 2^4 and 2^5 is 0.6086956522 at n = 23
Max between 2^5 and 2^6 is 0.5909090909 at n = 44
Max between 2^6 and 2^7 is 0.5760869565 at n = 92
Max between 2^7 and 2^8 is 0.5674157303 at n = 178
Max between 2^8 and 2^9 is 0.5594594595 at n = 370
Max between 2^9 and 2^10 is 0.5549374131 at n = 719
Max between 2^10 and 2^11 is 0.5501008742 at n = 1487
Max between 2^11 and 2^12 is 0.5474628926 at n = 2897
Max between 2^12 and 2^13 is 0.5441447479 at n = 5969
Max between 2^13 and 2^14 is 0.5424427088 at n = 11651
Max between 2^14 and 2^15 is 0.5400710975 at n = 22223
Max between 2^15 and 2^16 is 0.5387840206 at n = 45083
Max between 2^16 and 2^17 is 0.537043657 at n = 89516
Max between 2^17 and 2^18 is 0.5360200678 at n = 181385
Max between 2^18 and 2^19 is 0.5346454311 at n = 353683
Max between 2^19 and 2^20 is 0.53377923 at n = 722589
Mallow's number is 1490
->

Unicon

The Icon solution also works in Unicon.

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>

Java

Translation of: C

with corrections to 0 indexing.

<lang java>public class HofCon {

   public static void main(String[] args){
       doSqnc(1<<20);
   }
   public static void doSqnc(int m) {
       int[] a_list = new int[m + 1];
       int max_df = 0;
       int p2_max = 2;
       int v, n;
       int k1 = 2;
       int lg2 = 1;
       double amax = 0;
       a_list[0] = a_list[1] = 1;
       v = a_list[2];
       for (n = 2; 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)) {
               System.out.printf("Maximum between 2^%d and 2^%d was %f\n", lg2, lg2 + 1, amax);
               amax = 0;
               lg2++;
           }
           k1 = n;
       }
   }

}</lang> Output:

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

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

Perl 6

Works with: Rakudo version 2010.09

This is written in relatively low-level primitives (for Perl 6) because rakudo still has some performance issues with long lists. One interesting feature, however, is the use of the max infix operator as an assignment operator, and using it on a Pair so that the maximum is calculated on the key of the Pair, while the value of the Pair carries the index of the maximum through to the end without having to keep a second variable. <lang perl6>my $POW = 20; my $top = 2**$POW;

my @a = (0,1,1); @a[$top] = 0; # pre-extend array

my $n = 3; my $p = 1;

loop ($n = 3; $n <= $top; $n++) {

   @a[$n] = $p = @a[$p] + @a[$n - $p];

}

my $last55; for 1 ..^ $POW -> $power {

   my $beg = 2 ** $power;
   my $end = $beg * 2 - 1;
   my $max;
   my $ratio;
   loop (my $n = $beg; $n <= $end; $n++) {
       my $ratio = @a[$n] / $n;
       $last55 = $n if $ratio * 100 >= 55;
       $max max= $ratio => $n;
   }
   say $power.fmt('%2d'), $beg.fmt("%10d"), '..', $end.fmt("%-10d"), $max.key, " at ", $max.value;

} say "Mallow's number would appear to be ", $last55+1;</lang>

Output:

<lang> 1 2..3 0.666666666666667 at 3

2         4..7         0.666666666666667 at 6
3         8..15        0.636363636363636 at 11
4        16..31        0.608695652173913 at 23
5        32..63        0.590909090909091 at 44
6        64..127       0.576086956521739 at 92
7       128..255       0.567415730337079 at 178
8       256..511       0.559459459459459 at 370
9       512..1023      0.554937413073713 at 719

10 1024..2047 0.550100874243443 at 1487 11 2048..4095 0.547462892647566 at 2897 12 4096..8191 0.544144747863964 at 5969 13 8192..16383 0.542442708780362 at 11651 14 16384..32767 0.540071097511587 at 22223 15 32768..65535 0.538784020584256 at 45083 16 65536..131071 0.537043656999866 at 89516 17 131072..262143 0.536020067811561 at 181385 18 262144..524287 0.534645431078112 at 353683 19 524288..1048575 0.533779229963368 at 722589 Mallow's number would appear to be 1490</lang> Incidentally, this takes three hours or so currently. We expect this to get much faster.

Here is a more list-oriented version. Note that @a is a lazy array, and the Z variants are "zipwith" operators. <lang perl6>my $n = 3; my @a := (0,1,1, -> $p { @a[$p] + @a[$n++ - $p] } ... *);

my $last55; for 1..19 -> $power {

   my @range := 2**$power .. 2**($power+1)-1;
   my @ratios = (@a[@range] Z/ @range) Z=> @range;
   my $max = [max] @ratios;
   ($last55 = .value if .key >= .55 for @ratios) if $max.key >= .55;
   say $power.fmt('%2d'), @range.min.fmt("%10d"), '..', @range.max.fmt("%-10d"),
       $max.key, ' at ', $max.value;

} say "Mallow's number would appear to be ", $last55+1;</lang>

PicoLisp

<lang PicoLisp>(de hofcon (N)

  (cache '*HofCon (format (seed N))
     (if (>= 2 N)
        1
        (+
           (hofcon (hofcon (dec N)))
           (hofcon (- N (hofcon (dec N)))) ) ) ) )

(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')


PL/I

<lang PL/I> /* First part: */

declare L (10000) fixed static initial ((1000) 0); L(1), L(2) = 1; do i = 3 to 10000;

  k = L(i);
  L(i) = L(i-k) + L(1+k);

end; </lang>

PureBasic

<lang PureBasic>If OpenConsole()

 Define.i upperlim, i=1, k1=2, n=3, v=1
 Define.d Maximum
 Print("Enter limit (ENTER gives 2^20=1048576): "): upperlim=Val(Input())
 If upperlim<=0: upperlim=1048576: EndIf
 Dim tal(upperlim)
 If ArraySize(tal())=-1
   PrintN("Could not allocate needed memory!"): Input(): End
 EndIf
 tal(1)=1: tal(2)=1
 While n<=upperlim
   v=tal(v)+tal(n-v)
   tal(n)=v
   If Maximum<(v/n): Maximum=v/n: EndIf
   If Not n&k1
     PrintN("Maximum between 2^"+Str(i)+" and 2^"+Str(i+1)+" was "+StrD(Maximum,6))
     Maximum=0.0
     i+1
   EndIf
   k1=n
   n+1
 Wend
 
 Print(#CRLF$+"Press ENTER to exit."): Input()
 CloseConsole()

EndIf</lang>

Python

<lang python>from __future__ import division

def maxandmallows(nmaxpower2=20):

   # Note: The first hc number is returned in hc[1];
   # hc[0] is not part of the series.
   nmax = 2**nmaxpower2
   hc, mallows, mx, mxpow2 = [None, 1, 1], None, (0.5, 2), []
   for n in range(2, nmax + 1):
       ratio = hc[n] / n
       if ratio > mx[0]: mx = (ratio, n)
       if ratio >= 0.55: mallows = n + 1
       if ratio == 0.5:
           print("In the region %7i < n <= %7i: max a(n)/n = %s" %

((n)//2, n, "%6.4f at n = %i" % mx))

           mxpow2.append(mx[0])
           mx = (ratio, n)
       hc.append(hc[hc[n]] + hc[-hc[n]])
   return hc, mallows if mxpow2 and mxpow2[-1] < 0.55 and n > 4 else None

if __name__ == '__main__':

   hc, mallows = maxandmallows(20)
   if mallows:
       print("\nYou too might have won $1000 with the mallows number of %i" % mallows)

</lang>

Sample output

In the region       1 < n <=       2: max a(n)/n = 0.5000 at  n = 2
In the region       2 < n <=       4: max a(n)/n = 0.6667 at  n = 3
In the region       4 < n <=       8: max a(n)/n = 0.6667 at  n = 6
In the region       8 < n <=      16: max a(n)/n = 0.6364 at  n = 11
In the region      16 < n <=      32: max a(n)/n = 0.6087 at  n = 23
In the region      32 < n <=      64: max a(n)/n = 0.5909 at  n = 44
In the region      64 < n <=     128: max a(n)/n = 0.5761 at  n = 92
In the region     128 < n <=     256: max a(n)/n = 0.5674 at  n = 178
In the region     256 < n <=     512: max a(n)/n = 0.5595 at  n = 370
In the region     512 < n <=    1024: max a(n)/n = 0.5549 at  n = 719
In the region    1024 < n <=    2048: max a(n)/n = 0.5501 at  n = 1487
In the region    2048 < n <=    4096: max a(n)/n = 0.5475 at  n = 2897
In the region    4096 < n <=    8192: max a(n)/n = 0.5441 at  n = 5969
In the region    8192 < n <=   16384: max a(n)/n = 0.5424 at  n = 11651
In the region   16384 < n <=   32768: max a(n)/n = 0.5401 at  n = 22223
In the region   32768 < n <=   65536: max a(n)/n = 0.5388 at  n = 45083
In the region   65536 < n <=  131072: max a(n)/n = 0.5370 at  n = 89516
In the region  131072 < n <=  262144: max a(n)/n = 0.5360 at  n = 181385
In the region  262144 < n <=  524288: max a(n)/n = 0.5346 at  n = 353683
In the region  524288 < n <= 1048576: max a(n)/n = 0.5338 at  n = 722589

You too might have won $1000 with the mallows number of 1490

If you don't create enough terms in the sequence, no mallows number is returned.

Ruby

<lang ruby>class HofstadterConway10000

 def initialize
   @sequence = [nil, 1, 1]
 end
 attr_reader :sequence
 def [](n)
   raise ArgumentError, "n must be >= 1" if n < 1
   a = @sequence
   a.length.upto(n) {|i| a[i] = a[a[i-1]] + a[i-a[i-1]] }
   a[n]
 end

end

hc = HofstadterConway10000.new</lang>

I run the sequence backwards to make it easier to find the Mallows number.

The first call to the hc.[] method will pass the value 2^20, pre-populating the entire sequence.

<lang ruby>p = -1 20.downto(1).each_cons(2) do |j, i|

 max_n, max_v = -1, -1
 (2**j).downto(2**i).each do |n|
   v = hc[n].to_f / n
   max_n, max_v = n, v if v > max_v
   # Mallows number
   p = n if p == -1 and v >= 0.55
 end
 puts "maximum between 2^#{i} and 2^#{j} occurs at #{max_n}: #{max_v}"

end

puts "the mallows number is #{p}"</lang>

output

maximum between 2^19 and 2^20 occurs at 722589: 0.533779229963368
maximum between 2^18 and 2^19 occurs at 353683: 0.534645431078112
maximum between 2^17 and 2^18 occurs at 181385: 0.536020067811561
maximum between 2^16 and 2^17 occurs at 89516: 0.537043656999866
maximum between 2^15 and 2^16 occurs at 45083: 0.538784020584256
maximum between 2^14 and 2^15 occurs at 22223: 0.540071097511587
maximum between 2^13 and 2^14 occurs at 11651: 0.542442708780362
maximum between 2^12 and 2^13 occurs at 5969: 0.544144747863964
maximum between 2^11 and 2^12 occurs at 2897: 0.547462892647566
maximum between 2^10 and 2^11 occurs at 1487: 0.550100874243443
maximum between 2^9 and 2^10 occurs at 719: 0.554937413073713
maximum between 2^8 and 2^9 occurs at 370: 0.559459459459459
maximum between 2^7 and 2^8 occurs at 178: 0.567415730337079
maximum between 2^6 and 2^7 occurs at 92: 0.576086956521739
maximum between 2^5 and 2^6 occurs at 44: 0.590909090909091
maximum between 2^4 and 2^5 occurs at 23: 0.608695652173913
maximum between 2^3 and 2^4 occurs at 11: 0.636363636363636
maximum between 2^2 and 2^3 occurs at 6: 0.666666666666667
maximum between 2^1 and 2^2 occurs at 3: 0.666666666666667
the mallows number is 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