The Museum of HP Calculators

# The Golomb Sequence for the HP-41

This program is Copyright © 2004 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

## Overview

-The Golomb sequence is defined by  G(n) = the number of times n appears , assuming  G(1) = 1
-"GOL" uses the recursive formula G(n) = 1 + G(n-G(G(n-1)))   ( n>1 )  to calculate G(n). The first values are:

 n 1 2 3 4 5 6 7 8 9 10 11 12 G(n) 1 2 2 3 3 4 4 4 5 5 5 6

-However, since the HP-41 has only 319 registers, another approach is needed to obtain G(n) for n > 318
and "GOL2" can compute up to G(9,999,999,999) in a "reasonable" time with only 8 data registers.

1°) Recursive Formula

Data Registers:             R00 = 0 ; R01 = G(1) = 1 ; ............ ; Rnn = G(n)
Flags: /
Subroutines: /

01  LBL "GOL"
02  .1
03  %
04  0
05  STO 00
06  SIGN
07  STO 01
08  ST+ Y
09  LBL 01
10  RDN
11  STO Y
12  RCL IND T
13  -
14  X<>Y
15  RCL IND Y
16  RCL 01
17  +
18  STO IND Y
19  ISG Y
20  GTO 01
21  END

( 37 bytes / SIZE nnn+1 )

 STACK INPUTS OUTPUTS Y / 1+n.nnn X n G(n) L / 1

Example:     41  XEQ "GOL"  yields  G(41) = 12  ( in 14 seconds )  and G(1) to G(41) are in registers R01 to R41

2°)  Up to  G(9,999,999,999) = 1,820,598

-"GOL2" uses the following method:

-Let   am  = the greatest integer n such that  G( n ) = m
bm = ------------------------------   G(G( n )) = m
and    cm = ------------------------------    G(G(G( n ))) = m

-We have   a1 = b1 = c1 = 1 and   am = am-1 + gm ;   bm = bm-1 + m.gm ;  cm = cm-1 + m.gm( am-( gm-1 )/2 )
and let  ccm,k  defined by   ccm,0 =  cm  ;    ccm,k+1 =  ccm,k + (m+1)(am+k+1)

-First, we find the greatest integer m such that  cm < n
-Then, --------------------------  k  --------   ccm,k < n
-Finally,  G(n) = bm + k(m+1) + INT(( n - ccm,k + am +k )/( am + k +1 ))

-For n = 1010 , m = 330 and fortunately, a very simple formula yields G(m) for m < 423   i-e  G(m) = the nearest integer to 1.2 m0.61832

Data Registers:   R00 = n  ;   R01 = am and  am + k + 1  ;  R02 = bm and  bm + k(m+1)(k+1)  ; R03 = gm
R04 = m+1  ;  R05 = 0.61832  ;  R06 = 1.2  ;  R07 = 0.5
Flags: /
Subroutines: /

01  LBL "GOL2"
02  STO 00
03  .61832
04  STO 05
05  1.2
06  STO 06
07  RCL 00
08  .5
09  STO 07
10  SIGN
11  STO 01
12  STO 02
13  STO 03
14  STO 04
15  -
16  FIX 0
17  LBL 01
18  RCL 01
19  RCL 04
20  ENTER^
21  SIGN
22  +
23  STO 04
24  RCL 05
25  Y^X
26  RCL 06
27  *
28  RND
29  STO 03
30  RCL 07
31  ST* Y
32  +
33  +
34  RCL 03
35  ST+ 01
36  RCL 04
37  *
38  ST+ 02
39  *
40  -
41  X>0?
42  GTO 01
43  RCL 03
44  ST- 01
45  RCL 04
46  *
47  ST- 02
48  R^
49  LBL 02
50  STO Y
51  RCL 01
52  1
53  +
54  STO 01
55  RCL 04
56  ST+ 02
57  *
58  -
59  X>0?
60  GTO 02
61  CLX
62  RCL 01
63  ST+ Y
64  DSE Y
65  /
66  INT
67  RCL 02
68  +
69  RCL 04
70  -
71  FIX 4
72  END

( 102 bytes / SIZE 008 )

Examples:    9,999,999,999  XEQ "GOL2"  >>>>  1,820,598   ( in 6mn22s )
1,000,000          R/S          >>>>         6,137   ( in 43seconds )

References:    Guy Toublanc & Th. Alle Zoethout in "48SXTANT" issues#39&40
N.J.A. Sloane's On-Line Encyclopedia of Integer Sequences:      www.research.att.com/~njas/sequences