Threaded Mode | Linear Mode

fred_76 · (This post was last modified: 04-12-2019 12:44 PM by fred_76.)

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

Don Shepherd · (This post was last modified: 02-12-2019 01:09 AM by Don Shepherd.)

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

fred_76 · 02-12-2019, 09:28 AM

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+ !

Don Shepherd · 02-12-2019, 01:54 PM

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

fred_76 · 02-12-2019, 06:27 PM

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

Albert Chan · (This post was last modified: 02-15-2019 01:02 PM by Albert Chan.)

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.

fred_76 · 02-13-2019, 10:09 AM

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

Thomas Klemm · 02-13-2019, 11:53 AM

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

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 !

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

Albert Chan · (This post was last modified: 02-13-2019 07:23 PM by Albert Chan.)

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.

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

fred_76 · (This post was last modified: 02-14-2019 04:49 PM by fred_76.)

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

Albert Chan · 02-14-2019, 06:26 PM

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

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

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.

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

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

fred_76 · (This post was last modified: 02-15-2019 09:06 AM by fred_76.)

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

Albert Chan · (This post was last modified: 02-15-2019 05:51 PM by Albert Chan.)

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