Threaded Mode | Linear Mode

fred_76 · (This post was last modified: 04-11-2019 05:25 PM by fred_76.)

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