+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (29C) Prime numbers up to 10'000 (/thread-11319.html)
(29C) Prime numbers up to 10'000 - Thomas Klemm - 09-02-2018 11:01 AM
Inspired by Jürgen's video printing the Prime numbers up to 10'000 with an HP-19C I took a closer look at its sibling the HP-29C.
It has 30 registers which is enough to store the primes up to 100.
We can use indirect addressing with register 0 so I thought I give it a try.
The first part (lines 01-35) calculates the primes up to 109 and stores them in the registers 3 to 29.
We don't check numbers that are divisible by 2 and 3.
Thus we start with 11, 13, 17, 19, 23, 25, …
This is the sequence of numbers congruent to 1 or 5 mod 6: A007310.
The difference of consecutive elements is: 2, 4, 2, 4, 2, 4, …
We can calculate that sequence recursively (lines 10-13 and 39-42):
\(\Delta n'=6-\Delta n\)
For numbers smaller than 112 = 121 we only have to check divisibility by 5 and 7.
In the second part (lines 36-60) we continue from 109 with 113, 115, 119, … and check only if they are divisible by primes smaller than their square root.
Due to the lack of both the x<y and x≥y instruction we have to use the combination of two comparison instructions:
Code:
49: 14 41 : x≤y ; n ≤ π(i)² ?
50: 14 71 : x=y ; n = π(i)² ?
51: 13 03 : GTO 3 ; n ≥ π(i)²The ISZ command has no upper limit. Thus we have to check that separately (lines 32-34).
In the second part I've omitted that check. The program will stop with an error for the first prime bigger than 1092 = 11,881 which is 11,887.
If you want a decent termination you can add a check after the PAUSE command.
Here's what I came up with:
Code:
01: 05 : 5
02: 23 00 : STO 0 ; i = 5
03: 23 03 : STO 3 ; π(3) = 5
04: 02 : 2
05: 23 01 : STO 1 ; ∆n = 2
06: 07 : 7
07: 23 02 : STO 2 ; n = 7
08: 23 04 : STO 4 ; π(4) = 7
09: 15 13 00 : LBL 0 ; for n loop
10: 06 : 6 ; 6
11: 24 01 : RCL 1 ; ∆n 6
12: 41 : - ; 6-∆n
13: 23 01 : STO 1 ; ∆n' = 6-∆n
14: 23 51 02 : STO+ 2 ; n' = n + ∆n
15: 24 02 : RCL 2 ; n
16: 05 : 5 ; 5 n
17: 71 : ÷ ; n/5
18: 15 62 : FRAC ; {n/5}
19: 15 71 : x=0 ; 5 | n ?
20: 13 00 : GTO 0 ; next n
21: 24 02 : RCL 2 ; n
22: 07 : 7 ; 7 n
23: 71 : ÷ ; n/7
24: 15 62 : FRAC ; {n/7}
25: 15 71 : x=0 ; 7 | n ?
26: 13 00 : GTO 0 ; next n
27: 24 02 : RCL 2 ; n
28: 23 22 : STO i ; π(i) = n
29: 14 74 : PAUSE ; is prime
30: 15 24 : ISZ ; i = i+1
31: 24 00 : RCL 0 ; i
32: 03 : 3
33: 00 : 0 ; 30 i
34: 14 51 : x>y ; 30 > i ?
35: 13 00 : GTO 0 ; next n
36: 15 13 01 : LBL 1 ; for n loop
37: 06 : 6 ; 6
38: 24 01 : RCL 1 ; ∆n 6
39: 41 : - ; 6-∆n
40: 23 01 : STO 1 ; ∆n' = 6-∆n
41: 23 51 02 : STO+ 2 ; n' = n + ∆n
42: 03 : 3
43: 23 00 : STO 0 ; i = 3
44: 15 13 02 : LBL 2 ; for i loop
45: 24 22 : RCL i ; π(i)
46: 24 02 : RCL 2 ; n π(i)
47: 24 22 : RCL i ; π(i) n π(i)
48: 71 : ÷ ; n/π(i) π(i)
49: 14 41 : x≤y ; n ≤ π(i)² ?
50: 14 71 : x=y ; n = π(i)² ?
51: 13 03 : GTO 3 ; n ≥ π(i)²
52: 24 02 : RCL 2 ; n
53: 14 74 : PAUSE ; is prime
54: 13 01 : GTO 1 ; next n
55: 15 13 03 : LBL 3 ; n/π(i)
56: 15 62 : FRAC ; {n/π(i)}
57: 15 71 : x=0 ; π(i) | n
58: 13 01 : GTO 1 ; next n
59: 15 24 : ISZ ; i = i+1
60: 13 02 : GTO 2 ; next iFor those like me that don't have an HP-29C I recommend this online emulator.
You can just copy the listing to the program that pops up when you click the display.
This is the listing of the registers:
Code:
00: i
01: ∆n
02: n
03: 5
04: 7
05: 11
06: 13
07: 17
08: 19
09: 23
10: 29
11: 31
12: 37
13: 41
14: 43
15: 47
16: 53
17: 59
18: 61
19: 67
20: 71
21: 73
22: 79
23: 83
24: 89
25: 97
26: 101
27: 103
28: 107
29: 109This list of the primes up to 11933 was copied from The First 10,000 Primes:
Code:
2 3 5 7 11 13 17 19 23 29
31 37 41 43 47 53 59 61 67 71
73 79 83 89 97 101 103 107 109 113
127 131 137 139 149 151 157 163 167 173
179 181 191 193 197 199 211 223 227 229
233 239 241 251 257 263 269 271 277 281
283 293 307 311 313 317 331 337 347 349
353 359 367 373 379 383 389 397 401 409
419 421 431 433 439 443 449 457 461 463
467 479 487 491 499 503 509 521 523 541
547 557 563 569 571 577 587 593 599 601
607 613 617 619 631 641 643 647 653 659
661 673 677 683 691 701 709 719 727 733
739 743 751 757 761 769 773 787 797 809
811 821 823 827 829 839 853 857 859 863
877 881 883 887 907 911 919 929 937 941
947 953 967 971 977 983 991 997 1009 1013
1019 1021 1031 1033 1039 1049 1051 1061 1063 1069
1087 1091 1093 1097 1103 1109 1117 1123 1129 1151
1153 1163 1171 1181 1187 1193 1201 1213 1217 1223
1229 1231 1237 1249 1259 1277 1279 1283 1289 1291
1297 1301 1303 1307 1319 1321 1327 1361 1367 1373
1381 1399 1409 1423 1427 1429 1433 1439 1447 1451
1453 1459 1471 1481 1483 1487 1489 1493 1499 1511
1523 1531 1543 1549 1553 1559 1567 1571 1579 1583
1597 1601 1607 1609 1613 1619 1621 1627 1637 1657
1663 1667 1669 1693 1697 1699 1709 1721 1723 1733
1741 1747 1753 1759 1777 1783 1787 1789 1801 1811
1823 1831 1847 1861 1867 1871 1873 1877 1879 1889
1901 1907 1913 1931 1933 1949 1951 1973 1979 1987
1993 1997 1999 2003 2011 2017 2027 2029 2039 2053
2063 2069 2081 2083 2087 2089 2099 2111 2113 2129
2131 2137 2141 2143 2153 2161 2179 2203 2207 2213
2221 2237 2239 2243 2251 2267 2269 2273 2281 2287
2293 2297 2309 2311 2333 2339 2341 2347 2351 2357
2371 2377 2381 2383 2389 2393 2399 2411 2417 2423
2437 2441 2447 2459 2467 2473 2477 2503 2521 2531
2539 2543 2549 2551 2557 2579 2591 2593 2609 2617
2621 2633 2647 2657 2659 2663 2671 2677 2683 2687
2689 2693 2699 2707 2711 2713 2719 2729 2731 2741
2749 2753 2767 2777 2789 2791 2797 2801 2803 2819
2833 2837 2843 2851 2857 2861 2879 2887 2897 2903
2909 2917 2927 2939 2953 2957 2963 2969 2971 2999
3001 3011 3019 3023 3037 3041 3049 3061 3067 3079
3083 3089 3109 3119 3121 3137 3163 3167 3169 3181
3187 3191 3203 3209 3217 3221 3229 3251 3253 3257
3259 3271 3299 3301 3307 3313 3319 3323 3329 3331
3343 3347 3359 3361 3371 3373 3389 3391 3407 3413
3433 3449 3457 3461 3463 3467 3469 3491 3499 3511
3517 3527 3529 3533 3539 3541 3547 3557 3559 3571
3581 3583 3593 3607 3613 3617 3623 3631 3637 3643
3659 3671 3673 3677 3691 3697 3701 3709 3719 3727
3733 3739 3761 3767 3769 3779 3793 3797 3803 3821
3823 3833 3847 3851 3853 3863 3877 3881 3889 3907
3911 3917 3919 3923 3929 3931 3943 3947 3967 3989
4001 4003 4007 4013 4019 4021 4027 4049 4051 4057
4073 4079 4091 4093 4099 4111 4127 4129 4133 4139
4153 4157 4159 4177 4201 4211 4217 4219 4229 4231
4241 4243 4253 4259 4261 4271 4273 4283 4289 4297
4327 4337 4339 4349 4357 4363 4373 4391 4397 4409
4421 4423 4441 4447 4451 4457 4463 4481 4483 4493
4507 4513 4517 4519 4523 4547 4549 4561 4567 4583
4591 4597 4603 4621 4637 4639 4643 4649 4651 4657
4663 4673 4679 4691 4703 4721 4723 4729 4733 4751
4759 4783 4787 4789 4793 4799 4801 4813 4817 4831
4861 4871 4877 4889 4903 4909 4919 4931 4933 4937
4943 4951 4957 4967 4969 4973 4987 4993 4999 5003
5009 5011 5021 5023 5039 5051 5059 5077 5081 5087
5099 5101 5107 5113 5119 5147 5153 5167 5171 5179
5189 5197 5209 5227 5231 5233 5237 5261 5273 5279
5281 5297 5303 5309 5323 5333 5347 5351 5381 5387
5393 5399 5407 5413 5417 5419 5431 5437 5441 5443
5449 5471 5477 5479 5483 5501 5503 5507 5519 5521
5527 5531 5557 5563 5569 5573 5581 5591 5623 5639
5641 5647 5651 5653 5657 5659 5669 5683 5689 5693
5701 5711 5717 5737 5741 5743 5749 5779 5783 5791
5801 5807 5813 5821 5827 5839 5843 5849 5851 5857
5861 5867 5869 5879 5881 5897 5903 5923 5927 5939
5953 5981 5987 6007 6011 6029 6037 6043 6047 6053
6067 6073 6079 6089 6091 6101 6113 6121 6131 6133
6143 6151 6163 6173 6197 6199 6203 6211 6217 6221
6229 6247 6257 6263 6269 6271 6277 6287 6299 6301
6311 6317 6323 6329 6337 6343 6353 6359 6361 6367
6373 6379 6389 6397 6421 6427 6449 6451 6469 6473
6481 6491 6521 6529 6547 6551 6553 6563 6569 6571
6577 6581 6599 6607 6619 6637 6653 6659 6661 6673
6679 6689 6691 6701 6703 6709 6719 6733 6737 6761
6763 6779 6781 6791 6793 6803 6823 6827 6829 6833
6841 6857 6863 6869 6871 6883 6899 6907 6911 6917
6947 6949 6959 6961 6967 6971 6977 6983 6991 6997
7001 7013 7019 7027 7039 7043 7057 7069 7079 7103
7109 7121 7127 7129 7151 7159 7177 7187 7193 7207
7211 7213 7219 7229 7237 7243 7247 7253 7283 7297
7307 7309 7321 7331 7333 7349 7351 7369 7393 7411
7417 7433 7451 7457 7459 7477 7481 7487 7489 7499
7507 7517 7523 7529 7537 7541 7547 7549 7559 7561
7573 7577 7583 7589 7591 7603 7607 7621 7639 7643
7649 7669 7673 7681 7687 7691 7699 7703 7717 7723
7727 7741 7753 7757 7759 7789 7793 7817 7823 7829
7841 7853 7867 7873 7877 7879 7883 7901 7907 7919
7927 7933 7937 7949 7951 7963 7993 8009 8011 8017
8039 8053 8059 8069 8081 8087 8089 8093 8101 8111
8117 8123 8147 8161 8167 8171 8179 8191 8209 8219
8221 8231 8233 8237 8243 8263 8269 8273 8287 8291
8293 8297 8311 8317 8329 8353 8363 8369 8377 8387
8389 8419 8423 8429 8431 8443 8447 8461 8467 8501
8513 8521 8527 8537 8539 8543 8563 8573 8581 8597
8599 8609 8623 8627 8629 8641 8647 8663 8669 8677
8681 8689 8693 8699 8707 8713 8719 8731 8737 8741
8747 8753 8761 8779 8783 8803 8807 8819 8821 8831
8837 8839 8849 8861 8863 8867 8887 8893 8923 8929
8933 8941 8951 8963 8969 8971 8999 9001 9007 9011
9013 9029 9041 9043 9049 9059 9067 9091 9103 9109
9127 9133 9137 9151 9157 9161 9173 9181 9187 9199
9203 9209 9221 9227 9239 9241 9257 9277 9281 9283
9293 9311 9319 9323 9337 9341 9343 9349 9371 9377
9391 9397 9403 9413 9419 9421 9431 9433 9437 9439
9461 9463 9467 9473 9479 9491 9497 9511 9521 9533
9539 9547 9551 9587 9601 9613 9619 9623 9629 9631
9643 9649 9661 9677 9679 9689 9697 9719 9721 9733
9739 9743 9749 9767 9769 9781 9787 9791 9803 9811
9817 9829 9833 9839 9851 9857 9859 9871 9883 9887
9901 9907 9923 9929 9931 9941 9949 9967 9973 10007
10009 10037 10039 10061 10067 10069 10079 10091 10093 10099
10103 10111 10133 10139 10141 10151 10159 10163 10169 10177
10181 10193 10211 10223 10243 10247 10253 10259 10267 10271
10273 10289 10301 10303 10313 10321 10331 10333 10337 10343
10357 10369 10391 10399 10427 10429 10433 10453 10457 10459
10463 10477 10487 10499 10501 10513 10529 10531 10559 10567
10589 10597 10601 10607 10613 10627 10631 10639 10651 10657
10663 10667 10687 10691 10709 10711 10723 10729 10733 10739
10753 10771 10781 10789 10799 10831 10837 10847 10853 10859
10861 10867 10883 10889 10891 10903 10909 10937 10939 10949
10957 10973 10979 10987 10993 11003 11027 11047 11057 11059
11069 11071 11083 11087 11093 11113 11117 11119 11131 11149
11159 11161 11171 11173 11177 11197 11213 11239 11243 11251
11257 11261 11273 11279 11287 11299 11311 11317 11321 11329
11351 11353 11369 11383 11393 11399 11411 11423 11437 11443
11447 11467 11471 11483 11489 11491 11497 11503 11519 11527
11549 11551 11579 11587 11593 11597 11617 11621 11633 11657
11677 11681 11689 11699 11701 11717 11719 11731 11743 11777
11779 11783 11789 11801 11807 11813 11821 11827 11831 11833
11839 11863 11867 11887 11897 11903 11909 11923 11927 11933The program should run as is on the HP-19C as well. You may want to replace the PAUSE statement with a PRx statement.
Edit: Added comments and stack-diagrams to the code and hopefully some clarifications to the description.
RE: (29C) Prime numbers up to 10'000 - PedroLeiva - 09-03-2018 01:14 PM
(09-02-2018 11:01 AM)Thomas Klemm Wrote: Inspired by Jürgen's video printing the Prime numbers up to 10'000 with an HP-19C I took a closer look at its sibling the HP-29C. It has 30 registers which is enough to store the primes up to 100. We can use indirect addressing with register 0 so I thought I give it a try.Very nice. It also works for HP 67, physical device and emulators.
For those like me that don't have an HP-29C I recommend this online emulator. You can just copy this listing to the program that pops up when you click the display.
..............................
The program should run as is on the HP-19C as well. You may want to replace the PAUSE statement with a PRx statement.
Some simply instructions:
To Start: RTN R/S
To Stop: R/S
To re-Start: R/S
About Emulators, thanks for the links to Sydneysmith site, I did not know it, it must be recent, for HP 29 and HP 67 the links are:
http://www.sydneysmith.com/products/hp29u/run/index.html
http://www.sydneysmith.com/products/gss-hp67u/run/
Great, TYVM. Pedro
RE: (29C) Prime numbers up to 10'000 - Dieter - 09-03-2018 05:24 PM
(09-03-2018 01:14 PM)PedroLeiva Wrote: Very nice. It also works for HP 67, physical device and emulators.
An HP67/97 version requires some changes, especially as the index register is no 0 but "I" which, unlike on the 19C/29C, is located at the end of the address range. This means that the test for the largest possible index has to be adjusted and you have to make sure that the last register "I" is not overwritten by a prime. Which is an important point as the original program is designed to terminate with an error message because of the attempt to access a register beyond 29. For the 67/97 storing something in R25 (="I") will not generate an error but it will mess up the program. On the other hand all indexes may be decreased by 1 as R0 is no longer required.
Since I am not completely sure how the program works maybe Thomas can provide some more detailled hints for an HP67/97 version. I have a version running but I am not sure if it terminates correctly.
(09-03-2018 01:14 PM)PedroLeiva Wrote: About Emulators, thanks for the links to Sydneysmith site, I did not know it, it must be recent, for HP 29 and HP 67 the links are: (...)
I like the Panamatik emulators, both for the 19C (with printer) and the 67.
Dieter
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 09-03-2018 05:46 PM
I wondered how much better using A007310 and primes instead of odd numbers would be.
The most expensive operation is probably checking if n is divisible by a probe π.
Thus I counted how often that happens.
n ∈ A007310 and π ∈ primes
This is essentially the algorithm of my program from above:
Code:
primes = [ 5, 7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101 ]
N = 10000
n = 5
Δn = 2
count = 0
while n < N:
for π in primes:
if n < π**2:
print(n)
break
count += 1
if n % π == 0:
break
n += Δn
Δn = 6 - Δn
print(count)count = 28760
both n and π are odd numbers
This is the algorithm used by Jürgen:
Code:
N = 10000
n = 5
Δn = 2
count = 0
while n < N:
for π in range(3, 102, Δn):
if n < π**2:
print(n)
break
count += 1
if n % π == 0:
break
n += Δn
print(count)count = 55958
That is more by a factor of 1.9456.
We can now gradually improve the algorithm and keep track of the changes.
n ∈ A007310 and π is odd
Code:
N = 10000
n = 5
Δn = 2
count = 0
while n < N:
for π in range(5, 102, Δn):
if n < π**2:
print(n)
break
count += 1
if n % π == 0:
break
n += Δn
Δn = 6 - Δn
print(count)count = 43625
n is odd and π ∈ {3} ∪ A007310
Code:
A007310 = [ 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35,
37, 41, 43, 47, 49, 53, 55, 59, 61, 65, 67,
71, 73, 77, 79, 83, 85, 89, 91, 95, 97, 101 ]
N = 10000
n = 5
Δn = 2
count = 0
while n < N:
for π in [ 3 ] + A007310:
if n < π**2:
print(n)
break
count += 1
if n % π == 0:
break
n += Δn
print(count)count = 40350
both n and π are ∈ A007310
Code:
A007310 = [ 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35,
37, 41, 43, 47, 49, 53, 55, 59, 61, 65, 67,
71, 73, 77, 79, 83, 85, 89, 91, 95, 97, 101 ]
N = 10000
n = 5
Δn = 2
count = 0
while n < N:
for π in A007310:
if n < π**2:
print(n)
break
count += 1
if n % π == 0:
break
n += Δn
Δn = 6 - Δn
print(count)count = 35354
This is better by a factor 1.5827 compared to the original program.
And that's close to 3 ÷ 2 = 1.5 what I would have expected.
However just using primes as probes is still better by a factor of 1.2292.
That's probably since in the worst case of n being a prime all probes have to be checked.
The length of primes is 27 and that of A007310 is 33.
This gives a factor of 1.2222.
Cheers
Thomas
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 09-03-2018 07:05 PM
(09-03-2018 05:24 PM)Dieter Wrote: Since I am not completely sure how the program works maybe Thomas can provide some more detailed hints for an HP67/97 version.
If we have only 26 registers we can only store these 23 primes:
[5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
But since register I occupies register 25 instead of 0 the whole list is shifted.
Thus the primes start with register 2 instead of 3 and so on.
We can't reuse some of the numbers anymore when initialising the registers but that shouldn't be a problem.
Code:
01: 2
02: STO 0 ; ∆n = 2
03: 7
04: STO 1 ; n = 7
05: STO 3 ; π(3) = 7
06: 5
07: STO 2 ; π(2) = 5
08: 4
09: STO I ; i = 4You have to adjust the registers through out the whole program.
In the for i loop we have to initialise register I with 2 now.
And then we better switch this with calculating the next n:
Code:
36: LBL 1 ; for n loop
37: 6 ; 6
38: RCL 0 ; ∆n 6
39: - ; 6-∆n
40: STO 0 ; ∆n' = 6-∆n
41: STO+ 1 ; n' = n + ∆n
42: 2
43: STO I ; i = 2
44: LBL 2 ; for i loopWe have to make sure that the index is smaller than 25.
Thus lines 32-33 have to be replaced by:
Code:
30: ISZ ; i = i+1
31: RCL I ; i
32: 2
33: 5 ; 25 i
34: x>y ; 25 > i ?
35: GTO 0 ; next nThe line numbering will probably be off as well.
Since 97² = 9409 we can't go further than that. But the paper will run out anyway at about 9’300 so that's not a problem.
If you want to make really sure that register I isn't overwritten you can add the same check at the end of the program:
Code:
59: ISZ ; i = i+1
60: GTO 2 ; next iYou may have noticed the Python programs in my last post. They might help understanding the algorithm.
HTH
Thomas
PS: When I contains 25 the command RCL (i) will return 25. Thus it will be checked if 25 is a divisor of n.
But we already checked if 5 is and that wasn't the case.
Furthermore the previous probe was 97 which is bigger and thus n > 25².
Thus we will happily increment I once more and get an error the next time we try to execute RCL (i).
At least that's what I assume would happen. I haven't tested it yet.
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 09-03-2018 09:26 PM
Listing for the HP-67:
Code:
001: 02: 2
002: 33 00: STO 0 ; ∆n = 2
003: 07: 7
004: 33 01: STO 1 ; n = 7
005: 33 03: STO 3 ; π(3) = 7
006: 05: 5
007: 33 02: STO 2 ; π(2) = 5
008: 04: 4
009: 35 33: ST I ; i = 4
010: 31 25 00: LBL 0 ; for n loop
011: 06: 6 ; 6
012: 34 00: RCL 0 ; ∆n 6
013: 51: - ; 6-∆n
014: 33 00: STO 0 ; ∆n' = 6-∆n
015: 33 61 01: STO+ 1 ; n' = n + ∆n
016: 34 01: RCL 1 ; n
017: 05: 5 ; 5 n
018: 81: ÷ ; n/5
019: 32 83: FRAC ; {n/5}
020: 31 51: x=0 ; 5 | n ?
021: 22 00: GTO 0 ; next n
022: 34 01: RCL 1 ; n
023: 07: 7 ; 7 n
024: 81: ÷ ; n/7
025: 32 83: FRAC ; {n/7}
026: 31 51: x=0 ; 7 | n ?
027: 22 00: GTO 0 ; next n
028: 34 01: RCL 1 ; n
029: 33 24: STO (i) ; π(i) = n
030: 35 72: PAUSE ; is prime
031: 31 34: ISZ ; i = i+1
032: 35 34: RC I ; i
033: 02: 2
034: 05: 5 ; 25 i
035: 32 81: x>y ; 25 > i ?
036: 22 00: GTO 0 ; next n
037: 31 25 01: LBL 1 ; for n loop
038: 06: 6 ; 6
039: 34 00: RCL 0 ; ∆n 6
040: 51: - ; 6-∆n
041: 33 00: STO 0 ; ∆n' = 6-∆n
042: 33 61 01: STO+ 1 ; n' = n + ∆n
043: 02: 2
044: 35 33: ST I ; i = 2
045: 31 25 02: LBL 2 ; for i loop
046: 34 24: RCL (i) ; π(i)
047: 34 01: RCL 1 ; n π(i)
048: 34 24: RCL (i) ; π(i) n π(i)
049: 81: ÷ ; n/π(i) π(i)
050: 32 71: x≤y ; n ≤ π(i)² ?
051: 32 51: x=y ; n = π(i)² ?
052: 22 03: GTO 3 ; n ≥ π(i)²
053: 34 01: RCL 1 ; n
054: 35 72: PAUSE ; is prime
055: 22 01: GTO 1 ; next n
056: 31 25 03: LBL 3 ; n/π(i)
057: 32 83: FRAC ; {n/π(i)}
058: 31 51: x=0 ; π(i) | n
059: 22 01: GTO 1 ; next n
060: 31 34: ISZ ; i = i+1
061: 22 02: GTO 2 ; next iRegisters:
Code:
00: ∆n
01: n
02: 5
03: 7
04: 11
05: 13
06: 17
07: 19
08: 23
09: 29
10: 31
11: 37
12: 41
13: 43
14: 47
15: 53
16: 59
17: 61
18: 67
19: 71
20: 73
21: 79
22: 83
23: 89
24: 97
25: iRE: (29C) Prime numbers up to 10'000 - Archilog - 09-04-2018 12:04 AM
Hello,
Great job!
Jürgen Keller's video actually was inspired from a program published in Byte a long time ago. A guy called C.Ret did an awesome work which can be found there: Silicium ... In French. But the most interesting is easily readable.
RE: (29C) Prime numbers up to 10'000 - Dieter - 09-04-2018 06:56 AM
(09-03-2018 09:26 PM)Thomas Klemm Wrote: Listing for the HP-67:
Code:
...
043: 03: 3
044: 35 33: ST I ; i = 3
Shouldn't this be 2 STO I ?
Otherwise my HP67-version is virtually identical.
OK, I added LBL A at the beginning. ;-)
Dieter
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 09-04-2018 09:31 AM
(09-04-2018 06:56 AM)Dieter Wrote: Shouldn't this be 2 STO I ?
Correct. I've updated the listing accordingly.
Thanks a lot for the notification.
Kind regards
Thomas
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 09-04-2018 06:53 PM
(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.
RE: (29C) Prime numbers up to 10'000 - pier4r - 09-04-2018 07:12 PM
Amazing work, little question. How much time does it take on the 29C (real hw) ?
RE: (29C) Prime numbers up to 10'000 - Albert Chan - 09-04-2018 07:31 PM
(09-04-2018 06:53 PM)Thomas Klemm Wrote: A guy called C. Ret did an awesome work ...
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.
I had a similar (*) code in Lua to build list of primes.
Instead of using min-heap, it uses Lua build-in hash table.
https://github.com/achan001/PrimePi/blob/master/prime.lua
(*) main difference is how sequence overlap is treated.
Above sequence for 7 and 11 overlap at 231, 385, 539, 693, ...
C. Ret code carried the prime along, ignoring the overlap.
Lua code skip over the overlap, with the sequence carried the right step.
Both approaches are equivalent.
RE: (29C) Prime numbers up to 10'000 - C.Ret - 09-05-2018 08:12 PM
(09-04-2018 06:53 PM)Thomas Klemm Wrote: This took me a while to figure out.
Here's how to initialize 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.
Good catch! You are right Thomas!
I have now to modify the explanations I gave in the Silicium forum.
Yes, register 5 is "hard coded" to be the root of the min-heap, and I mistakenly indicate to initiate it to 5. I was developing several versions of this code, having great troubles to make it fit into the 99 steps of program memory. A previous version (not published) increase R1 after storing the new square and prime value into register – so I mess up starting value needed in register R1.
Sorry for that and the wrong label for register R3 (I wrote 1 in the text when the correct register number is 3 as draw in the flow-chart). The first prime has to be put in register R3.
You may have spent a few time to figure out the bug! Please accept my apologies. If we meet, please remind me to offer you a large glass of beer (or any local beverage you prefer).
You also perfectly understand the philosophy of this algorithm. I try to reduce the number of tests needed to print the whole list of prims up to 10’000. Detecting prims by trying divisibility factors need a bunch of divisions. As you propose storing prime factor in memory is a great enhancement; the HP-19/29 has enough register to store all prime factors to achieve testing odd integer candidate up to 10’000.
In the same time, filtering candidate to test may reduce the effort. That what you made, generating a sequence {+4 +2 +4 +2 +4 +2 … } drastically reduce the number of composite candidates not prime. Unfortunately, detecting prime by successive factor divisibility is most efficient for composite (multiple) numbers; as soon as a divisor is found, we jump to next candidate. For prime, we must test all the possible factors along (i.e. up to the square root); no shortcut to next candidate.
Using a more sophisticated sequence, such as the one I have used ( i.e. {+4 +2 +4 +2 +4 +6 +2 +6 …} ) don’t help much; it reduce only by a few the “composite candidates” which are easy to detect and keep all the “long-testing true prime candidates” for divisibility tests.
This observation lead me to try to find another way for detecting prims. A way which may be as simple as a test. If the HP-19/29 had have more than 10’000 registers a sieve of Eratosthene may have been a easy way. So I look around to see if other sieve or tricky algorithms may help.
The sieve of Sundaram appears to be especially adapted as far as an adapted data structure can be bring out. As you correctly explain, we can't store in registers all the multiples tables produced by Sundaram’s algorithm. We only need to store the first value of all the sequences generated for each factor. We don’t need more register than for the divisibility testing method. But we need a convenient way to keep track of progress and to always know what the minimal value is. This can be a challenge and can ruin the effort.
A data structure can be of great advantage here, since it allows keeping track of minimum value with an optimal number of tests. The min-Heap structure is maintained so that his root (register 5) always contains the minimal composite candidate generated by any of the already encounted primes. Please note that this data-structure was only invented a few years before the HP-19C was built. Implementing a min-Heap in the HP-19/29 in 1972 may have been consider as really top new technology (Tip top Hype ! What was the words at this time ?)
By luck, the extra work to maintain the min-Heap will not ruin the effort because of two happy-facts:
First, we don’t need a full implementation of the min-Heap; since the Sundaram algorithm makes that any new sequence entry (for a new discovered prime factor) start as far as the square of this factor. So we don’t need “pop up” procedure, any new sequence entry can be securely add at the very end of the heap (pointed by value in register R1).
Second, the heap is short and not fully sorted. It is a binary structure, each node have two children’s which are greater (or egual). Less than 31 values mean that this binary tree is less than 4 levels high. When the tested candidate reach the min-heap root value, this value is increased by two times the factor. Then less than 4 values in the heap must be tested in order to replace (“push down”) the update root value.
In fact, only the top of the min-heap is active, middle range is quite passive (each level divide by 2 activity due to binary structure) and the bottom is a long time in stand-by until square values are reaches.
This made the code quite competitive despite the cost of extra test and the operations needed to maintain the heap in the case of a sophisticated filtered candidate. In contrary to the divisibility test methods, sieve methods didn’t lose the advantage of a few composite candidate number.
If we test all integers or all odd integers in a row, clearly all the over-costs for min-Heap operations make this algorithm slow. With sophisticated filtered candidate, the min-Heap method gains advantage since primality is detected with only one test. The way the Sundaram sequences are built spare operating time of the min-Heap.
But a lot of time is lost since the HP-19/29 have only one indirect register. Having two indirect registers may have simplify the operations for testing and swapping nodes in the binary heap and may have consequently greatly reduce running time.
Another difficulty is that the heap didn’t start at R0: or R1:, the off-set make the computation of node's and children's indices a bit Strange: j = 2*i or j = 2*i+1 is transformed as: j = (i-2)*2 or j = (i-2)*2+1
As Albert Chan notice it, my code ignore overlapping value. This is true and done this way expressly; overlapping value (i.e. same value in the sequence of different factors) intersecting the sophisticated filtered candidate sequence are quite rare and the expenses of code and supplementary tests needed to detect them didn’t enter the program memory. Moreover, due to the fact than the candidate are unequally spaced (+2 +4 or +6), for each candidate, more than one pre-compute multiple may have to be update. That why after the “min-Heap maintenance loop” the program returns to label 5 where actual candidate is tested again versus the new min-Heap root (aka Register 5) which may be an overlapping value or not. This make no difference. In theory, avoiding overlapping value will reduce heap stack operations number. In practice, even with more elaborated code no significant gain was observed.
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 09-07-2018 08:07 PM
(09-05-2018 08:12 PM)C.Ret Wrote: You may have spent a few time to figure out the bug!
Finding the bug was easy once I was sure that I entered the code correctly into the emulator.
But I've long been convinced that I made a typo somewhere.
And then there's a weird behaviour of the emulator that wouldn't allow me to single-step into a subroutine (e.g. GSB 8).
Instead it would just step over it and return.
I assume that this is a bug of the emulator.
But maybe someone can verify if that happens with a real calculator as well?
This made me insert an R/S command directly after the LBL 8.
Quote:Please accept my apologies. If we meet, please remind me to offer you a large glass of beer (or any local beverage you prefer).
Are you planning to attend the Allschwil Meeting 2018?
Because that would be an opportunity to meet.
I'll happily accept a beer, although that's not really necessary.
These kind of errors happen to me all the time and I'm glad I have Dieter who draws my attention to them, as in post #8 of this thread.
I'm very impressed by your flowchart and listing of the program.
What do you use to create them?
It appears that the listing can be executed.
Or how else should I interpret the smaller grey numbers like 1227x to the left of the line number?
Thanks a lot for your explanations!
Kind regards
Thomas
RE: (29C) Prime numbers up to 10'000 - rprosperi - 09-07-2018 09:16 PM
(09-07-2018 08:07 PM)Thomas Klemm Wrote: I'm very impressed by your flowchart and listing of the program.
What do you use to create them?
It appears that the listing can be executed.
Or how else should I interpret the smaller grey numbers like 1227x to the left of the line number?
Indeed the color-coded system of documenting the flowchart and synchronized source code is quite interesting, I've not seen anything like this before except in very expensive high-end professional software design tools, and that was a long time ago, so no idea if that caught-on.
How were these made, I too am curious about the numbers embedded in the listings? Please explain this when time allows.
Thanks for sharing this.
RE: (29C) Prime numbers up to 10'000 - Archilog - 09-09-2018 02:22 AM
(09-07-2018 08:07 PM)Thomas Klemm Wrote: I'm very impressed by your flowchart and listing of the program.
What do you use to create them?
It appears that the listing can be executed.
Or how else should I interpret the smaller grey numbers like 1227x to the left of the line number?
Those numbers represent the number of times the instruction is executed (did I say it is an awesome work?).
RE: (29C) Prime numbers up to 10'000 - Jurgen Keller - 11-11-2018 09:35 AM
It looks like I missed a lot of interesting stuff during my involuntary timeout. I'm glad that my video inspired other people to going deeper into this topic. It's always a great pleasure to read Thomas' (in HP terms) "well-engineered" articles. That's the kind of brainfood I need to keep a clear and fresh mind :-)
RE: (29C) Prime numbers up to 10'000 - Dieter - 11-11-2018 05:15 PM
(09-07-2018 08:07 PM)Thomas Klemm Wrote: And then there's a weird behaviour of the emulator that wouldn't allow me to single-step into a subroutine (e.g. GSB 8).
Instead it would just step over it and return.
I assume that this is a bug of the emulator.
No, that's authentic HP behaviour. ;-)
Indeed there are calculators where SST at a subroutine call executes that routine as a whole. This may have advantages sometimes, on other occasions it requires an immediate [R/S] to continue step by step. The HP67 and 97 behave the same way.
It looks like this is one of the many cases where HP deliberately changed the behaviour of their calculators – consistency obviously what not the most important design goal.
Dieter
RE: (29C) Prime numbers up to 10'000 - Dieter - 11-11-2018 05:21 PM
(11-11-2018 09:35 AM)Jurgen Keller Wrote: It looks like I missed a lot of interesting stuff during my involuntary timeout.
Welcome back, Jürgen. ;-)
Dieter
RE: (29C) Prime numbers up to 10'000 - Thomas Klemm - 12-20-2018 06:08 AM
Meanwhile C.Ret adapted the program for the HP-41C:
Here's the same program for the HP-42S:
Code:
00 { 188-Byte Prgm }
01▸LBL "VPRIMS"
02 CLRG
03 FIX 00
04 CF 29
05 STO 25
06 7
07 STO 24
08 "2 3 5"
09 3
10 STO 23
11 XEQ 08
12▸LBL 01
13 XEQ 02
14 XEQ 02
15 4
16 XEQ 04
17 6
18 XEQ 04
19 XEQ 03
20 6
21 XEQ 04
22 GTO 01
23▸LBL 02
24 4
25 XEQ 04
26▸LBL 03
27 2
28▸LBL 04
29 STO+ 24
30▸LBL 05
31 RCL 25
32 RCL 24
33 RCL 01
34 IP
35 X=Y?
36 RTN
37 X>Y?
38 GTO 08
39 LASTX
40 FP
41 200
42 ×
43 STO+ 01
44 1
45 STO ST L
46 RCL 01
47▸LBL 06
48 RCL ST L
49 STO+ ST L
50 X<>Y
51 RCL 00
52 LASTX
53 -
54 X<0?
55 GTO 05
56 RCL IND ST L
57 X<>Y
58 X≠0?
59 ISG ST L
60 GTO 07
61 R↓
62 RCL IND ST L
63 X>Y?
64 DSE ST L
65 X≤Y?
66 X<>Y
67▸LBL 07
68 R↓
69 X>Y?
70 GTO 05
71 STO IND ST Z
72 X<>Y
73 STO IND ST L
74 GTO 06
75▸LBL 08
76 R↓
77 X≤Y?
78 ISG 23
79 GTO 00
80 RCL 23
81 5
82 MOD
83 X≠0?
84 GTO 09
85 PRA
86 CLA
87▸LBL 09
88 X≠0?
89 ├" "
90 ARCL ST Y
91 R↓
92 X↑2
93 X≤Y?
94 ISG 00
95 RTN
96 LASTX
97 1/X
98 %
99 +
100 STO IND 00
101 RTN
102▸LBL 00
103 PRA
104 VIEW 23
105 ADV
106 FIX 04
107 SF 29
108 STOP
109 ENDAnd this is the print-out after running:
1E4 XEQ "VPRIMS"
Cheers
Thomas
PS: Instead of ARCL here:
Code:
90 ARCL ST Y
91 R↓We could use AIP:
Code:
88 R↓
89 AIPAnd thus get rid of these lines:
Code:
03 FIX 00
04 CF 29
…
106 FIX 04
107 SF 29