[HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) (/thread-12406.html) Pages: 1 2 3 [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) - fred_76 - 02-11-2019 05:16 PM Hello My name is Fred, I am from France. I am used to HP calculators since I was young in the 80's. My father was working in a military-like industry and was supplied with new HP calc every two or three years. Therefore I was always fed with the previous calc he had. It all started to me with the grey HP65, then the HP67, the HP34, the HP41 (don't remember which one exactly). I don't have them anymore, they were trashed ... sigh ... I bought a HP32S and that was one of the very best I had. It was stollen during a company convention... I then asked the company to buy me a calc, it was a HP48S, again it was stollen when my office was broken ! The company bought me another calc, a HP48GII that I found desesperatly anoying and hard to program. My sons just offered me a HP35S and I rediscovered the pleasure of the RPN calculators. Here is a code that I find quite effective to determine if a number is prime or not, and that also calculates the factors. It make an extensive use of the stack to avoid unecessary STO/RCLs. It also makes use of fast tricks, such as storing repetitive constants and using included operations. For example [6] [+] is about 4.3 times slower than [RCL+S] (with 6 previously stored in S). The logic below is a brut force with elimination of the trivial factors that are multiples of 2, 3 and 5. From 7, you just loop by adding 4,2,4,2,4,6,2,6 to the trial factor. Here is the code : Code: ``` P001    LBL P    :Prime calc routine P002    ENTER    :stack N in the pile P003    ENTER    :to Z  P004    ENTER    :so that it is in T at the next line          P005    2     P006    STO T    :store 2 (t as two) P007    RCL+T    :add 2 P008    STO F    :store 4 (f as four) P009    RCL+T    :add 2 P010    STO S    :store 6 (S as six)          P011    RCL-F    :subtract 4 to get 2 P013    RMDR    :check if N is a mult of 2 P013    x=0?            :if rest is =0 then N is a mult of 2 P014    GTO P087    :so goto end of calc          P015    CLx     P016    3            :check if N is a mult of 3 P017    RMDR     P018    x=0?     P019    GTO P087              P020    CLx     P021    LAST x     P022    RCL+T    :check if N is a mult of 5 P023    RMDR     P024    x=0?     P025    GTO P087              P026    CLx     P027    LAST x     P028    RCL+T    :check if N is a mult of 7 P029    RMDR     P030    x=0?     P031    GTO P087                              :from that point, it will loop and add successively 4,2,4,2,4,6,2,6 to the                  :factor and check is N is a mult of the factor          P032    CLx     P033    LAST X     P034    RCL+F    :4 P035    RMDR     P036    x=0?     P037    GTO P087              P038    CLx     P039    LAST X     P040    RCL+T    :2 P041    RMDR     P042    x=0?     P043    GTO P087              P044    CLx     P045    LAST X     P046    RCL+F    :4 P047    RMDR     P048    x=0?     P049    GTO P087              P050    CLx     P051    LAST X     P052    RCL+T    :2 P053    RMDR     P054    x=0?     P055    GTO P087              P056    CLx     P057    LAST X     P058    RCL+F    :4 P059    RMDR     P060    x=0?     P061    GTO P087              P062    CLx     P063    LAST X     P064    RCL+S    :6 P065    RMDR     P066    x=0?     P067    GTO P087              P068    CLx     P069    LAST X     P070    RCL+T    :2 P071    RMDR     P072    x=0?     P073    GTO P087              P074    CLx     P075    LAST X     P076    RCL+S    :6 P077    RMDR     P078    x=0?     P079    GTO P087                                :check if the factor^2 is greater than N,                  :(SQR is 30% faster than SQRT)                 :if greater, then N is prime, otherwise loop P080    CLx     P081    LAST X     P082    x² P083    x n, you can pre-calculate number of prime-wheel turns needed. Example, N = 99999 999977 ceil( √N ) = 316228 Max wheel turns = ceil((316228-7)/30) = 10541, simple loop count to 0 suffice. RE: Hello and program for prime number (HP35s) - fred_76 - 02-13-2019 10:09 AM Albert Good try. I tried, replacing the lines : Code: ```CLx LAST X x² x 0 OR 1 as A is composite or probably prime. RABIN? 571 A, B -> 0 or 1, 1 if A is a strong pseudo-prime base B, 0 if not. Though I haven't tried but if PRIME? uses RABIN? (which I assume implements the Miller-Rabin Primality Test) it might be considerably faster for bigger numbers. I've implemented this for the HP-48GX: (48G) Miller-Rabin Primality Test This article gives a bit more details: Miller-Rabin Primality Test for the HP-48 Cheers Thomas RE: [HP35s] Hello and program for prime number - fred_76 - 02-13-2019 04:03 PM Thank you Thomas for this idea. But I'm definitely not friendly enough with HP48GX language. In fact I was never able to understand it, when I had this calculator, it wasn't sold with any programming manual, and it is hard to understand it without any explanation... It may however be possible to implement that on HP35S, who knows ! RE: [HP35s] Hello and program for prime number - fred_76 - 02-13-2019 04:44 PM Thomas, In the library, I have a noob question : In line 92, what is "A-1" ? How do you get this ? Quote:(...) 90 RTN 91 CLx 92 A-1 93 x<>y 94 XEQ Z008 (...) I can't find it on the manual... RE: [HP35s] Hello and program for prime number - Albert Chan - 02-13-2019 05:26 PM Just to show what is involved doing Primality SPRP Test. Example, Prove X = 56789 is composite X-1 = 56788 = 2^2 * 14197 = 2^s * d → s = 2 → d = 14197 = 0b11,011,101,110,101 Check if X is an a-SPRP, for a=2: Bits Mod X 1 ﻿ ﻿ ﻿ ﻿ 2^1 ≡ 2 1 ﻿ ﻿ ﻿ ﻿ 2^3 ≡ 2² * 2 ≡ 8 0 ﻿ ﻿ ﻿ ﻿ 2^6 ≡ 8² ≡ 64 1 ﻿ ﻿ ﻿ ﻿ 2^13 ≡ 64² * 2 ≡ 8192 1 ﻿ ﻿ ﻿ ﻿ 2^27 ≡ 8192² * 2 ≡ 25321 1 ﻿ ﻿ ﻿ ﻿ 2^55 ≡ 25321² * 2 ≡ 10462 0 ﻿ ﻿ ﻿ ﻿ 2^110 ≡ 10462² ≡ 21041 1 ﻿ ﻿ ﻿ ﻿ 2^221 ≡ 21041² * 2 ≡ 50063 1 ﻿ ﻿ ﻿ ﻿ 2^443 ≡ 50063² * 2 ≡ 13275 1 ﻿ ﻿ ﻿ ﻿ 2^887 ≡ 13275² * 2 ≡ 18716 0 ﻿ ﻿ ﻿ ﻿ 2^1774 ≡ 18716² ≡ 14104 1 ﻿ ﻿ ﻿ ﻿ 2^3549 ≡ 14104² * 2 ≡ 38687 0 ﻿ ﻿ ﻿ ﻿ 2^7098 ≡ 38687² ≡ 9874 1 ﻿ ﻿ ﻿ ﻿ 2^14197 ≡ 9874² * 2 ≡ 35115 → a^d ≠ ±1 0 ﻿ ﻿ ﻿ ﻿ 2^28394 ≡ 35115² ≡ 3668 → a^(2^r*d) ≠ -1, r < s -> 56789 is *not* 2-SPRP -> 56789 is *guaranteed* composite 2-SPRP may still be composite, but much less likely. Furthur testing may required. An issue is how to do the squaring and multiply without overflow. It may require splitting the numbers, and do in steps. RE: [HP35s] Hello and program for prime number - Thomas Klemm - 02-13-2019 07:20 PM (02-13-2019 04:44 PM)fred_76 Wrote:  In line 92, what is "A-1" ? How do you get this ? Not sure if I remember that correctly but this looks like an equation. So you would enter it: [EQN] [RCL] A - 1 [C] And since that's probably the next question, you can use [STO] to store intermediate results in a variable which is indicated by ►: Code: `182    C+1►C` This increments variable C. But I haven't used my HP-35s for some time and right now the batteries are dead so I can't verify. In the manual check chapter 6: "Entering and Evaluating Equations". HTH Thomas RE: [HP35s] Hello and program for prime number - fred_76 - 02-14-2019 01:10 PM (thank you Thomas, it works !) Here are the comparison of the versions of the program I tested so far. The version that is presented at the top of this thread is the version 3, i.e. the optimized one for "brute force". Next step would be to implement the AKS, Miller-Rabin or SPRP, but those algorithms are quite difficult to program to my level... Version 1 The principle was to use only the stack registers, i.e. no variable. The trial factors steps are such that the multiples of 2, 3 and 5 are discarded from the tests. These steps are 4,2,4,2,4,6,2,6 starting from 7. The main loop is : Code: ```CLx (trial factor step) LAST X + RMDR x=0? GTO :is_not_prime``` At the end of the cycle, there is a test to check if the trial factor is lower than the square root of the number to test. In fact, as the calculation of x² is faster than the calculation of √x, the comparison is made between the square of the trial factor and the number. Code: ```CLx LAST X x² x>=y? GTO :is_prime GTO :loop_again``` Run time 6 digits 999 983 0'18.42'' Version 2 The HP35S takes far more time to evaluate the constants than recalling them from a variable. For example the following code : Code: ```1 +``` takes about 2.2 more time than (1 is stored in variable I) : Code: ```RCL I +``` The 2nd version code is : Code: ```CLx RCL (trial_factor_step) LAST X + RMDR x=0? GTO :is_not_prime``` Run time 6 digits 999 983 0'12.16'' 1.5 faster than original version Version 3 Another optimization is done. The HP35S can make some operations together with recalling/storing variables. And this is quite effective. If 1 is stored in I, the following code Code: ```RCL I +``` takes about 1.9 more times than Code: `RCL+ I` . The 3rd version code is : Code: ```CLx LAST X RCL+ (trial_factor_step) RMDR x=0? GTO :is_not_prime``` A second optimization is done on the exit part. When a test is done, it runs faster when the test returns "true" than when it returns "false", and also when the condition is exclusive (< or >) than inclusive (<= or >=). Code: ```CLx LAST X x² x= n: break     r = n % (w + 2)          # last spoke, reset r to avoid overflow return n                     # no factors``` RE: [HP35s] Program for prime number (brut force) - Dave Britten - 02-14-2019 06:48 PM (02-12-2019 12:52 AM)Don Shepherd Wrote:  Fred, here is a prime factor finder I wrote for the 12c+, based on work done by Dave Britten. Like yours, it also excludes multiples of 2, 3, and 5 from the trial factor pool. It determines 999,863 is prime in 2 seconds. 9,999,991 in 5 seconds 99,999,989 in 16 seconds 999,999,937 in 50 seconds It won't work with your two larger examples. And my work is of course just stolen from the HP-67 Math Pac 2. http://www.hpmuseum.org/software/67pacs/67factor.htm It's a super simple (and useful) algorithm to port to pretty much any programmable calculator I happen to be using. RE: [HP35s] Program for prime number (brut force) - fred_76 - 02-14-2019 07:51 PM (02-14-2019 06:26 PM)Albert Chan Wrote:  What happens if you do not check each remainder goes to zero, but collect it all per turns ? If product of the turn remainders reached zero, it meant you hit a factor there. You do not even have to check 2,3,5,7 ! Like this: Code: ```r = (n % 2) * (n % 3) * (n % 5) for turn in range(1, n):     # actual termination by size of factors     # NOTE: wheel had only 7 spokes = (7,11,13,17,19,23,29) + 30*(turn-1)     for w in primewheel(turn): r *= n % w     if r == 0: return turn   # factor found, factor < 30*turn     if w * w >= n: break     r = n % (w + 2)          # last spoke, reset r to avoid overflow return n                     # no factors``` You are right, this would simplify the code. But unfortunately, at least on the HP35s, the instruction x=0? (about 3 ms) is 4x faster than STO*Var (about 12 ms). Thus you will end-up with a run time increased by a factor of about 4. RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-14-2019 08:50 PM (02-14-2019 06:26 PM)Albert Chan Wrote:  What happens if you do not check each remainder goes to zero, but collect it all per turns ? Here's a program for the HP-42S using this idea: Code: ```00 { 120-Byte Prgm } 01▸LBL "PRIME?" 02 STO 09 03 1 04 STO 00 05 7 06 STO 01 07 11 08 STO 02 09 13 10 STO 03 11 17 12 STO 04 13 19 14 STO 05 15 23 16 STO 06 17 29 18 STO 07 19 30 20 STO 08 21 RCL 09 22 RCL 09 23 RCL 09 24 2 25 MOD 26 X<>Y 27 3 28 MOD 29 × 30 X<>Y 31 5 32 MOD 33 × 34▸LBL 00 35 X<>Y 36 RCL 01 37 MOD 38 × 39 X<>Y 40 RCL 02 41 MOD 42 × 43 X<>Y 44 RCL 03 45 MOD 46 × 47 X<>Y 48 RCL 04 49 MOD 50 × 51 X<>Y 52 RCL 05 53 MOD 54 × 55 X<>Y 56 RCL 06 57 MOD 58 × 59 X<>Y 60 RCL 07 61 MOD 62 × 63 X=0? 64 RTN 65 CLX 66 RCL 07 67 X↑2 68 X>Y? 69 RTN 70 CLX 71 RCL 08 72 STO+ 00 73 STO+ 01 74 STO+ 02 75 STO+ 03 76 STO+ 04 77 STO+ 05 78 STO+ 06 79 STO+ 07 80 CLX 81 RCL 00 82 MOD 83 GTO 00 84 END``` Translation should be straight forward for the HP-35s or pretty much any calculator that supports MOD and register arithmetics. It works only for input > 29 and returns 0 for composite numbers. In this case one of the numbers in registers 00-07 is a factor. Otherwise it returns a number a bit bigger than the input. Examples: 24883 XEQ "PRIME?" y: 24883 x: 0 24883 = 149 × 167 99999999977 XEQ "PRIME?" y: 99999999977 x: 100000780441 99999999977 is a prime. Cheers Thomas RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-14-2019 09:17 PM (02-14-2019 07:51 PM)fred_76 Wrote:  You are right, this would simplify the code. But unfortunately, at least on the HP35s, the instruction x=0? (about 3 ms) is 4x faster than STO*Var (about 12 ms). Thus you will end-up with a run time increased by a factor of about 4. It also depends on how many steps are executed within the loop. Right now it's 49 vs. 53, so we gain a few steps using Albert's suggestion. You can also try numbered registers instead of variables. I could imagine that they are faster. But I could be completely wrong. Best regards Thomas RE: [HP35s] Program for prime number (brut force) - fred_76 - 02-15-2019 08:47 AM (02-14-2019 08:50 PM)Thomas Klemm Wrote:  Here's a program for the HP-42S using this idea: Code: ```00 { 120-Byte Prgm } 01▸LBL "PRIME?" 02 STO 09 03 1 04 STO 00 05 7 06 STO 01 07 11 08 STO 02 09 13 10 STO 03 11 17 12 STO 04 13 19 14 STO 05 15 23 16 STO 06 17 29 18 STO 07 19 30 20 STO 08 21 RCL 09 22 RCL 09 23 RCL 09 24 2 25 MOD 26 X<>Y 27 3 28 MOD 29 × 30 X<>Y 31 5 32 MOD 33 × 34▸LBL 00 35 X<>Y 36 RCL 01 37 MOD 38 × 39 X<>Y 40 RCL 02 41 MOD 42 × 43 X<>Y 44 RCL 03 45 MOD 46 × 47 X<>Y 48 RCL 04 49 MOD 50 × 51 X<>Y 52 RCL 05 53 MOD 54 × 55 X<>Y 56 RCL 06 57 MOD 58 × 59 X<>Y 60 RCL 07 61 MOD 62 × 63 X=0? 64 RTN 65 CLX 66 RCL 07 67 X↑2 68 X>Y? 69 RTN 70 CLX 71 RCL 08 72 STO+ 00 73 STO+ 01 74 STO+ 02 75 STO+ 03 76 STO+ 04 77 STO+ 05 78 STO+ 06 79 STO+ 07 80 CLX 81 RCL 00 82 MOD 83 GTO 00 84 END``` Translation should be straight forward for the HP-35s or pretty much any calculator that supports MOD and register arithmetics. (...) Cheers Thomas Version 8 I entered the program in the HP35s, it was straightforward. The measured run times are still slower than V3 : Run time 6 digits 999 983 = 13.4 s 7 digits 9 999 991 = 40,4 s The difference between the two algorithms at each cycle are : V3 : stack ops = 16 | additions = 8 | multiplications = 0 | MOD = 8 | tests = 9 V8 : stack ops = 9 | additions = 8 | multiplications = 7 | MOD = 8 | tests = 2 The 7 multiplications in V8 take longer to calculate than the additional 7 tests and stack operations in V3. (02-14-2019 09:17 PM)Thomas Klemm Wrote:  You can also try numbered registers instead of variables. I could imagine that they are faster. But I could be completely wrong. There is no "numbered" registers in the HP35S to my knowledge, just alpha registers and indirect registers. I measured the time taken (1000 loops) to calculate the execution time of indirect vs direct registers. The unit op exec time is then : RCL (I) = 10 ms RCL I = 5 ms We should better use direct registers to win time. RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-15-2019 03:07 PM (02-14-2019 01:10 PM)fred_76 Wrote:  Version 3 Run time 6 digits 999 983 0m 9.8s (about 2x faster than V0) 7 digits 9 999 991 0m 28s 8 digits 99 999 989 1m 28s 9 digits 999 999 937 4m 36s 10 digits 9 999 999 967 14m 37s 11 digits 99 999 999 977 46m 26s 12 digits 999 999 999 989 2h 32m 18s Curious about performance of checks vs. size of N Redo above tables, with cases checked = 4 + 8 * ceil((√(N) - 7)/30) Code: ```Digits Time(sec) Checks  Speed(msec per check) 6      9.8       276     35.5 7      28        852     32.9 8      88        2676    32.9 9      276       8434    32.7 10     877       26676   32.9 11     2786      84332   33.0 12     9138      266676  34.3``` It seems check speed unaffected by size of N. Timings depends mainly on number of checks required. It is possible to reduce required checks, by few applications of Fermat's method. Example, N = 999,999,999,989. Tried to complete the squares: Code: ```X          Y^2 = X^2-N 1e6        11 1e6 + 1    2e6 + 1 + 11 1e6 + 2    4e6 + 3 + 12 1e6 + 3    6e6 + 5 + 15 1e6 + 4    8e6 + 7 + 20 1e6 + 5    10e6 + 9 + 27 1e6 + 6    12e6 + 11 + 36 1e6 + 7    14e6 + 13 + 47 1e6 + 8    16e6 + 15 + 60 1e6 + 9    18e6 + 17 + 75``` If Y is integer, above can be factored as (X + Y)(X - Y) None of above Y's worked, but it reduce cases checked greatly. (BTW, only the Y^2 ends in 36 need checking, others will never be integer square) Max factors to search = floor(next(X) - next(Y)) = floor((1e6 + 10) - √(20e6 + 19 + 92)) = 995537 -> Cases to check = 4 + 8 * ceil((995537 - 7) / 30) = 4 + 8 * 33185 = 265484 -> Removed checks = 266676 - 265484 = 1192 (or, 1192/8 = 149 turns)