XProfan

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

Repeat

    print "\n Primzahltest [Gruppe = -1]: ";
    input w$ :case w$="-2":w$="90000103"'noch korrekt
    'Testprime "10,000,007" nicht mehr erkannt!

    if w$>90000103

        print " *** ZUVERLÄSSIGKEIT sinkt ab 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&,Von!,Bis!,Lauf!
Von!=90000000'
Bis!=90000500'bis ca. 90,000,000
case Bis!>100000000:print "\n Lange Dauer erwarten! \n"
Lauf!=Von!+(remodf(von!,2)=0)

REPEAT

    if MillerRabin(Lauf!)

        print Lauf!,
        inc Anz&

    endif

    Lauf!=Lauf!+2

UNTIL Lauf!>(Bis!+1)

print "\n Anzahl: ",Anz&
'Check: 2871 Primes von 20000 bis 50000
'Check: 764 Primes von 520000 bis 530000
'Check: 658 Primes von 5020000 bis 5030000
'Check: 577 Primes von 50020000 bis 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 , wie 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