hp41programs

Golomb The Golomb Sequence for the HP-41

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

STACK	INPUTS	OUTPUTS
X	n	G(n)

Examples:

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

References:

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