[HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
02-11-2019, 05:16 PM (This post was last modified: 04-12-2019 12:44 PM by fred_76.)
Post: #1
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
[HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
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<y?     P084    GTO P032              P085    Rv            :if N is prime P086    RTN            :stop and show X=Y=N          P087    CLx            :if N is not prime P088    LAST X    :show first factor P089 x=y?            :for 2-3-5-7 P090 GTO P085    :2-3-5-7 are prime  P091    STOP    :stop to show result X=factor, Y=N P092    INT/            :press R/S to P093    GTO P001    :start again to find the next factor P094    RTN

To use it, just enter the number to test, then press XEQ P.
• If the number is prime, X=Y=Number
• If the number is not prime, X=first factor found, Y=Number
• then press R/S to continue the calc and find the next factors until X=Y

It runs quite fast on my HP35S :

PHP Code:
Digits        Prime Number           Total Time6 digits       999 863               00:00:107 digits       9 999 991             00:00:288 digits       99 999 989            00:01:289 digits       999 999 937           00:04:3610 digits      9 999 999 967         00:14:3711 digits      99 999 999 977        00:46:2612 digits      999 999 999 989       02:32:18

All the best

Fred

--- edit : added time for 12 digits
--- edit : added lines P089/P090 to handle exceptions 2/3/5/7
02-12-2019, 12:52 AM (This post was last modified: 02-12-2019 01:09 AM by Don Shepherd.)
Post: #2
 Don Shepherd Senior Member Posts: 713 Joined: Dec 2013
RE: Hello and program for prime number (HP35s)
(02-11-2019 05:16 PM)fred_76 Wrote:  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 whiche one exactly). I don't have them anymore, they were trashed ... sight ...

I bought a HP32S and that was one on the very best I had. It was stollen...

I then asked the company to buy me a calc, it was a HP48S, again it was stollen ! 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 small 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 to chack all potential factors, eliminating 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 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<y?     P084    GTO P032              P085    Rv            :if N is prime P086    RTN            :stop and show X=Y=N          P087    CLx            :if N is not prime P088    LAST X    :show first factor P089    STOP    :stop to show result X=factor, Y=N P090    INT/            :press R/S to P091    GTO P001    :start again to find the next factor P092    RTN

To use it, just enter the number to test, then press XEQ P.
• If the number is prime, X=Y=Number
• If the number is not prime, X=first factor found, Y=Number
• then press R/S to continue the calc and find the next factors until X=Y

It runs quite fast on my HP35S :

PHP Code:
Digits        Prime Number           Total Time6 digits       999 863               00:00:09.537 digits       9 999 991             00:00:28.198 digits       99 999 989            00:01:27.709 digits       999 999 937           00:04:36.3510 digits      9 999 999 967         00:14:36.5911 digits      99 999 999 977        00:46:26.18

All the best

Fred

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.
02-12-2019, 09:28 AM
Post: #3
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: Hello and program for prime number (HP35s)
Hi Don,

Wow, the 12c+ is quite fast. I shall find a way to optimize the code of the 35s a lot to go as fast as the 12c+ !
02-12-2019, 01:54 PM
Post: #4
 Don Shepherd Senior Member Posts: 713 Joined: Dec 2013
RE: Hello and program for prime number (HP35s)
(02-12-2019 09:28 AM)fred_76 Wrote:  Hi Don,

Wow, the 12c+ is quite fast. I shall find a way to optimize the code of the 35s a lot to go as fast as the 12c+ !

The quickness of the 12c+ is due to the very fast processor in this case, not the optimization of RPN code (although I always like to optimize that as much as possible). Running the same code on an older 12c is much, much, much slower.
02-12-2019, 06:27 PM
Post: #5
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: Hello and program for prime number (HP35s)
I just found that archived discussion :

The optimization I propose makes the code running about twice as fast. I can’t get it run faster... (I could get about 3.5% reduction by excluding the multiples of 7 but it would add about 250 lines of code. Not worth the effort).
02-12-2019, 08:14 PM (This post was last modified: 02-15-2019 01:02 PM by Albert Chan.)
Post: #6
 Albert Chan Senior Member Posts: 943 Joined: Jul 2018
RE: Hello and program for prime number (HP35s)
Hi, fred_76

Termination code may be optimized a bit.

Instead of testing factor^2 > 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.
02-13-2019, 10:09 AM
Post: #7
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: Hello and program for prime number (HP35s)
Albert

Good try. I tried, replacing the lines :

Code:
CLx LAST X x² x<y? GTO :loop

by (after initialisation of I=(N^0.5-7)%30+1) :

Code:
 DSE I GTO :loop

But it is slower...

In fact, the first code takes about 25 ms to execute while the second code takes about 35 ms. The HP35s is fast when it calculates x², but slow to execute DSE.

For 999 999 937, the DSE modified code takes 4'57" compared to 4'36".
02-13-2019, 11:53 AM
Post: #8
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: [HP35s] Hello and program for prime number
You may want to have a look at (35S) Number Theory Library.
There are the following functions:

Quote:PRIME? 148
A -> 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

Miller-Rabin Primality Test for the HP-48

Cheers
Thomas
02-13-2019, 04:03 PM
Post: #9
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Hello and program for prime number
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 !
02-13-2019, 04:44 PM
Post: #10
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Hello and program for prime number
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...
02-13-2019, 05:26 PM (This post was last modified: 02-13-2019 07:23 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 943 Joined: Jul 2018
RE: [HP35s] Hello and program for prime number
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.
02-13-2019, 07:20 PM
Post: #12
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: [HP35s] Hello and program for prime number
(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
02-14-2019, 01:10 PM (This post was last modified: 02-14-2019 04:49 PM by fred_76.)
Post: #13
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Hello and program for prime number
(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<y? GTO :loop_again GTO :is_prime

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

Version 4
Instead of calculating x² at each end of cycle, store √N in a variable and compare with the trial factor. The problem is that the stack is quite tricky to restore and those operations take slightly more time than the calculation of x², therefore it is not a good idea.

Run time
8 digits 99 999 989 1m 30s (2 s more than V3)

Version 5
Based on V3, after idea of Albert Chan, precalculate the max number of cycles to be made C = INT[(N-7)/30]+1 and use of DSE. The problem is that DSE takes more time than calculating x².

Run time
8 digits 99 999 989 1m 34s (6 s more than V3)

Version 6
Instead of duplicating the same lines of code for each trial factor step, call a subroutine.
But calling subroutines takes more time. It is good for the sake of the HP35s memory, but not time effective.

Run time
8 digits 99 999 989 1m 42s (14 s more than V3)

Version 7
Just for fun, store trial factor steps in indirect registers and call subroutines. This is definitely not time effective.

Run time
8 digits 99 999 989 3m 12s (2.2x slower than V3)
02-14-2019, 06:26 PM
Post: #14
 Albert Chan Senior Member Posts: 943 Joined: Jul 2018
RE: [HP35s] Program for prime number (brut force)
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
02-14-2019, 06:48 PM
Post: #15
 Dave Britten Senior Member Posts: 1,492 Joined: Dec 2013
RE: [HP35s] Program for prime number (brut force)
(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.
02-14-2019, 07:51 PM
Post: #16
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Program for prime number (brut force)
(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.
02-14-2019, 08:50 PM
Post: #17
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: [HP35s] Program for prime number (brut force)
(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
02-14-2019, 09:17 PM
Post: #18
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: [HP35s] Program for prime number (brut force)
(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
02-15-2019, 08:47 AM (This post was last modified: 02-15-2019 09:06 AM by fred_76.)
Post: #19
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
RE: [HP35s] Program for prime number (brut force)
(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.
02-15-2019, 03:07 PM (This post was last modified: 02-15-2019 05:51 PM by Albert Chan.)
Post: #20
 Albert Chan Senior Member Posts: 943 Joined: Jul 2018
RE: [HP35s] Program for prime number (brut force)
(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)
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)