Superperfekte Zahlen

harryp · Beiträge: 52 Erhaltene Danke: 9 Win 8.1 Delphi 7, XE8

Hi,

ich bin gerade dabei mir (angefixt vom user profile icon

Mathematiker) ein paar zahlentheoretische Verfahren anzusehen. Dabei bin ich eben auf superperfekte Zahlen gestoßen (Zahlen, bei denen die Teilersumme der Teilersumme der Zahl gerade dem Doppelten der Zahl entspricht). Nach ein wenig Rumbasteln hab ich eine "naive" Variante der Ermittlung superperfekter Zahlen und eine etwas bessere Variante geschafft. Während die "naive" Variante knappe 5 Minuten zur Ermittlung aller superperfekten Zahlen bis zur 275.000 benötigt, sind es bei der etwas besseren noch knapp über 30 Sekunden.

Da ich gern noch etwas dazulernen möchte, meine Frage an euch: was kann ich an dieser Version noch effizienter gestalten oder gibt es evtl. eine generell andere Idee der Ermittlung, die effizienter ist als meine?

Anbei meine "etwas bessere" Version:

markieren

Delphi-Quelltext

1:
2:
3:

  for k := 1 to obergrenze do
    if Teilersumme2(Teilersumme2(k))=2*k then
      ListBox1.Items.Add(IntToStr(k));

und

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:

function TForm1.Teilersumme2(zahl: integer): integer;
 var
  summe, k: integer;
  faktor, vielfachheit: Array of integer;
begin
  // Summe ist bei der Berechnung ein Produkt -> beginnt bei 1
  summe := 1;

  // Primfaktorzerlegung von zahl ermitteln
  setLength(faktor,0);
  setLength(vielfachheit,0);
  k := 2;
  while zahl <> 1 do
   begin
     if zahl mod k = 0 then
      begin
        setLength(faktor,Length(faktor)+1);
        setLength(vielfachheit,Length(vielfachheit)+1);
        faktor[Length(faktor)-1] := k;
        vielfachheit[Length(vielfachheit)-1] := 1;
        zahl := zahl div k;
      end;
      while zahl mod k = 0 do
       begin
        vielfachheit[Length(vielfachheit)-1] := vielfachheit[Length(vielfachheit)-1]+1;
        zahl := zahl div k;
       end;
      k := k + 1;
   end;

  // Teilersumme ermitteln
  for k := 0 to length(faktor)-1 do
   begin
     summe := summe * (Potenz(faktor[k],vielfachheit[k]+1)-1) div (faktor[k]-1);
   end;

  // Rückgabe der Teilersumme
  Teilersumme2 := summe;
end;

Randbemerkung: Ich hab für round(Power(...)) einen Square-and-Multiply-Ansatz versucht, der scheitert aber schlagartig irgendwo zwischen 46.000 und 47.000.

edit: Square-and-Multiply-Problem lag am Wertebereich von Integer; mit Int64 gehts jetzt und hat mir nochmal 3 Sekunden gebracht.

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:

function TForm1.Potenz(basis, exponent: integer): int64;
 var
  b: String;
  k: integer;
  ergebnis: int64;
begin
  ergebnis := 1;

  // Binärdarstellung vom Exponent ermitteln
  b := '';
  while exponent > 0 do
   begin
     b := IntToStr(exponent mod 2)+ b;
     exponent := exponent div 2;
   end;

  // Square&Multiply
  for k := 1 to Length(b) do
   begin
     ergebnis := ergebnis * ergebnis;                   // Square
     if b[k] = '1' then ergebnis := ergebnis * basis;   // Multiply
   end;

  // Rückgabe der Potenz
  Potenz := ergebnis;
end;

Danke schon mal im Vorfeld für eure Ideen und eine schöne Weihnachtszeit noch (an die die noch feiern).

Mathematiker · Verfasst: Fr 26.12.14 00:04

Hallo,
zuerst zur Theorie: Eine Zahl n ist genau dann superperfekt, wenn n eine Zweierpotenz 2^k $2^k$ und 2^(k+1)-1 $2^(k+1)-1$ eine Mersennesche Primzahl ist. Der Beweis steht in www.digizeitschrifte...850199_0024/log7.pdf

Ungerade superperfekte Zahlen kennt man bis etwa 10^24 $10^24$ nicht. Bei dieser Größenordnung wird das Ermitteln einer superperfekten ungeraden Zahl mit einem "normalen" Programm aussichtslos. Es scheitert am schnellen Faktorisieren der Zahl.

Zum Programm: Das ständige Erhöhen der dynamischen Felder

markieren

Delphi-Quelltext

1:
2:

        setLength(faktor,Length(faktor)+1);
        setLength(vielfachheit,Length(vielfachheit)+1);

kostet sehr viel Zeit. Ich würde das Feld der Faktoren von Anfang an auf z.B. 20 Elemente setzen. Die kleinste Zahl mit 20 verschiedenen Primfaktoren 2*3*5*7*... $2*3*5*7*...$ dürfte die Int64-Grenze überschreiten. (Ich hab's nicht ausgerechnet :wink:

) Das Feld der Vielfachheit wird mit 64 Elementen groß genug sein.

Außerdem prüfst du, ob jede Zahl k ab 2 aufwärts Teiler sein kann. Ich würde nur die Primzahlen für k wählen. Ein schneller Primzahltest ist z.B.

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:

function IstPrime(N: Cardinal): Boolean;
const
  MaxPrime = 137;
  MaxTrial = MaxPrime * MaxPrime; // maximal trial div bounds
  MinPrime = 137;                 // trial div bounds if N >= MaxTrial, must be <= MaxPrime
  Primes: array[0..31] of Byte =
   (   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,
      61,  67,  71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137);

  InvPrimes: array[0..31] of Cardinal =  // InvPrimes[i] = Primes[i]^-1 mod 2^32
   ($AAAAAAAB,$CCCCCCCD,$B6DB6DB7,$BA2E8BA3,$C4EC4EC5,$F0F0F0F1,$286BCA1B,$E9BD37A7,
    $4F72C235,$BDEF7BDF,$914C1BAD,$C18F9C19,$2FA0BE83,$677D46CF,$8C13521D,$A08AD8F3,
    $C10C9715,$07A44C6B,$E327A977,$C7E3F1F9,$613716AF,$2B2E43DB,$FA3F47E9,$5F02A3A1,
    $7C32B16D,$D3431B57,$8D28AC43,$DA6C0965,$0FDBC091,$EFDFBF7F,$C9484E2B,$077975B9);
asm
       TEST   AL,1
       JZ     @@0                          // even ??
       CMP    EAX,7
       JA     @@1                          // > 7 ??
       DEC    EAX
       SETNZ  AL
       RET
@@0:   CMP    EAX,2                        // even numbers
       SETZ   AL
       RET
@@1:   PUSH   EBP                          // do trial divsion to small primes
       PUSH   EBX
       PUSH   ESI
       PUSH   EDI
       MOV    EBX,EAX
       CMP    EAX,MaxTrial
       MOV    EBP,MinPrime
       JAE    @@2
       PUSH   EAX
       FILD   DWord Ptr [ESP]
       FSQRT
       FISTP  DWord Ptr [ESP]
       POP    EBP                          // EBP = Sqrt(N)
@@2:   MOV    EDI,Offset Primes
       MOV    ESI,Offset InvPrimes
       XOR    ECX,ECX
@@3:   MOVZX  EDX,Byte Ptr [EDI + ECX]     // take care, InvPrimes[] MUST
                                          // be after Primes[] declared
       MOV    EAX,EBX
       CMP    EDX,EBP
       JA     @@5
       IMUL   EAX,[ESI + ECX * 4]
       INC    ECX
       MUL    EDX
       JC     @@3
       TEST   EDX,EDX
@@4:   POP    EDI
       POP    ESI
       POP    EBX
       POP    EBP
       SETNZ  AL
       RET
@@5:   CMP    EBX,MaxTrial                 // N <= MaxPrime^2 ??
       MOV    EAX,EBX
       JBE    @@4
       IMUL   EAX,EBX                      // compute domain param U = -N^-1 mod 2^32
       SUB    EAX,2                        // Lookuptable can reduce from 72 donwto 32 cycles
       IMUL   EAX,EBX
       MOV    EDX,EAX
       IMUL   EAX,EBX
       ADD    EAX,2
       IMUL   EAX,EDX
       MOV    EDX,EAX
       IMUL   EAX,EBX
       ADD    EAX,2
       IMUL   EAX,EDX
       MOV    EBP,EAX
       IMUL   EBP,EBX
       ADD    EBP,2
       IMUL   EBP,EAX                      // U = -N^-1 mod^2^32 in EBP
       MOV    EDI,EBX
       MOV    EAX,EBX
       DEC    EDI                          // N -1
       NEG    EAX
       BSF    ECX,EDI
       MUL    EAX
       PUSH   ECX                          // bits remain         [ESP + 20]
       MOV    ESI,EAX
       BSR    ECX,EDI
       MOV    EAX,EDX
       XOR    EDX,EDX
       NEG    ECX
       DIV    EBX                          // div, can't be removed
       MOV    EAX,ESI
       ADD    ECX,32
       DIV    EBX                          // div
       SHL    EDI,CL
       MOV    EAX,EDX
       IMUL   EAX,EBP
       MOV    ESI,EDX                      // C = -N^2 mod N, to fast convert into
       MUL    EBX                          // montgomery domain
       PUSH   ESI                          // C                    [ESP + 16]
       ADD    EAX,ESI
       ADC    EDX,0
       PUSH   EDX                          // +1 in montgomery     [ESP + 12]
       PUSH   EDI                          // bit mask exponent    [ESP +  8]
       NEG    EDX
       ADD    EDX,EBX
       CMP    EBX,$08A8D7F                 // N < $08A8D7F, do SPP to bases 31,73
       PUSH   EDX                          // -1 in montgomery     [ESP +  4]
       JAE    @@6
       MOV    EAX,31
       CALL   @@9
       MOV    EAX,73
       PUSH   Offset @@7
       JMP    @@9
@@6:   MOV    EAX,2                        // do SPP to bases 2,7,61
       CALL   @@9
       MOV    EAX,7
       CALL   @@9
       MOV    EAX,61
       CALL   @@9
@@7:   INC    EAX
@@8:   LEA    ESP,[ESP + 4 * 5]            // don't change flags !!
       JMP    @@4
@@9:   MUL    DWord Ptr [ESP + 16]         // convert base in montgomery
       MOV    EDI,EAX                      // Base' = Base * C mod N
       IMUL   EAX,EBP                      // montgomery REDC
       MOV    ESI,EDX
       MUL    EBX
       ADD    EAX,EDI
       ADC    EDX,ESI
       MOV    ECX,[ESP + 8]                // bit mask of exponent N -1
       MOV    EAX,EDX
       PUSH   EDX
@@A:   MUL    EAX                          // X = X^2 mod N
       MOV    EDI,EAX
       IMUL   EAX,EBP
       MOV    ESI,EDX
       MUL    EBX
       ADD    EAX,EDI
       ADC    EDX,ESI
       JNC    @@B
       SUB    EDX,EBX
@@B:   ADD    ECX,ECX
       MOV    EAX,EDX
       JNC    @@D
       MUL    DWord Ptr [ESP]              // X = X * Base mod N
       MOV    EDI,EAX
       IMUL   EAX,EBP
       MOV    ESI,EDX
       MUL    EBX
       ADD    EAX,EDI
       ADC    EDX,ESI
       JNC    @@C
       SUB    EDX,EBX
@@C:   TEST   ECX,ECX
       MOV    EAX,EDX
@@D:   JNZ    @@A
       CMP    EAX,EBX
       JB     @@E
       SUB    EAX,EBX
@@E:   CMP    EAX,[ESP + 16]  
       MOV    ECX,[ESP + 24]           
       JE     @@J
@@F:   CMP    EAX,[ESP + 8]           
       JE     @@J
       DEC    ECX
       JNG    @@I
       MUL    EAX
       MOV    EDI,EAX
       IMUL   EAX,EBP
       MOV    ESI,EDX
       MUL    EBX
       ADD    EAX,EDI
       ADC    EDX,ESI
       JC     @@G
       CMP    EDX,EBX
       JB     @@H
@@G:   SUB    EDX,EBX
@@H:   CMP    EDX,[ESP + 16]        
       MOV    EAX,EDX
       JNE    @@F
@@I:   ADD    ESP,8
       XOR    EAX,EAX
       JMP    @@8
@@J:   POP    EDX
end;

Diese Routine ist nicht von mir. Ich muss noch einmal nachsehen, wer der Urheber war. :nixweiss:

Ich glaube user profile icon

Gammatester.

Beste Grüße und Gratulation zum Gewinn :zustimm:

Mathematiker

_________________
Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein

harryp · Beiträge: 52 Erhaltene Danke: 9 Win 8.1 Delphi 7, XE8

Lieber Mathematiker,

erstmal danke für die schnelle Hilfe.

Die Einschränkung der Feldgrößen auf 16 (das Produkt der ersten 16 Primzahlen beträgt etwa 3,2*10^19 und ist damit größer als der Int64-Bereich) und damit einhergehende Änderung auf ein statisches Feld hat mir nochmal etwas gebracht.

Ein Problem hab ich noch mit der istPrime-Funktion. Die liefert bei mir (von dem Quelltext hier einfach komplett kopiert) sowohl bei 3 als auch bei 4, genauer sogar generell, den Rückgabewert false. Da ich mich mit Assembler nicht auskenne - kann das am System liegen?

Die Idee mit der Mersenneschen Primzahl werd ich nachher mal noch testen. Schon erstaunlich, auf wie viele zahlentheoretische Inhalte man bei dem scheinbar kleinen Thema "Superperfekte Zahlen" stößt (Teilersumme, Primfaktorzerlegung, Mersenne-Primzahl, schnelles Potenzieren, Primzahltest, ...).

Mathematiker · Verfasst: Fr 26.12.14 16:48

Hallo,

harryp hat folgendes geschrieben :

Ein Problem hab ich noch mit der istPrime-Funktion. Die liefert bei mir (von dem Quelltext hier einfach komplett kopiert) sowohl bei 3 als auch bei 4, genauer sogar generell, den Rückgabewert false. Da ich mich mit Assembler nicht auskenne - kann das am System liegen?

Ich hänge einmal ein Minitestprogramm für die istprime-Routine an. Bei mir funktioniert es (Win 8.1).
Wenn es nicht geht, tut mir leid, würde ich eine andere Primzahltestroutine zusammenbasteln. Ich denke, dass das Ausschließen der zusammengesetzten Zahlen doch etwas beschleunigt.

Der Profi für das Beschleunigen von derartigen Berechnungen ist aber user profile icon

Horst_H.

Beste Grüße
Mathematiker

Einloggen, um Attachments anzusehen!

_________________
Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein

harryp · Beiträge: 52 Erhaltene Danke: 9 Win 8.1 Delphi 7, XE8

Vielen Dank für den Test. Das istPrime-Problem hat sich mehr oder weniger geklärt (auch wenn ich noch nicht verstehe, warum es als Methode von TForm1 nicht funktioniert, wohl aber wenn man die Methode frei in die Unit reinschreibt).

Der Tipp mit den Mersenne-Primzahlen war jedoch Gold wert. Anstelle alle natürlichen Zahlen durchzutesten werden jetzt Mersenne-Primzahlen erzeugt und getestet. Das vereinfacht das Problem natürlich vollständig und so kommen alle mit Int64 erreichbaren jetzt innerhalb von einer Sekunde.

markieren

Delphi-Quelltext

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

  
  // Ermittlung über Mersenne-Primzahlen
  r := 1;
  repeat
    k := Potenz(2,r+1)-1;
    if istPrime(k) then ListBox1.Items.Add(IntToStr(Potenz(2,r)));
    r := r + 1;
  until Potenz(2,r)>obergrenze;

Danke

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme