The Museum of HP Calculators

The Golomb Sequence for the HP-41

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   a_m = the greatest integer n such that G( n ) = m
          b_m = ------------------------------   G(G( n )) = m
and    c_m = ------------------------------    G(G(G( n ))) = m

-We have a₁ = b₁ = c₁ = 1 and a_m = a_m-1 + g_m ; b_m = b_m-1 + m.g_m ; c_m = c_m-1 + m.g_m( a_m-( g_m-1 )/2 )
and let cc_m,k defined by cc_m,0 = c_m ; cc_m,k+1 = cc_m,k + (m+1)(a_m+k+1)

-First, we find the greatest integer m such that c_m < n
-Then, -------------------------- k -------- cc_m,k < n
-Finally, G(n) = b_m + k(m+1) + INT(( n - cc_m,k + a_m +k )/( a_m + k +1 ))

-For n = 10¹⁰ , m = 330 and fortunately, a very simple formula yields G(m) for m < 423 i-e G(m) = the nearest integer to 1.2 m^0.61832

Data Registers: R00 = n ; R01 = a_m and a_m + k + 1 ; R02 = b_m and b_m + k(m+1)(k+1) ; R03 = g_m
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

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall