Here is the code of the MR primality check.
I've been inspired by the post of
Thomas Klemm. I haven't used the
Gerald H arithmetic library (far too complex to me to understand with all the imbricated XEQs... sorry Gerald). That leaves place for further optimisation.
Modular addition
Here we reduce first A to A≡N and B to B≡N then calculate
A-N+B (which can't overflow) ≡ N
we could save a calculation of A-N if A+B does not overflow, but the test takes longer than the substraction A-N.
Code:
Line Code Comment
A001 LBL A X=N, Y=B, Z=A returns (A+B)≡N
A002 RMDR
A003 x<>y
A004 Last X
A005 RMDR
A006 Last X
A007 -
A008 Last X
A009 Rv
A010 +
A011 R^
A012 RMDR
A013 RTN
Exemple of use :
RCL A
RCL B
RCL N
XEQ A
returns (A+B)≡N in X and N in Last X (y z t contain rubbish)
Modular multiplication
Based on the fact that if you write :
A = xk+y
B = uk+u
where k is a base number, I took k = (999 999 999 999)^0.5 = 1 000 000
Then
A*B = (xk+y)*(uk+v) = xuk²+xvk+yuk+yv or better = k(xuk+xv+yu)+yv
A*A = x²k²+xyk+xyk+y² = k(x²k+xy+xy)+y² Note that you can't use directly 2xyk as it sometimes overflows.
The code is slightly more complex as it treats differently the case (AxB)≡N and (AxA)≡N for speed.
Code:
Line Code Comment
B001 LBL B X=N, Y=B, Z=A
B002 CF 0
B003 STO N N=module
B004 Rv
B005 x=y?
B006 SF 0
B007 ENTER
B008 ENTER
B009 1000000
B010 STO K k=10…
B011 INT/
B012 STO U
B013 CLx
B014 Last X
B015 RMDR
B016 STO V
B017 FS? 0
B018 GTO a=b
B019 CLx
B020 Last X
B021 INT/
B022 STO X
B023 CLx
B024 Last X
B025 RMDR
B026 STO Y
B027 RCL X A<>B
B028 RCL*U
B029 RCL N
B030 RMDR
B031 RCL*K xuk
B032 RCL X
B033 RCL*V xv
B034 RCL N
B035 XEQ A xv+xuk
B036 RCL Y
B037 RCL*U yu
B038 RCL N
B039 XEQ A yu+xv+xuk
B040 RCL*K k(yu+xv+xuk)
B041 RCL Y
B042 RCL*V yv
B043 RCL N
B044 XEQ A k(yu+xv+xuk)+yv
B045 RTN
B046 RCL U A=B
B047 RCL*V uv
B048 ENTER
B049 ENTER
B050 RCL N
B051 XEQ A 2uv
B052 RCL U
B053 x² u²
B054 RCL N
B055 RMDR
B056 RCL*K ku²
B057 RCL N
B058 XEQ A ku²+2uv
B059 RCL*K k(ku²+2uv)
B060 RCL V
B061 x² v²
B062 RCL N
B063 XEQ A k(ku²+2uv)+v²
B064 CF 0
B065 RTN
Exemple of use :
RCL A
RCL B
RCL N
XEQ B
returns (A*B)≡N in X and N in Last X (y z t contain rubbish)
Modular exponentiation
This one is the regular "right-to-left" modular exponentiation.
Code:
Line Code Comment
C001 LBL C N in X, B in Y, A in Z
C002 STO M
C003 Rv
C004 STO B
C005 Rv
C006 STO A
C007 1
C008 STO R
C009 RCL B
C010 x<=0?
C011 GTO C032
C012 2
C013 RMDR =1 if odd, 0 if even
C014 X>0? if odd
C015 XEQ C026
C016 RCL A
C017 RCL A
C018 RCL M
C019 XEQ B call modular *
C020 STO A A = (A*A) ≡ M
C021 RCL B
C022 2
C023 INT/
C024 STO B
C025 GTO C010
C026 RCL R
C027 RCL A
C028 RCL M
C029 XEQ B call modular *
C030 STO R R = (A*R) ≡ M
C031 RTN
C032 RCL R
C033 RTN
Exemple of use :
RCL A
RCL B
RCL N
XEQ C
returns (A^B)≡N in X (yzt and last x contain rubbish)
Check if composite
As explained
here.
Code:
Line Code Comment
D001 LBL D X=N, Y=D, Z=S, T=A
D002 STO P stores N in P
D003 Rv
D004 STO D stores D in D
D005 Rv
D006 STO S stores S in S
D007 x<>y
D008 RCL P
D009 RMDR A≡ N
D010 RCL D
D011 Last X
D012 XEQ C (A^D) ≡ N
D013 STO Z stores in Z
D014 1
D015 x=y?
D016 GTO D031
D017 1
D018 RCL-P
D019 RCL+Z
D020 x=0?
D021 GTO D031
D022 RCL Z
D023 RCL Z
D024 RCL P
D025 XEQ B X²=N
D026 STO Z
D027 DSE S
D028 GTO D017
D029 1 if true
D030 RTN
D031 0 if false
D032 RTN
Prime routine
I used the bases shown
here (thank you Albert Chan).
Code:
Line Code Comment
E001 LBL E
E002 STO L
E003 1
E004 -
E005 STO E
E006 0
E007 STO T
E008 RCL E
E009 2
E010 RMDR
E011 x>0?
E012 GTO E019
E013 RCL E
E014 2
E015 INT/
E016 STO E
E017 ISG T increment T
E018 GTO E008
E019 RCL L
E020 132239 trigger for 1 base
E021 x>y?
E022 GTO E088
E023 Rv
E024 360018361 trigger for 2 bases
E025 x>y?
E026 GTO E075
E027 Rv
E028 154639673381 trigger for 3 bases
E029 x>y?
E030 GTO E056
E031 SF 4 4 bases
E032 SF 3
E033 SF 2
E034 SF 1
E035 2
E036 XEQ E095
E037 CF 4
E038 x>0?
E039 GTO E100
E040 2570940
E041 XEQ E095
E042 CF 3
E043 x>0?
E044 GTO E100
E045 211991001
E046 XEQ E095
E047 CF 2
E048 x>0?
E049 GTO E100
E050 3749873356
E051 XEQ E095
E052 CF 1
E053 x>0?
E054 GTO E100
E055 GTO E108
E056 SF 3 3 bases
E057 SF 2
E058 SF 1
E059 15
E060 XEQ E095
E061 CF 3
E062 x>0?
E063 GTO E100
E064 176006322
E065 XEQ E095
E066 CF 2
E067 x>0?
E068 GTO E100
E069 4221622697
E070 XEQ E095
E071 CF 1
E072 x>0?
E073 GTO E100
E074 GTO E108
E075 SF 2 2 bases
E076 SF 1
E077 1143370
E078 XEQ E095
E079 CF 2
E080 x>0?
E081 GTO E100
E082 2350307676
E083 XEQ E095
E084 CF 1
E085 x>0?
E086 GTO E100
E087 GTO E108
E088 SF 1 1 base
E089 814494960528
E090 XEQ E095
E091 CF 1
E092 x>0?
E093 GTO E100
E094 GTO E108
E095 RCL T
E096 RCL E
E097 RCL L
E098 XEQ D
E099 RTN
E100 RCL L
E101 CF 1
E102 CF 2
E103 CF 3
E104 CF 4
E105 STO C
E106 VIEW C
E107 RTN
E108 RCL L
E109 STO P
E110 VIEW P
E111 RTN
Exemple of use :
RCL N
XEQ E
shows P= N if N is Prime
shows C= N if N is Composite
the flags 1-2-3-4 show the status of the calculations (4 bases max) so you see where it is
the flag 0 is used by the modular multiplication, it blinks so that you know it did not freeze :-)
One could go slightly faster by storing the constants that are used several times in variables (like 1000000, 1 and 2).
Adding brut force to find factors of for "small numbers"
1) I have also merged the "brut force" with the "MR" so that numbers smaller than 20 359 000 will be analysed by brut force, and greater by MR, because MR is slower than brut force for numbers smaller than 20 359 000.
2) When a number is found to be composite, it goes to brut force to find factors. When a factor is found, the quotient is sent back to MR to check if it is composite or prime, and so on. Brut force is however very slow to find factors if they are big. For example it takes only 34s for MR to say 999 999 998 479 is composite, but 2h30 for Brut Force to find its factors...
Fred