Naja St.Pauli, dann hab ich mal deine Hausaufgaben gemacht.
In der Funktion Witness heißt es
Zitat:
for i <-- 1 to t
do blablabla
if blablabal
then return true
Ich kenne mich mit Pseudocodes nicht aus, interpretiere das aber so, daß wenn "if blablabla" True ergibt, True zurückgegeben werden soll. D.h. es müßte result=true gesetzt und die Funktion verlassen werden. Du setzt zwar result=true verläßt aber nicht die Funktion. Ich hab das mal probiert, brachte aber auch nichts. Herr MillerRabin sagte immer noch 23 sei composite.
Also weitergesucht.
In der Funktion ModularExponentiation heit es im Pseudo Code
Zitat:
Let (bk, bk-1,...b0) be the binary representation of b
for i <-- k downto 0
Bei Dir wird das so umgesetzt daß, b als Binärstring in b_bin gestellt wird.
b_bin[Length(b_bin)] ist dabei Bit0.
Auf den ersten Blick sieht das so aus, wie im Pseudo Code, aber beim zweiten Blick eben nicht.
Die Schleife for i:=k downto 0 faßt zuerst hast höchstwertige Bit an und zuletzt Bit0.
Bei Deiner Schleife for i:=Length(b_bin) downto 1 ist es genau umgekehrt.
Ich hab das korrigiert und dann für alle (ungeraden) Zahlen von 3 bis 65535 probiert. Bis 46337 bringt das Programm korrekte Ergebnisse. Für höhere Werte sagt es nur noch composite.... Hab nicht weiter geprüft, warum.
Und weil Du Dich für Primzahlen zu interessieren scheinst ist hier noch eine kleine Routine,
die Extendend-Werte daraufhin prüft, ob sie prim sind. Die zu prüfende Zahl darf natürlich nicht größer als 2^64-1 sein, also max 64 Bit Werte.
Delphi-Quellcode:
{------------------------------------------------------------------------------}
{ FIsPrimeNPA }
{ Returns true if p is prime. NO ABORT-FUNCTIONALITY, NO PROGRESSBAR }
{ Argument : p = the number to test. }
{ Result : 0 = p is not prime - then the smallest divisor is returned in r }
{ 1 = p is prime }
{ p can be any extended number. }
{ the routine tests if it is a positive integer and returns false if not. }
{------------------------------------------------------------------------------}
FUNCTION FIsPrimeNPA(p:extended;
var r:extended):integer;
asm
mov edx,eax
// save pointer to r
sub esp,12
// reserve space on stack
fnstcw [ESP].Word
// save control word
fnstcw [ESP+2].Word
// and once more
fwait
or [ESP+2].Word,$0F00
// trunc decimals, full precision
fldcw [ESP+2].Word
// load new control word
fldz
// ST0=0
fstp tbyte [edx]
// as default store 0 as largest divisor
fldz
// ST0=0
fld p
// ST0=p ST1=0
fld1
// ST0=1 ST1=p ST2=0
fadd st,st
// ST0=2 ST1=p ST2=0
fcom st(1)
// compare 2 with p
fstsw ax
// status word into AX
sahf
// copy ah into flag register
je @True
// p=2, is prime
jnc @False2
// p<2, is not prime
// we now know that p>2 i.e. positive
// let's check if it is integer and odd
fld st(1)
// ST0=p ST1=2 ST2=p ST3=0
frndint
// truncate decimals
fcom st(2)
// compare with p
fstsw ax
// status word into AX
sahf
// copy ah into flag register
jne @False2
// p is no integer
fdiv st,st(1)
// divide p by 2
frndint
// truncate decimals
fmul st,st(1)
// multiply by 2
fcomp st(2)
// compare with p and pop from stack
fstsw ax
// status word into AX
sahf
// copy ah into flag register
je @False1
// p is even
// we know now that p is an odd integer > 2
// now let's find the maximum divisior for prime test
fld st(1)
// ST0=p ST1=2 ST2=p
fsqrt
// ST0=Sqrt(p) ST1=2 ST2=p
frndint
// truncate decimals
fistp qword ptr [ESP+4]
// save int(sqrt(p)) and remove from stack
fwait
mov ecx,[ESP+4]
// ECX=int(sqrt(p))
sub ecx,1
// subtract (smallest divisor - 2)
shr ecx,1
// ECX=number of divisors
je @True
// p is prime (3, 5 or 7)
// now we can start the prime test
fld1
// ST0=1(divisor) ST1=2 ST2=p ST3=0
// --------------------------------------------------------------
// INCREASE DIVISOR BY 2
// --------------------------------------------------------------
@Loop: fadd st,st(1)
// ST0=ST0+2
// --------------------------------------------------------------
// TEST IF P IS DIVISIBLE BY DIVISOR
// --------------------------------------------------------------
fld st(2)
// ST0=p ST1=divisor ST2=2 ST3=p ST4=0
fprem
// ST0=p mod d ST1=divisor ST2=2 ST3=p ST4=0
fcomp st(4)
// ST0=0 ? - ST0=divisor ST1=2 ST2=p ST3=0
fstsw ax
// control word into AX
sahf
// copy ah into flag register
je @False1
// p is divisible by d, i.e. p is not prime
sub ecx,1
// counter - 1
jne @Loop
// loop until counter = 0
@True: mov eax,1
// result=true
jmp @
End
@False1: fstp tbyte [edx]
@False2:
xor eax,eax
// result=false
//---------------------------------------------------------------
// finally clear fpu-stacks, reset control word and ESP
// --------------------------------------------------------------
@
End: ffree st(3)
ffree st(2)
ffree st(1)
ffree st(0)
fldcw [ESP].Word
add esp,12
end;
Gruß, Klaus