(29C) Prime numbers up to 10'000

[Thomas Klemm]

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 11² = 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 109² = 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 i

For 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: 109

This 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  11933

The 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.

[PedroLeiva]

(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.

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.

Very nice. It also works for HP 67, physical device and emulators.
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/hp29...index.html
http://www.sydneysmith.com/products/gss-hp67u/run/

Great, TYVM. Pedro

[Dieter]

(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

[Thomas Klemm]

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

[Thomas Klemm]

(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 = 4

You 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 loop

We 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 n

The 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 i

You 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.

[Thomas Klemm]

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 i

Registers:

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: i

[Archilog]

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.

[Dieter]

(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

[Thomas Klemm]

(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

[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.

[pier4r]

Amazing work, little question. How much time does it take on the 29C (real hw) ?

[Albert Chan]

(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.../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.

[C.Ret]

(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.

[Thomas Klemm]

(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

[rprosperi]

(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.

[Archilog]

(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?).

[Jurgen Keller]

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 :-)

[Dieter]

(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

[Dieter]

(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

[Thomas Klemm]

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 END

And 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 AIP

And thus get rid of these lines:

Code:

03 FIX 00
04 CF 29
…
106 FIX 04
107 SF 29