XProfan

Window Title "Statistischer Miller-Rabin Primzahltest"
' (CL) CopyLeft 2018-05 P.woodpecker, Wien; without each Gewähr!
Windowstyle 24:font 2:set("decimals",17)
Declare w$,x!,y!,m!,n!

Repeat

    print "\n Primzahltest [group = -1]: ";
    input w$ :case w$="-2":w$="90000103"'yet correctly
    'Testprime "10,000,007" not any more recognized!

    if w$>90000103

        print " *** ZUVERLÄSSIGKEIT sink ex 90,000,000 ! ***"
        sound 100,100':continue

    endif

    case w$="-1":Break
    n!=val(w$)
    locate %csrlin-1,50:print " MillerRabin_7: ";
    print if(MillerRabin(n!),"PRIMZAHL","NICHT prim.")

Until %key=27

Cls:set("decimals",0)
declare Anz&,From!,To!,running!
Von!=90000000'
Bis!=90000500'To ca. 90,000,000
case To!>100000000:print "\n long duration expect! \n"
Lauf!=From!+(remodf(of!,2)=0)

REPEAT

    if MillerRabin(running!)

        print running!,
        inc Anz&

    endif

    running!=running!+2

UNTIL running!>(To!+1)

print "\n amount: ",Anz&
'Check: 2871 Primes of 20000 To 50000
'Check: 764 Primes of 520000 To 530000
'Check: 658 Primes of 5020000 To 5030000
'Check: 577 Primes of 50020000 To 50030000
beep:sound 2000,200
waitinput
End

proc MillerRabin :parameters n!

    case n!<2:return 0
    case n!<>intf(n!):return 0
    case n!<4:return 1
    casenot remodf(n!,2):return 0
    declare s&,d!,x!,res&
    s&=0
    d!=n!-1

    whilenot remodf(d!,2)

        d!=d!/2
        inc s&

    endwhile

    res&=1

    WHILELOOP 7

        x! = modpotf(2+intf(rnd()*(n!-3)),d!,n!)
        Case (x!=1) OR (x!=n!-1):continue

        whileloop s&-1,1,-1

            x!=remodf(sqr(x!),n!)

            if x!=1:res&=0:Break:endif

                if x!=n!-1:res&=2:Break:endif

                endwhile

                case res&=0:break

                if res&=2:res&=1:continue:endif

                    case x!=n!-1:continue
                    res&=0

                ENDWHILE

                return res&

            ENDPROC

            Proc modpotf :parameters x!,y!,m!

                var pot! = 1

                while y!>0

                    if remodf(y!,2)=1

                        pot!=remodf(pot!*x!,m!)
                        y!=y!-1

                    else

                        x!=remodf(sqr(x!),m!)
                        y!=y!/2

                    endif

                endwhile

                return pot!

            endproc

            proc floor :parameters x!

                case abs(x!)<(10^-35):return 0
                case x!>0:return intf(x!)
                return (abs(x!-intf(x!)) < 10^-35)-intf(abs(x!-1))

            endproc

            proc remodf :parameters x!,y!

                ' Q: https://de.wikipedia.org/wiki/Modulo , How in ADA
                case abs(x!)<(10^-35):return 0
                case abs(y!)<(10^-35):return x!
                return ((x!>0)-(x!<0))*abs(x!-y!*floor(x!/y!))

            endproc

            proc frac :parameters x!

                var s!=(x!>0)-(x!<0)
                x!=abs(x!)
                x!=x!-round(x!,0)
                case x!<0:x!=1+x!
                return s!*x!

            endproc

            proc intf :parameters x!

                var s!=(x!>0)-(x!<0)
                x!=abs(x!)
                x!=x!-frac(x!)
                return s!*x!

            endproc