Primfaktorzerlegung optimieren

Horst_H · Verfasst: Mo 21.11.05 17:17

Hallo,
es dauert mit +2,+4 bei mir

Quelltext

1:
2:

9223372036854775806  = 2 * 3 * 715827883 * 2147483647
00:01:22.796

wenn ich das DiffFeld mod 2,3,5,7,11,13 einsetze sind es:

markieren

Quelltext

1:
2:

9223372036854775806  = 2 * 3 * 715827883 * 2147483647
00:00:46.031 Sekunden

das ist doch ein wenig besser, aber eben nicht optimal.

Gruss Horst

Born-to-Frag · Beiträge: 1094 Win XP SP2, Win 2000 SP4 Delphi 7, 2k5

Horst_H hat folgendes geschrieben:

Hallo,
ist das Ergebnis:
9223372036854775805 = 5*23 * 53301701 * 1504703107?
Dass dauert 9.5 Sekunden.(Mit DiffFeld Mod 2,3,5,7,11,13 nur 3.7 Sekunden)

Gute Guete wer rechnet denn da staendig eine Wurzel ausserhalb der while Schleife, dass ist mir beim lesen nicht aufgefallen.

markieren

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:

function PrimFaktorZerlegung (Zahl: Int64): String;
var
  Teiler,
  Root: Cardinal;
  Index : integer;
  begin
  Result := '';
  While Zahl and 1 = 0 do
    Begin
    Result := Result + IntToStr(2) + ' * ';
    Zahl := Zahl shr 1;
    end;
  {
   For Index := 0 to BIS do
     While Zahl Mod and Prim[Index] = 0 do
       Begin
       Result := Result + IntToStr(Prim[Index]) + ' * '; 
       Zahl := Zahl DIV Prim[index];
       end;
   }
  if zahl >= 2 then
    Begin
    Root := Trunc(SQRT(Zahl));
    Teiler := 1;
    Index := 0;
    repeat
      inc(Teiler,2);
      //inc(Teiler,DiffFeld[index]);
      IF Zahl mod Teiler = 0 then
        begin
        repeat
          Zahl := Zahl div Teiler;
          Result := Result + IntToStr(Teiler) + ' * ';
        until Zahl mod Teiler <> 0; 
        Root := Trunc(SQRT(Zahl));//Nur dann wenn sich wirklich was aendert
        end;
//      inc(index);
//      If index > High(DiffFeld) then  
//         index := 0;
    until Teiler > Root;
    end;
    //dann jetzt aber doch
    if zahl <> 1 then
      begin
      Result := Result + IntToStr(Zahl) + ' * ';;
      end;
    Delete(Result, Length(Result) - 2, 3);
end;

Gruss Horst

Reden wir hier vom gleichen? Also ich komme mit dieser funktion auf ~190sek

Und mit der optimieren auf ~131sekunden..

_________________
Theorie ist wenn man alles weiß, aber nichts funktioniert. Praxis ist wenn alles funktioniert, aber niemand weiß warum.
Microsoft vereint Theorie und Praxis: Nichts funktioniert und niemand weiß warum.

BenBE · Verfasst: Mo 21.11.05 18:03

k, hab jetzt n Algo mit folgender Ausgabe:

markieren

Quelltext

1:
2:
3:
4:
5:
6:
7:

Zu faktorisierende Zahl: 9223372036854775806
2
3
715827883
2147483647
Faktorisierung abgeschlossen!
Benötigte Zeit: 85,599047 s

Source dafür sieht so aus:

markieren

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:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:

Program WheelCalc;

{$APPTYPE CONSOLE}

Uses
    Windows,
    SysUtils,
    Classes,

    //Omorphia Specific
    ODbgInterface, ODbgMapfile, OIncTypes;

Var
    Factorize: Int64;

Const
    PrimeCount = 8;
    //    Primes: Array[0..PrimeCount - 1] of Integer = (2,3,5,7,11,13,17,19,23,29,31);

Var
    Wheel: Array Of Integer;
    WheelLen: Integer;
    WheelPeriod: Int64;

    X, X2: Integer;
    Y, Z: Int64;
    CurrPrime: Integer;

    WheelPtrSrc, WheelPtrDest: Integer;

//Omorphia Specific:
Procedure HandleError(Level: Integer; Place: TODbgMapLocation; Desc: String; Var ErrorObj: Exception; Var Handled: Boolean);
Var
    SL: TStringList;
Begin
    WriteLn('Exception Level ', Level, ' occured at the following location:');
    WriteLn('  ', PlaceToLocationStr(Place));
    WriteLn('The Stacktrace is as follows:');
    SL := GetStackTrace(10);
    While SL.Count > 0 Do
    Begin
        WriteLn('- ', SL[0]);
        SL.Delete(0);
    End;
End;

Var
    T1, T2, T3: Int64;

Begin
    AddErrorHandler(HandleError); //Omorphia Specific

    Write('Zu faktorisierende Zahl: ');
    ReadLn(Factorize);

    If Factorize < 0 Then
    Begin
        WriteLn('-1');
        Factorize := -Factorize;
    End;
    If Factorize < 2 Then
    Begin
        If Factorize = 0 Then
            Write('0');
        Exit;
    End;

    QueryPerformanceCounter(T1);
    
    X2 := 0;
    While Factorize And 1 = 0 Do
    Begin
        Factorize := Factorize Shr 1;
        Inc(X2);
    End;
    If X2 <> 0 Then
    Begin
        If X2 = 1 Then
            WriteLn('2')
        Else
            WriteLn('2 ^ ', X2);
    End;
    If Factorize = 1 Then
        Break;

    //Initialize the wheel
    WheelLen := 1;
    SetLength(Wheel, WheelLen);
    Wheel[0] := 2;
    WheelPeriod := 2;

    //Build up the wheel
    For X := 0 To PrimeCount - 1 Do
    Begin
        //Get the current prime for extension
        //CurrPrime := Primes[X];
        //Alternatively you could also write
        CurrPrime := 1 + Wheel[0];

        X2 := 0;
        While Factorize Mod CurrPrime = 0 Do
        Begin
            Factorize := Factorize Div CurrPrime;
            Inc(X2);
        End;
        If X2 <> 0 Then
        Begin
            If X2 = 1 Then
                WriteLn(CurrPrime)
            Else
                WriteLn(CurrPrime, ' ^ ', X2);
        End;
        If Factorize = 1 Then
            Break;

        WheelPeriod := WheelPeriod * CurrPrime;

        //Extend the whell to include as many field as required at first.
        SetLength(Wheel, WheelLen * CurrPrime);
        For X2 := 1 To CurrPrime - 1 Do
            Move(Wheel[0], Wheel[WheelLen * X2], WheelLen * SizeOf(Wheel[0]));

        //Start from the very beginning ;-)
        WheelPtrSrc := 0;
        WheelPtrDest := 0;
        Y := 1;

        //Until we finished checking all numbers of the current wheel
        While WheelPtrSrc < WheelLen * CurrPrime Do
        Begin
            Z := 0;
            Repeat
                Y := Y + Wheel[WheelPtrSrc];
                Z := Z + Wheel[WheelPtrSrc];
                Inc(WheelPtrSrc);
            Until (Y Mod CurrPrime <> 0) Or (WheelPtrSrc >= WheelLen * CurrPrime);
            Wheel[WheelPtrDest] := Z;
            Inc(WheelPtrDest);
        End;

        WheelLen := WheelPtrDest;
        SetLength(Wheel, WheelLen);
    End;

    WheelPtrSrc := 1;
    Y := 1 + Wheel[0] + Wheel[1];
    Z := Trunc(Sqrt(1.0 * Factorize));
    While Y < Z Do
    Begin
        X2 := 0;
        While Factorize Mod Y = 0 Do
        Begin
            Factorize := Factorize Div Y;
            Inc(X2);
        End;
        If X2 <> 0 Then
        Begin
            If X2 = 1 Then
                WriteLn(Y)
            Else
                WriteLn(Y, ' ^ ', X2);

            Z := Trunc(Sqrt(1.0 * Factorize));
        End;
        If Factorize = 1 Then
            Break;

        Inc(WheelPtrSrc);
        If WheelPtrSrc >= WheelLen Then
            WheelPtrSrc := 0;
        Y := Y + Wheel[WheelPtrSrc];
    End;

    If Factorize <> 1 Then
        WriteLn(Factorize);

    QueryPerformanceCounter(T2);
    QueryPerformanceFrequency(T3);

    WriteLn('Faktorisierung abgeschlossen!');
    WriteLn(Format('Benötigte Zeit: %.6f s', [(T2 - T1) / T3]));
    ReadLn;
End.

Die mit //Omorphia Specific gekennzeichneten Teile hab ich aus Debugging-Gründen aufgenommen, damit ich im Fehlerfalle die genaue Position des Fehlers leichter finden kann ;-)

Z.B. darf PrimeCount MAXIMAL 9 sein, weil das Programm ansonsten mit einem RangeError abstürzt (weil das Wheel-Array mehr als 4GB benötigen würde).

//Nachtrag:

markieren

Quelltext

1:
2:
3:
4:
5:
6:
7:

Zu faktorisierende Zahl: 9223372036854775805
5
23
53301701
1504703107
Faktorisierung abgeschlossen!
Benötigte Zeit: 10,967501 s

_________________

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.

Born-to-Frag · Beiträge: 1094 Win XP SP2, Win 2000 SP4 Delphi 7, 2k5

Kann ich nicht testen

Hab ja die ganzen uses nicht..

greetz

_________________
Theorie ist wenn man alles weiß, aber nichts funktioniert. Praxis ist wenn alles funktioniert, aber niemand weiß warum.
Microsoft vereint Theorie und Praxis: Nichts funktioniert und niemand weiß warum.

Horst_H · Verfasst: Mo 21.11.05 19:04

Hallo,

das ist naturidentisch mit der Version die ich Dir gesendet habe, nur wesentlich uebersichtlicher.Damals wollte ich auch Platz sparen (29 Mb fuer die Zahlen bis 1e9).
Gruss Horst

BenBE · Verfasst: Mo 21.11.05 19:09

Ich schrieb doch: Die durch Kommentare markierten Teile einfach auskommentieren. Das ist nur mein Error-Handler sowie eine Unit, die sich drt mit reinhängt.

_________________

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.

Born-to-Frag · Beiträge: 1094 Win XP SP2, Win 2000 SP4 Delphi 7, 2k5

Ok, mein fehler Ben

Ich sehs mir noch mal an.. auch deine Version Horst

_________________
Theorie ist wenn man alles weiß, aber nichts funktioniert. Praxis ist wenn alles funktioniert, aber niemand weiß warum.
Microsoft vereint Theorie und Praxis: Nichts funktioniert und niemand weiß warum.

alzaimar · Verfasst: Mo 21.11.05 20:37

Wohin die Reise gehen kann, zeigt dieser Link:
www.thorstenreinecke.de/qsieve/

Ein schnelles Programm in C mit Sourcen (starker Tobak).

_________________
Na denn, dann. Bis dann, denn.

BenBE · Verfasst: Mo 21.11.05 20:46

alzaimar hat folgendes geschrieben:

Ein schnelles Programm in C mit Sourcen (starker Tobak).

Weil das Programm schnell ist, oder weil's mit Sources ist (SCNR)??? :mrgreen:

_________________

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.

Born-to-Frag · Beiträge: 1094 Win XP SP2, Win 2000 SP4 Delphi 7, 2k5

Sind beide extrem kompliziert.. Noch ne Frage: Meine IsPrimzahl function ist doch Sinnlos geworden?! Denn wenn die Zahl keine Primzahl ist, wurde die function doppelt durchlafen. Und selbst wenn sie eine ist, würde das meine PrimFaktorZelegung - function schneller herausfinden!

Oder täusche ich mich da?

greetz

PS: Ich probiers morgen nochmal genauer aus, heut min ich zu müde

_________________
Theorie ist wenn man alles weiß, aber nichts funktioniert. Praxis ist wenn alles funktioniert, aber niemand weiß warum.
Microsoft vereint Theorie und Praxis: Nichts funktioniert und niemand weiß warum.

Horst_H · Verfasst: Di 22.11.05 11:48

Hallo,

IsPrimzahl ist nun obsolet geworden, da es doppelte Rechnerei ist

Wie ich schon weiter oben gezeigt habe, ist die Berechnung des Restes extrem langsam mit ueber 700 CPU Takten.
Wenn man das Programm in der CPU Ansicht verfolgt erkennt man das @_llDIV tatsaechlich Bitweise arbeitet und deshalb alles sooooo laaaangsaaaaamm.
Damit es etwas schneller wird, nicht umsonst benutzt das qsieve Programm auch die GMP - Bibliotheken, etwas problemspezifischeres.

markieren

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:

function Int64divModCar(Teiler : Cardinal;var Dividend,Quotient:Int64):Cardinal;assembler;
//berechnet Int64 divmod Cardinal
//Gibt Rest zurueck
//EAX= Teiler ; EDX =Adresse von Zahl1; ECX=Adresse von Zahl2
asm
  Push EBX
  MOV  EBX,EAX
  MOV  EAX,[EDX+4]
  Push EDI
  MOV  EDI,EDX
  XOR  EDX,EDX
  TEST EAX,EAX
  MOV  [ECX+4],EDX
  JE   @einstellig
  DIV  EBX
  MOV  [ECX+4],EAX
@einstellig:
  MOV  EAX,[EDI]
  POP  EDI
  DIV  EBX
  POP  EBX
  MOV  [ECX],EAX
  MOV  EAX,EDX;
end;

//und angewandt in der Primfaktorzerlegung:
  
  if Int64divModCar(Teiler,TestZahl,Quotient) = 0 then
//    If Zahl mod Teiler = 0 then
      begin
      repeat
        Result := Result + IntToStr(Teiler) + ' * ';
        Zahl := Quotient;
      until INt64divModCar(Teiler,TestZahl,Quotient) <> 0;
      Root := Trunc(SQRT(Zahl));
      end;

//zum testen:
procedure TForm1.Button2Click(Sender: TObject);
var
  i,k : cardinal;
  Dividend,
  Quotient : int64;
  T : tDateTime;
begin
  Dividend := 1001*65535;
  Dividend := Dividend *65537;
  i := 1;
  t := time;
  repeat
    k := INt64divModCar(i,Dividend,Quotient);
    inc(i);
  until i > 100000000;
  t := time-t;
  memo1.Lines.Add(FormatdateTime('hh:mm:ss.zzz',t));
end;

Damit braucht es bei mir Primfaktorzerlegung von ????...06 statt 46 Sekunden noch 6,6.
Die 1e8 Durchlaeufe von Button2Click dauerten 4.61 Sekunden =>
4.61/1e8* 1.953e9 = 90 Takte fuer einen Schleifendurchlauf statt 734 bringt dass schon einen Faktor 8

Gruss Horst
P.S.:
9223372036854775806 = 2 * 3 * 715827883 * 2147483647 Anzahl Divisionen: 137.301.661 max Teiler: 715827883
Das bedeutet bei 6.6 Sekunden == 94 Takte pro Durchlauf, sagenhaft nah an der reinen minimalste Schleife.
Bis 715827883 sind es 37.029.521 Primzahlen
Damit lohnt sich mein vorgeschlagenes Verfahren, innerhalb einer Primzahlsiebschleife den Test durchzufuehren, auch erst, wenn die Bestimmung pro Primzahl unter 2.7(=137Mio/37mio-1) *90 Takte= 240 sinkt.
Mit alzaimar's Sieb ala Atkin's sind es 57 Takte ! (1.5 Sekunden fuer die ersten 50.8 Mio Primzahlen).
Theoretische minimalste Zeit ~ 37e6*(57(finden)+90(testen))(Takte)= 5.439e9 Takte oder 2.78 Sekunden (div f-CPU)
Das ist nochmals eine Beschleunigung um das 2.37 fache.
Damit ist fuer mich das Ende der Fahnenstange fuer einfache Probedivision erreicht.
Jetzt hilft nur noch der Uebergang auf andere Algorithmen
de.wikipedia.org/wiki/Pollard-Rho-Methode ,da steht zwar etwas von wenigen Durchlaeufen, aber GCD bekommt man auch nicht geschenkt.

Born-to-Frag · Beiträge: 1094 Win XP SP2, Win 2000 SP4 Delphi 7, 2k5

Ja gut, mit Assemblern kenn ich mich noch nich so gut aus :oops:

Muss mich mal reinarbeiten

greetz

_________________
Theorie ist wenn man alles weiß, aber nichts funktioniert. Praxis ist wenn alles funktioniert, aber niemand weiß warum.
Microsoft vereint Theorie und Praxis: Nichts funktioniert und niemand weiß warum.

Horst_H · Verfasst: Do 24.11.05 14:39

Hallo,

Hier gibt es eine Unit BigInt, die fast komplett in Assembler geschrieben ist.Anbei ist dort auch eine Faktorisierung nach Pollard-Rho.
Das Prinzip ist einfach. statt GGT(TestZahl, kleineZahl)<>1 wird GGT(TestZahl,RiesigeTestzahl aus moeglichst vielen unterschiedlichen Faktoren) gebildet.
Wobei die Funktion, die die Testzahlen erzeugt, moeglichst GGT(f(x_i),f(x_j))=1 fuer i<>j liefern sollte.
Das heisst es werden immer neue Faktoren getestet.

Auch fuer die Bestimmung von Pi , kann man diese Bibliothek nutzen.

Gruss Horst

F34r0fTh3D4rk · Verfasst: Fr 21.04.06 16:04

Es gibt doch noch dieses Gittersieb-Verfahren, welches als eines der schnellsten gilt, ich hab nichts genaueres dazu gefunden, vielleicht hat jemand mal den algo (von mir aus auch als formulierung) ? [Zahlkörpersieb]

die version die ich momentan in meiner unit verwende ist diese hier:

markieren

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:

function p_primfaktor(a: int64): TPrimfaktor;
var
  Teiler: Int64;
begin
  if a < 0 then
  begin
    SetLength(result, 1);
    result[0] := -1;
    a := -a;
  end;
  while (a and 1) = 0 do
  begin
     setlength(result, length(result) + 1);
     Result[high(Result)] := 2;
     a := a shr 1;
  end;
  if a >= 2 then
  begin
    Teiler := 3;
    while a mod Teiler = 0 do
    begin  
      setlength(result, length(result) + 1);  
      Result[high(Result)] := Teiler;
      a := a div Teiler;
    end;
    Inc(Teiler, 2);
    while Teiler <= Trunc(SQRT(1.0 * a)) do
    begin
      while a mod Teiler = 0 do
      begin  
        setlength(result, length(result) + 1);  
        Result[high(Result)] := Teiler;
        a := a div Teiler;
      end;
      Inc(Teiler, 2);
      while a mod Teiler = 0 do
      begin  
        setlength(result, length(result) + 1);  
        Result[high(Result)] := Teiler;
        a := a div Teiler;  
      end;
      Inc(Teiler, 4);
    end;
  end;  
  if a <> 1 then
  begin
    setlength(result, length(result) + 1);
    Result[high(Result)] := a;
  end;
end;

wobei ich mir nicht sicher bin, ob Inc(Teiler, 4); richtig ist, es ist imgrunde nur die array variante ohne die variable root

Horst_H · Verfasst: Fr 21.04.06 16:37

Hallo,

Siehe Seite 2 diese Threads ab :
Verfasst am: Mo 21.11.05 12:15
Es geht nur um wheel-sieving, bei Dir nur ohne die Vielfachen von 2,3.
Ich glaube nicht das Du damit unter 6,6 Sekunden kommst.(1.8 Ghz Ahtlon64 , Duron 1800 ist 10% schneller, da schneller beim dividieren)

Gruss Horst

F34r0fTh3D4rk · Verfasst: Fr 21.04.06 17:27

ist wheel sieving das gleiche ?

Horst_H · Verfasst: Fr 21.04.06 18:02

Hallo,

im Prinzip ja.
Es ist einfach das Ausstreichen mit zu einer Reihe von Primzahlen (2,3,5,7,11,13...)teilerfremden Zahlen.
Also erst mit den kleinen Primzahlen probedividieren und anschliessend nur mit Zahlen die zu diesen teilerfremd sind.
Dabei ergibt sich ein regelmaessiges Muster von Differenzen.
Dein Beispiel benutzt eben nur zu 2,3 teilerfremde Zahlen.mit dem Muster +2,+4 (2+4=6=2*3)

Gruss Horst

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme