Post Reply 
[HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
04-11-2019, 03:18 PM (This post was last modified: 04-11-2019 05:25 PM by fred_76.)
Post: #53
RE: [HP35s] Program for prime number (brut force)
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.

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


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.

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

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 :
returns (A^B)≡N in X (yzt and last x contain rubbish)

Check if composite

As explained here.

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

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

Find all posts by this user
Quote this message in a reply
Post Reply 

Messages In This Thread
RE: [HP35s] Program for prime number (brut force) - fred_76 - 04-11-2019 03:18 PM

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