[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
 fred_76 Junior Member Posts: 33 Joined: Feb 2019
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.

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
 « Next Oldest | Next Newest »

 Messages In This Thread [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) - fred_76 - 02-11-2019, 05:16 PM RE: Hello and program for prime number (HP35s) - Don Shepherd - 02-12-2019, 12:52 AM RE: [HP35s] Program for prime number (brut force) - Dave Britten - 02-14-2019, 06:48 PM RE: Hello and program for prime number (HP35s) - fred_76 - 02-12-2019, 09:28 AM RE: Hello and program for prime number (HP35s) - Don Shepherd - 02-12-2019, 01:54 PM RE: Hello and program for prime number (HP35s) - fred_76 - 02-12-2019, 06:27 PM RE: Hello and program for prime number (HP35s) - Albert Chan - 02-12-2019, 08:14 PM RE: Hello and program for prime number (HP35s) - fred_76 - 02-13-2019, 10:09 AM RE: [HP35s] Hello and program for prime number - Thomas Klemm - 02-13-2019, 11:53 AM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 10:45 AM RE: [HP35s] Hello and program for prime number - fred_76 - 02-13-2019, 04:03 PM RE: [HP35s] Hello and program for prime number - fred_76 - 02-13-2019, 04:44 PM RE: [HP35s] Hello and program for prime number - Albert Chan - 02-13-2019, 05:26 PM RE: [HP35s] Hello and program for prime number - Thomas Klemm - 02-13-2019, 07:20 PM RE: [HP35s] Hello and program for prime number - fred_76 - 02-14-2019, 01:10 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-15-2019, 03:07 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-14-2019, 06:26 PM RE: [HP35s] Program for prime number (brut force) - fred_76 - 02-14-2019, 07:51 PM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-14-2019, 08:50 PM RE: [HP35s] Program for prime number (brut force) - fred_76 - 02-15-2019, 08:47 AM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-14-2019, 09:17 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 11:16 AM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-16-2019, 11:36 AM RE: [HP35s] Program for prime number (brut force) - fred_76 - 02-16-2019, 11:48 AM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-16-2019, 11:57 AM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 12:54 PM RE: [HP35s] Program for prime number (brut force) - fred_76 - 02-16-2019, 02:22 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-16-2019, 02:29 PM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-16-2019, 02:12 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 05:39 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 05:43 PM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-16-2019, 05:57 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 05:59 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-16-2019, 10:32 PM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-16-2019, 06:07 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 06:18 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-16-2019, 07:09 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-16-2019, 07:57 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-17-2019, 05:22 AM RE: [HP35s] Program for prime number (brut force) - Thomas Klemm - 02-16-2019, 08:39 PM RE: [HP35s] Program for prime number (brute force) - Thomas Klemm - 02-17-2019, 04:15 AM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-17-2019, 05:46 AM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-17-2019, 04:33 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 02-17-2019, 10:58 PM RE: [HP35s] Program for prime number (brut force) - Gerald H - 02-18-2019, 06:09 AM RE: [HP35s] Program for prime number (brut force) - fred_76 - 04-09-2019, 03:52 PM RE: [HP35s] Program for prime number (brut force) - ttw - 04-09-2019, 05:22 PM RE: [HP35s] Program for prime number (brut force) - fred_76 - 04-10-2019, 07:26 AM RE: [HP35s] Program for prime number (brut force) - fred_76 - 04-10-2019, 03:56 PM RE: [HP35s] Program for prime number (brut force) - Albert Chan - 04-11-2019, 02:52 AM RE: [HP35s] Program for prime number (brut force) - ttw - 04-10-2019, 09:24 PM RE: [HP35s] Program for prime number (brut force) - ttw - 04-11-2019, 07:42 AM RE: [HP35s] Program for prime number (brut force) - fred_76 - 04-11-2019 03:18 PM RE: [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) - Albert Chan - 04-13-2019, 06:24 PM RE: [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) - Albert Chan - 04-15-2019, 03:55 PM RE: [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin) - Gerald H - 08-22-2022, 10:25 AM

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