Autor Beitrag
Boldar
ontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic starofftopic star
Beiträge: 1555
Erhaltene Danke: 70

Win7 Enterprise 64bit, Win XP SP2
Turbo Delphi
BeitragVerfasst: Fr 28.08.09 23:14 
Hallo,
Ich habe versucht den Solovay-Strassen-Test mittels Delphi und Bignum2 zu implementieren.
Allerding kriege ich eine AV in der Prozedur BM_GCD.(Das Benbe das nich abgefangen hat...)
Woran liegt das?
Und: Ich wäre auch ganz glücklich, wenn mal einer von den experten hier sich den Code in Bezug auf Optimierungsmöglichkeiten ansehen könnte...
mfg und Vielen Dank Boldar

Quelltext:

ausblenden volle Höhe Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
unit primes;

interface

uses bignum2, dialogs, sysutils;

function isprime (z: TBignumber; rounds: integer):boolean;

implementation
               
function bigrandom (a: TBignumber):TBignumber;
var s: string;
  I: Integer;
var stellen: integer;
begin
  setlength(s, 0);
  stellen := length(BMD_bigNumtostr(a, false));
  for I := 0 to stellen - 1 do
  begin
    s := s + inttostr(random(10));
  end;
end;

function jacobi (a, b: TBignumber):TBignumber;forward;

function isprime (z: TBignumber; rounds: integer):boolean;
var a, b, c, j, exp, g, _1, _2, _3: TBignumber;
var _1kb, _bkn_1: boolean;
var i: integer;
begin
  result := true;
  _1 := BMD_strtobignum('1', false);
  _2 := BMD_strtobignum('2', false);
  _3 := BMD_strtobignum('3', false);
  i:=0;
  repeat
  inc (i);
  a := bigrandom (z);
  g := BM_GCD (a, z);
  if (BM_CompareNE (_1,g)) then
  begin
    result := false;
    exit;
  end;
//////////////////////////////////////////////////////////// 1 < b < n-1,
  exp := BM_sub (z, _1);
  exp := BM_divide (exp, _2);
  c := BM_Power(a, exp);
  b := BM_modulo (c, z);
  _1kb := BM_CompareC(_1, b);
  _bkn_1 := BM_CompareC(b, BM_sub (z, _1));
  if _1kb and _bkn_1 then
  begin
    result := false;
    exit;
  end;
  j:= jacobi (a, z);
  if BM_CompareNe(j, BM_Modulo(b, z)) then
  begin
    result := false;
    exit;
  end;
  until i>rounds;

end;

function jacobi (a, b: TBignumber):TBignumber;
function even (a: TBignumber):boolean;inline;
var _0, _1, _2, _3: TBignumber;
begin
  result := false;
  _0 := BMD_strtobignum('0', false);
  _1 := BMD_strtobignum('1', false);
  _2 := BMD_strtobignum('2', false);
  _3 := BMD_strtobignum('3', false);
  if BM_CompareZ(BM_modulo(a, _2), _0) then result := true;
end;
var _m1, _0, _1, _2, _3, _4, _5, _6, _7, _8: TBignumber;
var g: TBignumber;
begin
  _m1:= BMD_strtobignum('-1', false);
  _0 := BMD_strtobignum( '0', false);
  _1 := BMD_strtobignum( '1', false);
  _2 := BMD_strtobignum( '2', false);
  _3 := BMD_strtobignum( '3', false);
  _4 := BMD_strtobignum( '4', false);
  _5 := BMD_strtobignum( '5', false);
  _6 := BMD_strtobignum( '6', false);
  _7 := BMD_strtobignum( '7', false);
  _8 := BMD_strtobignum( '8', false);
  if even(b) then
  begin
    showmessage ('Fehler: b grade');
    exit;
  end;
  if BM_CompareNC (a, b) then a := BM_Modulo(a, b);
  if BM_CompareZ  (a, _0) then
  begin
    result := _0;
    exit;
  end;
  if BM_CompareZ  (a, _1) then
  begin
    result := _1;
    exit;
  end;

  if BM_CompareC (a, _0) then if BM_Comparez(BM_Modulo(BM_divide(BM_sub(b, _1), _2), _2), _0) then
  begin
    result := (jacobi (BM_Multiply(a, _m1), b));
    exit;
  end
  else
  begin
    result :=BM_Multiply(jacobi (BM_Multiply(a, _m1), b), _m1);
    exit;
  end;

  if even(a) then
  begin
    if BM_CompareZ(BM_Modulo(BM_divide(BM_Multiply(b, bm_sub(b, _1)), _8), _2), _0) then result := jacobi(BM_divide(a, _2),b) else result := BM_Multiply(jacobi(BM_divide(a, _2),b), _m1);
  end;
  g := BM_GCD (a, b);
  if BM_CompareZ(g, a) then
  begin
    result := _0;
    exit;
  end
  else
  begin
  if BM_CompareNE(g, _1) then
  begin
    result := BM_Multiply(jacobi(g, b), jacobi(bm_divide(a, g), b));
    exit;
  end
  else
  begin
  if BM_CompareZ(BM_Modulo(BM_Divide(BM_Multiply(BM_sub(a, _1), BM_Sub(b, _1)), _4), _2), _0) then
  begin
    result := jacobi(b, a);
    exit;
  end
  else
  begin
    result := BM_Multiply(jacobi(b, a), _m1);
  end;
  end;
end;
end;

end.
Gammatester
ontopic starontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic star
Beiträge: 328
Erhaltene Danke: 101



BeitragVerfasst: Sa 29.08.09 00:08 
Habe allerdings selten einen so ineffektiven Code gesehen. Zum Bignum2-GCD-Problem kann und will ich nichts sagen
(eigentlich solltest Du es doch mit dem Debugger raus finden können).

Zur Optimierung:

1. Wer benutzt denn noch Solovay-Strassen? Miller-Rabin ist viel effektiver. Empfehlung: nimm BPSW (Baillie-Pomerance-Selfridge-Wagstaff).

2. Jacobi kann nur die Werte -1,0,1 annehmen. Warum man dann dafür TBignumber als Result-Typ hernimmt, ist mit völlig schleierhaft.

3. Jacobi rekursiv ist höchstens was für die Demo-Vorlesung

4. Wenn man even() nicht schneller berechnen kann, solltest Du vielleicht eine andere MP-Arith-Bibliotek benutzen (was sollen die Variable
_0, _1, _2, _3 in even, die jedesmal erzeugt, aus einem String berechnet, vernichtet werden. Eventuell falls es wirklich nichts besseres gibt, _2 einmal
erzeugen in jacobi nicht in even und natürlich jacobi nicht rekursiv.
BenBE
ontopic starontopic starontopic starontopic starontopic starontopic starhalf ontopic starofftopic star
Beiträge: 8721
Erhaltene Danke: 191

Win95, Win98SE, Win2K, WinXP
D1S, D3S, D4S, D5E, D6E, D7E, D9PE, D10E, D12P, DXEP, L0.9\FPC2.0
BeitragVerfasst: Sa 29.08.09 16:53 
Miller-Rabin geht über die BM_ExpMod-Funktion recht simpel zu realisieren

Odd\Even abfragen einfach im Datenbereich das niedrigste Bit abfragen ... (Wenn Exp=0; ansonsten Even)

_________________
Anyone who is capable of being elected president should on no account be allowed to do the job.
Ich code EdgeMonkey - In dubio pro Setting.