(29C) Prime numbers up to 10'000

[Thomas Klemm]

(09-04-2018 12:04 AM)Archilog Wrote:  A guy called C.Ret did an awesome work which can be found there: Silicium

Indeed it's amazing.

Here's the flow-chart:


And here's the program for the HP-19C:

Quote:Pour imprimer les nombres premiers jusqu'à m (m compris entre de 50 et11449):
* saisir le programme (97 pas),
* vérifier l'imprimante, éventuellement y mettre un rouleau de papier neuf,
* initialiser les registres mémoires R1, R2 , R3 et R4 avec les valeurs suivantes :
5 STO1 m STO2 7 STO1 et 200 STO4
* lancer le programme :
GSB0

This took me a while to figure out.
Here's how to initialise the registers, say for m = 10,000.

4 STO 1 ; 4 instead of 5
10000 STO 2
7 STO 3 ; STO 3 instead of STO 1
200 STO 4

Otherwise register 1 is incremented in line 88 to 6 and then the heap starts there.
But in line 29 it's hard-coded to use register 5.

I'm currently running both programs side by side in the aforementioned emulator.
His program is happily outrunning my program.

Here's a listing for those who want to run it:

Code:

01: 15 13 00 :   LBL 0      ; Minimal initialization
02:       02 :   2          ; 2
03:    14 74 :   PAUSE      ; » 2 «
04:       03 :   3          ; 3
05:    14 74 :   PAUSE      ; » 3 «
06:       51 :   +          ; 5
07:    14 74 :   PAUSE      ; » 5 «
08:    12 08 :   GSB 8      ; 7 is prime
09: 15 13 01 :   LBL 1      ; Main Loop
10:    12 02 :   GSB 2      ; 4 2
11:    12 02 :   GSB 2      ; 4 2
12:       04 :   4          ; 
13:    12 04 :   GSB 4      ; 4
14:       06 :   6          ; 
15:    12 04 :   GSB 4      ; 6
16:    12 03 :   GSB 3      ; 2
17:       06 :   6          ; 
18:    12 04 :   GSB 4      ; 6
19:    13 01 :   GTO 1      ; Main Loop
20: 15 13 02 :   LBL 2      ; Advance pseudo-prime list
21:       04 :   4          ; 
22:    12 04 :   GSB 4      ; 
23: 15 13 03 :   LBL 3      ; 
24:       02 :   2          ; 
25: 15 13 04 :   LBL 4      ; d
26: 23 51 03 :   STO+ 3     ; n ← n+d
27: 15 13 05 :   LBL 5      ; Test n
28:    24 03 :   RCL 3      ; n 
29:    24 05 :   RCL 5      ; H.p_1 n
30:    14 62 :   INT        ; H_1   n
31:    14 71 :   x=y        ; H_1 = n ?
32:    15 12 :   RTN        ; n is composite
33:    14 51 :   x>y        ; H_1 > n ?
34:    13 08 :   GTO 8      ; n is prime
35:    14 73 :   LSTx       ; H.p_1
36:    15 62 :   FRAC       ; .p_1
37:    24 04 :   RCL 4      ; 200       .p_1
38:       61 :   ×          ; 2p_1
39: 23 51 05 :   STO+ 5     ; H_1 ← H_1 + 2p_1
40:       05 :   5          ; 5
41: 15 13 06 :   LBL 6      ; « min-HEAP » loop
42:    23 00 :   STO 0      ; i ← x
43:    24 22 :   RCL i      ; H_i       i
44:       02 :   2          ; 2         H_i       i
45: 23 41 00 :   STO- 0     ; i ← i-2   H_i       i
46: 23 61 00 :   STO× 0     ; j ← 2i    H_i       i
47:       22 :   R↓         ; H_i       i
48:    24 00 :   RCL 0      ; j         H_i       i
49:    24 01 :   RCL 1      ; k         j         H_i       i
50:       41 :   -          ; j-k       H_i       i
51:    15 71 :   x=0        ; k = j ?
52:    13 07 :   GTO 7      ; 0         H_i       i
53:    15 51 :   x>0        ; k < j ?
54:    13 05 :   GTO 5      ; Heap OK
55:       22 :   R↓         ; H_i       i
56:    24 22 :   RCL i      ; H_j       H_i       i
57:    15 24 :   ISZ        ; j ← j+1
58:    24 22 :   RCL i      ; H_j+1     H_j       H_i       i
59:       41 :   -          ; H_j-H_j+1 H_i       i
60:    15 41 :   x<0        ; H_j < H_j+1 ?
61:    15 23 :   DSZ        ; j ← j-1
62: 15 13 07 :   LBL 7      ; _         H_i       i
63:       22 :   R↓         ; H_i       i
64:    24 22 :   RCL i      ; H_j       H_i       i
65:    14 51 :   x>y        ; H_j > H_i ?
66:    13 05 :   GTO 5      ; Exit only
67:    24 00 :   RCL 0      ; j         H_j       H_i       i
68:       21 :   x<>y       ; H_j       j         H_i       i
69:       22 :   R↓         ; j         H_i       i         H_j
70:       22 :   R↓         ; H_i       i         H_j       j
71:    23 22 :   STO i      ; H_j ← H_i
72:       22 :   R↓         ; i         H_j       j         H_i
73:    23 00 :   STO 0      ; i ← i
74:       22 :   R↓         ; H_j       j         H_i       i
75:    23 22 :   STO i      ; H_i ← H_j
76:       21 :   x<>y       ; j         H_j       H_i       i
77:    13 06 :   GTO 6      ; « min-HEAP » loop
78: 15 13 08 :   LBL 8      ; Is Prime - Store square and factor in HEAP
79:    24 02 :   RCL 2      ; m
80:    24 03 :   RCL 3      ; n         m
81:    14 51 :   x>y        ; n > m ?
82:       74 :   R/S        ; End of PRGM
83:    14 74 :   PAUSE      ; » n «
84:    15 63 :   x²         ; n²        m
85:    14 51 :   x>y        ; n² > m ?
86:    15 12 :   RTN        ; 
87:       01 :   1          ; 1         n²
88: 23 51 01 :   STO+ 1     ; k ← k+1
89:    24 01 :   RCL 1      ; k         1         n²
90:    23 00 :   STO 0      ; i ← k
91:       22 :   R↓         ; 1         n²
92:    14 73 :   LSTx       ; n         1         n²
93:       71 :   ÷          ; 1/n       n²
94:    15 21 :   %          ; n/100     n²
95:       51 :   +          ; n².n
96:    23 22 :   STO i      ; H.p_k ← n².n
97:    15 12 :   RTN        ;

What's so cool about it?
Instead of checking if a number is divisible by a prime the odd multiples of the primes are calculated starting from the square of the prime.

For 7 that would be: 49, 63, 77, 91, 105, …
Or for 11 that would be: 121, 143, 165, 187, …

Of course we can't keep all of them since we have only 30 registers.
But we don't have to. All we need is keeping the smallest among them.

These numbers are kept in a min-heap and thus we only have to check the smallest of them which is in register 5.

If the number is smaller it is a prime.
If they are the same it's not a prime and we can proceed with the next number.
If it is bigger we have to update the multiples of primes in the min-heap.

And while doing so make sure it's still a min-heap.
Thus some of the elements have to be rearranged.

The prime is added as a decimal value to the multiples:
For 7 that would be: 49.07, 63.07, 77.07, 91.07, 105.07, …
This allows to calculate the next number.

We can start with 7 since multiples of 2, 3 and 5 are omitted from the list of numbers.
The main loop with the multiple GSB commands generates sequence A007775 (apart from 1 of course).

I don't have a full understanding of the code for the min-heap yet.
Thus if C.Ret is reading this post I'm eager to get some explanations from the master himself.

Kind regards
Thomas

PS: Meanwhile we're at 1777 vs. 2137. The 2nd is of course C.Ret's program.

---

Edit: Added comments and stack diagrams to the listing of the program.