Post Reply 
[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
[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 Time
6 digits       999 863               00
:00:10
7 digits       9 999 991             00
:00:28
8 digits       99 999 989            00
:01:28
9 digits       999 999 937           00
:04:36
10 digits      9 999 999 967         00
:14:37
11 digits      99 999 999 977        00
:46:26
12 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
Find all posts by this user
Quote this message in a reply
02-12-2019, 12:52 AM (This post was last modified: 02-12-2019 01:09 AM by Don Shepherd.)
Post: #2
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 Time
6 digits       999 863               00
:00:09.53
7 digits       9 999 991             00
:00:28.19
8 digits       99 999 989            00
:01:27.70
9 digits       999 999 937           00
:04:36.35
10 digits      9 999 999 967         00
:14:36.59
11 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.
Find all posts by this user
Quote this message in a reply
02-12-2019, 09:28 AM
Post: #3
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+ !
Find all posts by this user
Quote this message in a reply
02-12-2019, 01:54 PM
Post: #4
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.
Find all posts by this user
Quote this message in a reply
02-12-2019, 06:27 PM
Post: #5
RE: Hello and program for prime number (HP35s)
I just found that archived discussion :

http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=169332

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).
Find all posts by this user
Quote this message in a reply
02-12-2019, 08:14 PM (This post was last modified: 02-15-2019 01:02 PM by Albert Chan.)
Post: #6
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.
Find all posts by this user
Quote this message in a reply
02-13-2019, 10:09 AM
Post: #7
RE: Hello and program for prime number (HP35s)
Albert

Good try. I tried, replacing the lines :

Code:
CLx
LAST 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".
Find all posts by this user
Quote this message in a reply
02-13-2019, 11:53 AM
Post: #8
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

This article gives a bit more details:
Miller-Rabin Primality Test for the HP-48

Cheers
Thomas
Find all posts by this user
Quote this message in a reply
02-13-2019, 04:03 PM
Post: #9
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 !
Find all posts by this user
Quote this message in a reply
02-13-2019, 04:44 PM
Post: #10
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...
Find all posts by this user
Quote this message in a reply
02-13-2019, 05:26 PM (This post was last modified: 02-13-2019 07:23 PM by Albert Chan.)
Post: #11
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.
Find all posts by this user
Quote this message in a reply
02-13-2019, 07:20 PM
Post: #12
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
Find all posts by this user
Quote this message in a reply
02-14-2019, 01:10 PM (This post was last modified: 02-14-2019 04:49 PM by fred_76.)
Post: #13
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>=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<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)
Find all posts by this user
Quote this message in a reply
02-14-2019, 06:26 PM
Post: #14
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
Find all posts by this user
Quote this message in a reply
02-14-2019, 06:48 PM
Post: #15
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. Smile

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.
Visit this user's website Find all posts by this user
Quote this message in a reply
02-14-2019, 07:51 PM
Post: #16
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.
Find all posts by this user
Quote this message in a reply
02-14-2019, 08:50 PM
Post: #17
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
Find all posts by this user
Quote this message in a reply
02-14-2019, 09:17 PM
Post: #18
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
Find all posts by this user
Quote this message in a reply
02-15-2019, 08:47 AM (This post was last modified: 02-15-2019 09:06 AM by fred_76.)
Post: #19
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.
Find all posts by this user
Quote this message in a reply
02-15-2019, 03:07 PM (This post was last modified: 02-15-2019 05:51 PM by Albert Chan.)
Post: #20
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)
Find all posts by this user
Quote this message in a reply
Post Reply 




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