befreundete Zahlen

Fiete · Beiträge: 617 Erhaltene Danke: 364 W7 Delphi 6 pro

Das Programm berechnet in einem Intervall befreundete Zahlen.
user defined image

Die Teilersumme wird so berechnet:

Delphi-Quelltext

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

procedure TFreundschaft.TeilerSumme(N,Grenze:DWord;var Summe:DWord);
  var L:DWord;
  begin
   Summe:=1;
   for L:=2 to Grenze do
    if N mod L=0 then inc(Summe,L+N div L);
  end;

Die Prozedur hab ich in Assembler implementiert:

markieren

Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:

procedure TFreundschaft.TeilerSumme(N,Grenze:DWord;var Summe:DWord);
     asm
       PUSHAD         // alle RegisterInhalte auf den Stapel
       MOV ESI,Grenze // Wert von Grenze
       MOV EBX,N      // Wert von N
       MOV EDI,1      // Summe:=1
       MOV ECX,2      // Divisor, Startwert L:=2
@FORL: XOR EDX,EDX    // EDX:=0
       MOV EAX,EBX    // Dividend
       DIV ECX        // EDX:=N mod L UND EAX:=N div L !!!
       TEST EDX,EDX   // EDX=0?
       JNZ @RAUS
       ADD EDI,EAX    // inc(Summe,L+N div L);
       ADD EDI,ECX
@Raus: INC ECX
       CMP ECX,ESI    // L<=Grenze
       JBE @FORL
       MOV ESI,Summe  // 32-BitZeiger auf Summe
       MOV [ESI],EDI
       POPAD          // Werte vom Stapel holen
     end;

Viel Spaß beim Testen
Gruß Fiete

Einloggen, um Attachments anzusehen!

_________________
Fietes Gesetz: use your brain (THINK)

Mathematiker · Verfasst: Di 05.08.14 19:24

Hallo Fiete,
wie immer eine interessante Sache.
Allerdings ist die Berechnung der Teilersummen schneller möglich. Nach der Theorie ist

markieren

Quelltext

1:	Teilersumme s = (p[1]^(a[1]+1) -1) / (p[1] -1) * ... * (p[n]^(a[n]+1) -1) / (p[n] -1)

wobei die p[i] die Primfaktoren der Zahl sind und die a[i] die Häufigkeit des Auftretens.
Gerade für größere Zahlen spart man mit diesem Algorithmus eine Menge Rechenzeit.

Für alle befreundeten Zahlen bis 1 Million brauche ich damit 1,6 s; mit Deinem Programm sind es 4,6 s.
Ich hänge ein Beispiel für die Berechnung der Teilersummen mit Exe und Quelltext an. Vielleicht nützt es was.

Beste Grüße
Mathematiker

Einloggen, um Attachments anzusehen!

Fiete · Beiträge: 617 Erhaltene Danke: 364 W7 Delphi 6 pro

Moin Mathematiker,
habe Deinen Mathetipp eingebaut, brauche jetzt 1,57s für alle befreundeten Zahlen bis 999999.
Algorithmen sind doch besser als bloßes Rechnen :wink:

Gruß Fiete

Einloggen, um Attachments anzusehen!

_________________
Fietes Gesetz: use your brain (THINK)

Horst_H · Verfasst: Sa 09.08.14 12:23

Hallo,

eine Winzigkeit bringt noch was.Statt 1,7 dann 1,3 Sekunden.
Der Test auf befreundet nur, wenn die Teilersumme zuvor größer als die Zahl war.
Das ist nur bei ~25% der Zahlen.

markieren

Delphi-Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:

procedure TFreundschaft.ButtonRechneClick(Sender: TObject);
....
  for K := von to bis do
  begin
    TeilerSumme(K, Summe);
    //Nur bei ausreichendem Verdacht testen
    if K < Summe then
    Begin
      TeilerSumme(Summe, Freund);
      if Freund = K then
      begin
        Inc(Nr);
        Ausgabe.Lines.Add(Format('%4d', [Nr]) + ') ' + IntToStr(K) + ' - ' + IntToStr(Summe));
      end;
    end;
  end;

Gruß Horst

Fiete · Beiträge: 617 Erhaltene Danke: 364 W7 Delphi 6 pro

Moin Horst_H,
es gibt nichts was man nicht verbessern könnte,
hätte eigentlich selbst drauf kommen müssen(betriebsblind).
Gruß Fiete

_________________
Fietes Gesetz: use your brain (THINK)

Mathematiker · Verfasst: So 10.08.14 19:53

Hallo,
durch Pedersen wurden im Dezember 2000 über 733950 befreundete Zahlenpaare veröffentlicht. Bis 10^13 ist die Liste vollständig. Alle Paare waren entweder gerade-gerade oder ungerade-ungerade. Allerdings gibt es bis heute keinen Beweis, dass es keine gerade-ungerade Paare gibt; wenn es auch vermutet wird.

Vielleicht könnte man durch einen Test auf Parität noch etwas Zeit gewinnen. Mir bringt es nichts, da ich auch Paare suche, deren Teilersummen um einen gewissen Wert voneinander abweichen.
Außerdem könnte man auch die (wenigen) vollkommenen Zahlen im Suchintervall anzeigen. Immerhin sind diese zu sich "selbst befreundet".

Und dann hätte ich noch eine Idee:
Ketten sozialer Zahlen, Gesellige Zahlen

Ist a[2] die Teilersumme von a[1], a[3] die von a[2], …, und a[1] die von a[r], so bilden diese r Zahlen eine r-gliedrige Kette von sozialen Zahlen. (englisch "sociable numbers")
1969 entdeckte Henri Cohen sieben Ketten der Ordnung 4 und später fand Steve Root sechs weitere dieser Ketten. Insgesamt kennt man heute nach Moews 95 Ketten geselliger Zahlen: 88 der Länge 4, 1 der Länge 5, 2 der Länge 6, 2 der Länge 8, 1 der Länge 9, 1 der Länge 28.
Die 5gliedrige Kette (Poulet, 1918) ist: 12496, 14288, 15472, 14536, 14264, 12496 …
die 28gliedrige Kette: 14316, 19116, 31704, 47616, 83328, 177792, 295488, 629072, 589786, 294896, 358336, 418904, 366556,
274924, 275444, 243760, 376736, 381028, 285778, 152990, 122410, 97946, 48976, 46946, 22976, 22744,
19916, 17716, 14316, …

Die Berechnung derartiger Ketten finde ich sehr reizvoll, gestehe aber, dass ich mich da noch nicht herangewagt habe. Mir ist noch nicht klar, wie dies effektiv gemacht werden kann.
Vielleicht will sich ja einer von Euch versuchen.

Beste Grüße
Mathematiker

Horst_H · Verfasst: Mo 11.08.14 10:47

Hallo,

Wie soll man eine 51 gliedrige Kette finden, wenn man bei 50 aufhört zu testen :-)

Heutzutage könnte man ja zu extrem vielen Zahlen die Teilersumme bestimmen und im Speicher halten. Bei 2 GB könnten es schon 5e8 Zahlen sein.
Auf der erste Anschauen, scheinen die Zahlen sich in der gleichen Größenordnung zu bewegen.
Dann könnte man von Zahl zu Zahl springen in eine Liste eintragen und markieren / negativ machen, bei wieder auftreffen in der Liste zurückgehen.

Dann braucht man ein Verfahren extrem schnell die Teilersummen zu bestimmen 1 Sekunde für die ersten 1e6 mittels Probedivision dauert dann wohl zu lange.Man kann aber Dank user profile icon

Mathematiker Hinweis einfach Primzahlzerlegungen vorgeben und schnell berechnen, da man rekursiv wohl das Teilersummenprodukt kennt und nicht alles neu rechnen muss.
Also erst mal nur Potenzen Primzahlen 2..2^n, 3..3^m,.. Primzahl <= sqrt(Größe Zahl) ^2 (es bleiben große Primzahlen nicht berücksichtigt)
Dann Primzahl Paare ( a^n*b^m ), Drillinge.. Zehnlinge (29# = 6.469..e9 ) {Primzahlen a<b<c..}
user profile icon

Fiete benutzt ja eine Tabelle der Primzahlen bis sqrt(Maxlongint).
Da müßte man eine schöne Rekursion zu finden.

Vielleicht könnte man auch das Sieb des Erathotenes anwenden, aber der sagt ja nur: durch zwei Teilbar, aber nicht in welcher Potenz oder doch?
Dann könnte man einfach sieben und erst statt der Teilersumme den Anteil des Produktes speichern.

Zitat:

p[1]^(a[1]+1) -1) / (p[1] -1) *{ZinsesZins Formel (q^n-1)/(q-1) = Summe i=(0..n-1) [q^i]}

Die Potenz könnte man erhalten, wenn man die Vielfachheit als Zahl zur Basis der Primzahl merkt und einfach bei jedem Ausstreichen um 1 erhöht. Bei enem Überlauf in eine größere Z Potenz dort wäre das die aktuell maximale Potenz der Primzahl dort.

markieren

Quelltext

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

Zum Beispiel ausstreichen mit 2:
Vielfachheit 1  zur Basis 2= 1;  ( Ziffer n..1 )
Vielfachheit 2  zur Basis 2= 10; // Hier Überlauf auf Ziffer 2 == 2^2 = 4
Vielfachheit 3  zur Basis 2= 11; 
Vielfachheit 4  zur Basis 2= 100; // Hier höchster Überlauf auf Ziffer 3 == 2^3 = 8 
Vielfachheit 5  zur Basis 2= 101; 
Vielfachheit 6  zur Basis 2= 110; // Hier Überlauf auf Ziffer 2 == 2^2 = 4 ( 6*2 = 12 = 4*3 )
Vielfachheit 7  zur Basis 2= 111; 
Vielfachheit 8  zur Basis 2= 1000; // Hier höchster Überlauf auf Ziffer 4 == 2^4 = 16

Zu den Potenzen der Basis kann man ja einmalig p[1]^(a[1]+1) -1) / (p[1] -1) berechnen
Naja, schneller sieht das nicht aus, die Vielfachheit zur Basis p immer mit zu führen :-(

Gruß Horst
P.S. Geht es um:
de.wikipedia.org/wiki/Inhaltskette

Mathematiker · Verfasst: Mo 11.08.14 15:35

Hallo Horst,

Horst_H hat folgendes geschrieben :

P.S. Geht es um:
de.wikipedia.org/wiki/Inhaltskette

Streng genommen geht es um diese Ketten. Nur soll eben hier die Folge nicht auf 1 enden, sondern zyklisch wieder zur Startzahl zurückkehren.
Und wie Du schon geschrieben hast, liegt genau hier das Problem. Man kennt die Länge der Kette ja nicht. Ich werde auch einmal intensiver nachdenken.

Beste Grüße
Mathematiker

Horst_H · Verfasst: Mo 11.08.14 16:40

Hallo,

das Vorgehen habe ich doch kurz skizziert.
Ein riesiges Feld, im zu der betreffenden Zahl die Teilersumme steht.
Jede Zahl nehme ich dann mal als Startzahl

markieren

Quelltext

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

 
Wiederhole:
  Nehme Startzahl, packe diese in eine Ketten-Liste.
  Nehme mir dort ihre Teilersumme
    ist sie positiv
     trage diese nun negativ ein
    sonst 
      Abbruch und Kettenlänge bestimmen, (Ketten-Liste wieder abwickeln( wenn nicht negativ lassen-> war zuvor schon hier)
  Mache Teilersumme zu Startzahl
Mache bei nächster Zahl weiter

EDIT:
Hier mal ein Code-Schnipsel, der alle Teilersummen von 2..524*1000*1000 und anschliessend die befreundeten Zahlen, bestimmt
Man kann das Verfahren enorm beschleunigen, wenn man für kleine Primzahlen bis p^n < z.B. 10000 die Abfolge der Ergebnisse von IncP(p) abspeichert.{ 0..8191 würde nur einmal berechnet und nicht 65535 mal ), kommt vielleicht noch.

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:

program IncTest;
{$IFDEF FPC}
   {$MODE DELPHI}
   {$OPTIMIZATION ON}
   {$OPTIMIZATION PeepHole}
   {$OPTIMIZATION CSE}
   {$OPTIMIZATION ASMCSE}
   {$CODEALIGN loop=8,proc=32}
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}
uses
  sysutils;
const
  MAX = 524*1000*1000;
var
  Potenzen: array[0..31] of LongInt;
  Faktor : array[0..31] of LongWord;
  TeilerFeld : array[0..MAX] of LongWord;

procedure AusgabePot;
var
  i: NativeInt;
  s,a : string;
begin
  s:= '';
  For i := High(Potenzen) Downto Low(Potenzen) do
  begin
    str(Potenzen[i],a);
    s := s+a+':';
  end;
  Writeln(s);
end;

Procedure BestimmeTeilerFaktor(p:NativeInt);
var
  i: NativeInt;
  Pot,probe : Int64;
begin
  Pot := p;
  Probe := 1;
  For i := 0 to High(Faktor) do
  begin
    Probe := Probe+Pot;
    IF (POT > MAX){ OR (Probe > High(Faktor[0]))} then
      BREAK;
    Faktor[i] := Probe;
    Pot := Pot*p;
  end;
end;

function IncPStelle(p,s: NativeInt):NativeInt;
var
  i : NativeInt;
begin
  result := s;
  repeat
    i := Potenzen[result];
    Inc(i);
    IF i < p then
      BREAK
    else
    begin
      i := 0;
      Potenzen[result] := i;
      inc(result);
    end;
  until false;
  Potenzen[result] := i;
end;

function IncP(p: NativeInt):NativeInt;
var
  i : NativeInt;
begin
  result := 0;
  repeat
    i := Potenzen[result];
    Inc(i);
    IF i < p then
      BREAK
    else
    begin
      i := 0;
      Potenzen[result] := i;
      inc(result);
    end;
  until false;
  Potenzen[result] := i;
end;

procedure Check;
var
  i,s,n : NativeInt;
begin
  n := 0;
  For i := 2 to MAX do
  begin
    s := TeilerFeld[i]-i;
    // Damit s nie > MAX
    IF i > s then
    begin
    try
      IF TeilerFeld[s]= TeilerFeld[i] then
      begin
        inc(n);
        writeln(i:10,s:10);
      end;
    except
      writeln('Fehler ',i:10,s:10);
    end;
    end;
  end;
  writeln('Anzahl befreundet ',n, ' bis ',Max);
end;
type
  tpLW = ^LongWord;
var
  T1,T0: TDatetime;
  p,
  sum,p1 : NativeInt;
  pLW:tpLW;

begin
  T0:= time;
  TeilerFeld[0]:= 0;
  pLW := @TeilerFeld[1];
  For p := 1 to MAX do
  begin
    pLW^:= 1;
    inc(pLW);
  end;

  p := 2;
  while p*p <= MAX do
  begin
    fillchar(Potenzen,SizeOf(Potenzen),#0);
    BestimmeTeilerFaktor(p);
    sum := 0;
    repeat
      inc(sum,p);// Aktuelle Zahl
      IF sum > max then
        BREAK;
      TeilerFeld[sum] := TeilerFeld[sum]*Faktor[IncP(p)]
    until false ;
    repeat
      inc(p);
    until (TeilerFeld[p] = 1);
  end;

  //Ab hier nur einfache Potenzen
  //BestimmeFaktor kann wegfallen
  IF (TeilerFeld[p] <> 1) then
    repeat
      inc(p);
    until (TeilerFeld[p] = 1) OR (p > max);
  while p <= MAX do
  begin
    sum := 0;
    p1 := p+1;
    repeat
      inc(sum,p);
      IF sum > max then
        BREAK;
      TeilerFeld[sum] := TeilerFeld[sum]*p1
    until false ;
    repeat
      inc(p);
    until (p > max) OR (TeilerFeld[p] = 1);
  end;

  T1:= time;
  Check;
  writeln(FormatDateTime('HH:NN:SS.ZZZ' ,T1-T0));
end.

ergibt:

markieren

Quelltext

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:

sh-4.2# ./IncTest
       284       220
      1210      1184
      2924      2620
      5564      5020
 .....
 505757684 474179596
 511419856 491373104
 514780497 426191535
 514823985 475838415
 523679536 509379344
Anzahl befreundet 442 bis 524000000
00:00:14.722

15 Sekunden ist schon mal besser als ca 500 s

Gruß Horst

Fiete · Beiträge: 617 Erhaltene Danke: 364 W7 Delphi 6 pro

Hier mein erster Versuch:

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:

procedure TGesell.ButtonRechneClick(Sender: TObject);
  var Start,Nr,Summe,von,bis,K,KLaenge:DWord;
      TSek:Extended;
  begin
   if EditVon.Text='1' then EditVon.Text:='2';
   if EditVon.Text='' then EditVon.Text:='2';
   if EditBis.Text='' then EditBis.Text:='12496';
   von:=StrToInt(EditVon.Text);
   bis:=StrToInt(EditBis.Text);
   KLaenge:=StrToInt(EditKL.Text);
   LabelZeit.Caption:='';
   Ausgabe.Clear;
   Screen.Cursor:=crHourGlass;
   TSek:=TimeSekunden;
   for K:=von to bis do
    begin
     Start:=K;
     Nr:=0;
     repeat
      TeilerSumme(Start,Summe);
      Start:=Summe;
      inc(Nr);          // Min            Max
     until (Start=K) or (Start=1) or (Start>1000000000) or (Nr=KLaenge);
     if Start=K then Ausgabe.Lines.Add('Kettenlänge='+Format('%4d',[Nr])+' für '+IntToStr(K));
    end;
   TSek:=TimeSekunden-TSek;
   Screen.Cursor:=crDefault;
   LabelZeit.Caption:=Format('%0.2f',[TSek])+' sek.';
  end;

Das ganze Projekt ist im Anhang
Gruß Fiete

Einloggen, um Attachments anzusehen!

_________________
Fietes Gesetz: use your brain (THINK)

Horst_H · Verfasst: Di 12.08.14 20:22

Hallo,

das scheint ja schon gut zu funktionieren.
505757684 bis 505757684+1E6 =506757684

markieren

Quelltext

1:
2:
3:

Kettenlänge=   2 für 505757684
Kettenlänge=   4 für 506040584
Kettenlänge=   2 für 506468950

dauerte bei mir 56 Sekunden.
//506468950 findet mein Programm für befreundete Zahlen nicht, denn das Gegenstück ist 535495610 und ausserhalb von 524E6
Man könnte Die Kette sehr lang machen, wenn man Hashtabelle benutzt.Die Teilersumme scheint doch so zufällig , das MOD Primzahl ( > sqrt(MaxZahl) ) ein schnelles Auffinden bewerkstelligen könnte.
Also:

markieren

Delphi-Quelltext

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

tListe = record
               Zahl,
               Teilersumme : LongWord;
               FolgeWertMitGleichemHash : longWord;// 0 für keinen Nachfolger
          end;
Liste : array of tList;
HashListe : array of integer // size = PrimZahl;

Aus der Teilersumme einen Hashwert berechnen.
In der Hashliste steht die Position des Listenelementes.Sollte der Hashwert schon belegt sein, trägt man bei aktuellen Listenelement den dortigen Wert ein und ersetzt in der Hashliste den Wert mit der Position des neuen Listenelementes.
Dann kann sich die Liste entlanghangeln bis das Folgeelement eben 0 ist.

Gruß Horst
EDIT:
Teilersumme brauche ich ja nicht speichern, denn das nächste Listenelement wird ja genau dort weitermachen und dort steht die Teilersumme in Zahl.Zahl in Liste[a].Zahl und Teilersumme = Liste[a+1].zahl...

markieren

Delphi-Quelltext

1:
2:
3:
4:
5:

tListe = record
               Zahl : LongWord,
//               Teilersumme : LongWord;
               FolgeWertMitGleichemHash : longWord;// 0 für keinen Nachfolger
          end;

Bei jedem finden einer Kette muss Liste und HashListe wieder löschen

Jetzt sehe ich erst die Tabelle, der gefundenen:
djm.cc/sociable.txt da kann man mit simplen Brute Force nicht gegen anstinken ;-)

Horst_H · Verfasst: Do 14.08.14 00:37

Hallo,

hat jemand überhaupt eine Vorstellung, wie schnell die Suche ist
Aus amicable.homepage.dk/knwncx.htm

Zitat:

Update History:
13-Nov-2006
Added a new 6-cycle (13 digits) from Andre Needham. He has brought the exhaustive search up to 4×1012
08-Sep-2006
Added a new 6-cycle (13 digits) from Andre Needham. He has brought the exhaustive search up to 3×1012
15-May-2006
Added a new 6-cycle (13 digits) from Andre Needham. He has brought the exhaustive search up to 2×1012
[29-Dec-2005
Andre Needham has brought the exhaustive search up to 1012. No new cycles.]

Also 1x10E12 abklappern vom 8.9.2006 bis 13.11.2006 ~ 56 Tage
Pro Tag = 1,8E10= 207000 Bestimmung der Teilersumme von bis zu Zahlen/Sekunde.
Bis 2e6 gibt es ~2*78498 Primzahlen.Jetzt müßte man wissen wie stark im Durchschnitt die Zahlen zusammengesetzt sind ( jede 2.te durch 2 ;-)

)und wie lang die Kette ist, die gefunden werden soll. user profile icon

Fiete Verfahren ist sehr einfach und effektiv.
Man kann ja als Bedingung einführen, dass alle folgenden Teilersummen > Startzahl sind.
Man will schließlich mit dem kleinsten Startglied einer Kette anfangen und nicht mit nachfolgenden Kettengliedern die Kette nochmals entdecken.
Mein Siebverfahren scheitert in den Regionen, weil ich alle Primzahlen bis Max/2 kennen muss.
Aber muss ich das wirklich?
Wenn ich zu jeder Zahl neben der momentanen Teilersumme auch deren Primzahlfaktorprodukt merke ( wenn ich mit 8 streichen dann 8 wenn dann mit 3 dann eben 8*3 ) Kann ich, wenn ich mit der letzten Primzahl< sqrt(max) gesiebt habe, bei Zahlen, deren Primzahlfaktorprodukt < Aktuelle Zahl ist, dieses durch eine einzige Division bestimmen :-D

Statt nun bei 2e9 mit allen Primzahlen < 1e9( 50 Mio Zahlen) sieben/markieren zu müsse, brauche ich es nur für Primzahlen bis 44721 und dann nur noch den Test durchführen und passend dividieren.
Das ist sogar schneller, denn bei Abschnitten mit Startzahl <> 0 muss man ja per MOD Berechnung die Zahl an die richtige Position bringen.
Lange Rede keinen Sinn.
Ich kann die Bestimmung der ersten Teilersumme in der Kette wohl stark beschleunigen.Aber das kostet Platz Teilersumme und Produkt je 64 Bit dann noch 157000 Primzahlen mit der Darstellung der Startzahl des Abschnittes zur Basis der Primzahl (1e12 zur Basis 3 hat 25 Stellen, 23 nur 9 )
Also könnte ich abschnittsweise etwa 1e8 Zahlen ( = 1.6e9 Byte ) verarbeiten.

Gruß Horst

Horst_H · Verfasst: Do 14.08.14 13:04

Hallo,

ein kleines Update zu user profile icon

Fietes Version.
TeilerSumme ist jetzt fast doppelt so schnell, weil auch nur etwa halb so oft dividiert wird, weil MOD jetzt durch Multiplikation mit dem zuvor berechneten Quotienten ersetzt wird.
Die Bedingung (Start<k) in ButtonRechneClick spart ganz enorm.
Statt 56 Sekunden sind es jetzt nur 7 Sekunden für
505757684 bis 505757684+1E6 =506757684 { 2..999999 dauerte 5,8 Sekunden)}
Dafür wird auch nur eine befreundete Zahl bei 506468950 gefunden.
Das macht aber nur dann Sinn, wenn man große Bereiche hintereinander vollständig abgrast.

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:

procedure TGesell.TeilerSumme(N: DWord; var Summe: DWord);
var
  i, j, Nr, Primzahl: DWord;
  PrimPotSum, PrimPot, S: int64;
begin
  i := N;
  S := 1;
  Nr := 1;
  Primzahl := PrimListe[Nr];

  repeat
    PrimPotSum := 1;
    PrimPot := Primzahl;
    j := i div Primzahl;
    if (DWord(i) - DWord(j * Primzahl) = 0) then
    begin
      repeat
        i := j;
        j := j div Primzahl;
        PrimPotSum := PrimPotSum + PrimPot;
        PrimPot := PrimPot * Primzahl;
      until (i - j * Primzahl <> 0);
      S := S * PrimPotSum;
    end;
    Inc(Nr);
    Primzahl := PrimListe[Nr];
    if sqr(Primzahl) > i then
      Primzahl := i;
  until i = 1;

  Summe := S - N;
end;


procedure TGesell.ButtonRechneClick(Sender: TObject);
var
  Start, Nr, Summe, von, bis, K, KLaenge: DWord;
  TSek: extended;
begin
  if EditVon.Text = '1' then
    EditVon.Text := '2';
  if EditVon.Text = '' then
    EditVon.Text := '2';
  if EditBis.Text = '' then
    EditBis.Text := '999999';
  von := StrToInt(EditVon.Text);
  bis := StrToInt(EditBis.Text);
  KLaenge := StrToInt(EditKL.Text);
  LabelZeit.Caption := '';
  Ausgabe.Clear;
  Screen.Cursor := crHourGlass;

  TSek := TimeSekunden;
  for K := von to bis do
  begin
    Start := K;
    Nr := 0;
    repeat
      TeilerSumme(Start, Start);
      Inc(Nr);          // Min            Max
    until (Start<k) OR (Start = K) or (Start = 1) or (Start > 1000000000) or (Nr = KLaenge);
    if Start = K then
      Ausgabe.Lines.Add('Kettenlänge=' + Format('%4d', [Nr]) + ' für ' + IntToStr(K));
  end;
  TSek := TimeSekunden - TSek;

  Screen.Cursor := crDefault;
  LabelZeit.Caption := Format('%0.2f', [TSek]) + ' sek.';
end;

Vielleicht sollte man mal mit 64-Bit rechnen und dann bei 4e12 mal die Geschwindigkeit testen.

Gruß Horst
EDIT:
Was enorm was bringt, die maximale Abweichung vom Startwert begrenzen.
Statt 6 dann 1,75 Sekunden.

markieren

Delphi-Quelltext

1:
2:
3:
4:

// Statt 
or (Start > 1000000000) or
//mal 
 or (Start > 44*K) or

Ich weiß nicht, ob es allgemein gültig ist, wie ich oben bemerkte, sehen die Zahlen immer im der gleichen Größenordnung aus.
Nur der 28-er reißt weit aus mit 14316 .. 629072 ( 1:43,94 )
Von Wert zu Wert war die Abweichung maximal 2.13 (83328 -> 177792), aber das reduziert überhaupt keine Laufzeit.
Könnte es sein, das es dort eine Konstante gibt, scheinbar nicht...
1074427200 4368081600 = 4.065497969523
3161077920 12933375840 4.091444806903

Edit2:
Bei 4e12 dauerte es 00:00:30.871 Sekunden für 1E5 Zahlen abzuklappern.(EDIT: Ganz vergessen Linux 64Bit FPC 2.4.2 dividiert einfach schneller... )
Das ist viel zu langsam.Nur die erste Teilersumme zu berechnen dauerte schon 15 Sekunden. Das sind 6666/s.
Das ist ganz weit weg von 207000/s, da ist ein Faktor 31 drin.

EDIT3:
Mal für recht große Zahlen ein kleiner Test:

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:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:
228:
229:
230:
231:
232:
233:
234:
235:
236:
237:
238:
239:
240:
241:
242:

{$IFDEF FPC}
  {$MODE DELPHI}
  {$ASMMODE INTEL}
  {$OPTIMIZATION ON}
  {$OPTIMIZATION REGVAR}
  {$OPTIMIZATION PEEPHOLE}
  {$OPTIMIZATION CSE}
  {$OPTIMIZATION ASMCSE}
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}
//{$DEFINE USEASSEMBLER}
uses
  sysutils,crt;
type
  tPrimRec  =  LongWord;
  tPrimList = array of tPrimRec;
var
  PrimListe :tPrimList;
  cntDiv : NativeInt;
function GuessPrimCnt(Max:NativeUint):NativeUint;
begin
  result := Trunc(Max/(ln(Max)-1.08366))+200;
end;
procedure PrimlisteErstellen(Bis:LongWord;Var PL: TPrimList);
var
  sieve : array of byte;
  UpLim : LongWord;
  sieveprime,
  PrimePos,
  actPos,
  primeCnt : LongWord;
  procedure InsertPrime(inPrime:LongWord);
  begin
    PL[primeCnt] := inPrime;
    Inc(primeCnt);
  end;
begin
  Setlength(PL,0);
  IF (Bis > 200*1000*1000) or (Bis < 2) then
    exit;
  // nur ungerade 2n+1
  UpLim  := (Bis-1) SHR 1;
  Setlength(sieve,UpLim+1);  // automatisch 0 initialisiert
  Setlength(PL,GuessPrimCnt(Bis));
  primeCnt := 0;
  InsertPrime(2);
  sieveprime := 1;
  PrimePos   := 0; // == (sieveprime-1) shr 1
  while sqr(sieveprime) <= Bis do
  begin
    repeat
      inc(PrimePos);
      inc(SievePrime,2);
    until sieve[PrimePos] = 0;
    InsertPrime(sievePrime);
    //Pos of square
    actPos := (sqr(SievePrime) -1) shr 1;
    IF actPos>UpLim then
      BREAK;
    // streichen
    repeat
      sieve[actPos] := 1;
      inc(actPos,SievePrime);
    until actPos>UpLim;
  end;
  // die restlichen
  inc(PrimePos);
  inc(SievePrime,2);
  while PrimePos<=UpLim do
  begin
    IF sieve[PrimePos] = 0 then
      InsertPrime(sievePrime);
    inc(PrimePos);
    inc(SievePrime,2);
  end;
  Setlength(sieve,0);
  setlength(PL,primecnt);
end;
{$IFDEF USEASSEMBLER}
function UInt64divModCar(Teiler : Cardinal;var Dividend,Quotient:UInt64):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;
{$ENDIF}
procedure TeilerSumme(n: Uint64;var Summe:Uint64);
var
  Divi,Quot,
  S,PrimPotSum, PrimPot: UInt64;
  Primzahl,Rest : LongWord;
  PPrimZahl : ^tPrimRec;
begin
  Divi := N;
  Summe := Divi;
  S := 1;
  PPrimZahl := @PrimListe[0];
  PrimZahl := PPrimZahl^;
  repeat
    inc(PPrimZahl);
    PrimPotSum := 1;
    PrimPot := Primzahl;
    inc(cntDiv);
    {$IFDEF  USEASSEMBLER}
    Rest := UInt64divModCar(PrimZahl,Divi,Quot);
    {$ELSE}
    Quot := Divi div Primzahl;
    Rest := (Divi - Quot * Primzahl);
    {$ENDIF}
    if Rest = 0 then
    begin
      repeat
        Divi:= Quot;
        inc(cntDiv);
        {$IFDEF  USEASSEMBLER}
        Rest := UInt64divModCar(PrimZahl,Divi,Quot);
        {$ELSE}
        Quot := Divi div Primzahl;
        Rest := (Divi - Quot * Primzahl);
        {$ENDIF}
        PrimPotSum := PrimPotSum + PrimPot;
        PrimPot := PrimPot * Primzahl;
      until Rest<> 0;
      S := S * PrimPotSum;
      IF Divi = 1 then
        BREAK;
    end;
    Primzahl := pPrimZahl^;
    if sqr(Primzahl) > Divi then
    begin
      S := S * (Divi+1);
      Break;
    end;
  until Divi = 1;
  Summe := S-Summe;
end;

procedure ButtonRechneClick;//(Sender: TObject);
var
  Start, von, bis, K,Max: Uint64;
  KLaenge,Nr,cnt,cnt2,cntTS,dcntTs,LowCntTs,MaxCntTs : NativeInt;
  TSek: extended;
begin
{
  if EditVon.Text = '1' then
    EditVon.Text := '2';
  if EditVon.Text = '' then
    EditVon.Text := '2';
  if EditBis.Text = '' then
    EditBis.Text := '999999';
  von := StrToInt(EditVon.Text);
  bis := StrToInt(EditBis.Text);
  KLaenge := StrToInt(EditKL.Text);
  LabelZeit.Caption := '';
  Ausgabe.Clear;
  Screen.Cursor := crHourGlass;
}
  cnt := 0;
  cnt2 := 0;
  cntTS := 0;
  LowCntTs := 0;
  MaxCntTs := 0;
//  von := 1*5*1000*1000*1000*1000+1;
  von:=  1820741168916950;
  Bis := von+1000;

  writeln(von:20,bis:20);
  Max:= sqr(PrimListe[High(PrimListe)]);
  IF MAX < Bis then
  Begin
    writeln(' Geht nicht!');
    HALT;
  end;
  KLaenge := 100;
  TSek := Time;
  K :=von;
  while K<= BIS do
  begin
    Start := K;
    Max := 3*K shr 1;
    dcntTs := cntDiv;
    Nr := 0;
    repeat
      inc(cntTS);
      TeilerSumme(Start, Start);
      Inc(Nr);          // Min            Max
    until  (Start<=K)  or (Start = 1)  or (Start > Max) or (Nr = KLaenge);
    IF (Start< K)  AND (Nr=1) then
    begin
      inc(cnt);
      inc(LowCntTs,cntDiv-dcntTs);
    end
    else
      IF (Start> MAX) then
      begin
        inc(cnt2);
        inc(MaxCntTs,cntDiv-dcntTs);
      end;
    if (Nr >= 0) AND (K = Start) then
    begin
      writeln('Kettenlänge=' + Format('%4d', [Nr]) + ' für ' + IntToStr(K));
     // BREAK;
    end;
    inc(K);
  end;
  TSek := Time - TSek;
{
  Screen.Cursor := crDefault;
  LabelZeit.Caption := Format('%0.2f', [TSek]) + ' sek.';
 }
  writeln(FormatDateTime(' HH:NN:SS.ZZZ',TSek));
  Nr := Round(ln(cntDiv)/ln(10)+2);
  writeln('Zu klein       ', cnt:Nr,' Divs ',LowCntTs:Nr);
  writeln('Zu groß        ', cnt2:Nr,' Divs ',MaxCntTs:Nr);
  writeln('Aufrufe TS     ', cntTS:Nr);
  writeln('Anzahl  Div    ', cntDiv:Nr);
  writeln('Anzahl  Div/TS ', cntDiv/cntTS:Nr:2);
end;

begin
  PrimListeErstellen(200*1000*1000,PrimListe);
  ButtonRechneClick;
end.

WIeder linux 64-Bit FPC 2.6.4

markieren

Quelltext

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

sh-4.2# ./TeilersummenketteII
    1820741168916950    1820741168917950
Kettenlänge=   4 für 1820741168916950
 00:00:02.843
Zu klein              749 Divs  211012750
Zu groß               204 Divs   45869404
Aufrufe TS           1373
Anzahl  Div     278512353
Anzahl  Div/TS  202849.49

Man sieht ein Großteil der Divisionen bei der Bestimmung der ersten Teilersumme. Zweihunderttausend Divisionen...

Gruß Horst

Fiete · Beiträge: 617 Erhaltene Danke: 364 W7 Delphi 6 pro

Hier meine aktuelle Version, bin gestern erst aus dem Harz zurück.
user defined image

Habe

Horst_Hs Anregungen eingebaut.
Die Ketten werden explizit angezeigt, ebenso die befreundeten und vollkommenen Zahlen.
Die Abbruchbedingung Start>50*K lässt sich vielleicht noch besser einschränken.
Für einen geeigneten Maximalwert ist mir noch nichts geistreiches mathematisches eingefallen :oops:

Gruß Fiete

Edit1: Die Vorschläge von user profile icon

Horst_H sind umgesetzt (z.B.Start>50.0*K),
der Datentyp ist jetzt Int64 und nicht mehr Cardinal(DWord). Die Grenzen können bis zu 16-stellig sein.
Die Testläufe produzierten überwiegend 4-er Ketten, warum muss noch user profile icon

Mathematiker untersuchen :wink:

Kettenlänge=4 für 1820741168916950
2051026009623850
2310219512619350
2050935748436650
1820741168916950

Viel Spaß beim Testen
Fiete

Einloggen, um Attachments anzusehen!

_________________
Fietes Gesetz: use your brain (THINK)

Horst_H · Verfasst: Mo 25.08.14 08:02

Hallo,

das klappt ja schon ganz gut.
Weil mein neuer Rechner nicht mehr starten will, habe ich wieder den alten ausgemottet und siehe, dort ist Delphi7 installiert.
Gute Güte, der braucht dann direkt 2,6 Sekunden, obwohl nur 10% langsamerer Takt.
Intel-Haswell dividiert erheblich schneller ( 4 Bit/Takt statt 1 Bit/Takt+Overhead ). Zum Glück hat ja sqrt(1e6) nur 10 Bit.
Ich die Ausgabe bei Freunde in der Reihenfolge geändert, klein vor groß ;-)

Zeile 163:...IntToStr(Kette[2])+' - '+IntToStr(Kette[1])

Gruß Horst

Mathematiker · Verfasst: Di 26.08.14 09:37

Hallo,

Fiete hat folgendes geschrieben :

Hier meine aktuelle Version, bin gestern erst aus dem Harz zurück.
... Die Ketten werden explizit angezeigt, ebenso die befreundeten und vollkommenen Zahlen.

Das Programm funktioniert schon richtig gut, und vor allem ist es schnell!
Allerdings findet es bei der Startzahl 805984760 nicht die Kette der Länge 9:

markieren

Quelltext

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

1  805984760
2  1268997640
3  1803863720
4  2308845400
5  3059220620
6  3367978564
7  2525983930
8  2301481286
9  1611969514
10  805984760

Ich habe selbst schon gesucht, habe aber noch keine Ursache gefunden.

Beste Grüße
Mathematiker

Horst_H · Verfasst: Di 26.08.14 10:30

Hallo,

Ich rate mal, das es an der Bestimmung der Primzahlen liegen könnte.
Du hast ja zwischendurch Zahlen bis 3,367x10^9, da ist die Frage, ob die Obergrenze in Edit2, die für die maximale Suchprimzahl verwendet wird, ausreichend groß gewählt ist.
Innerhalb von Sieben , hilft vielleicht die Erhöhung von "bis" auf 44*Edit2 ( Das Maximum, was ich bei der 28er Kette gesehen habe ).

Gruß Horst

Fiete · Beiträge: 617 Erhaltene Danke: 364 W7 Delphi 6 pro

Moin Mathematiker,
der Haken ist die Abbruchbedingung Start>50*K, mit Start>4000000000 (Datentyp Cardinal)ergibt sich die 9-er Kette.
Der Maximalwert muss noch mathematisch gefunden werden!
Schöne Aufgabe für Dich :wink:

Gruß Fiete

_________________
Fietes Gesetz: use your brain (THINK)

Horst_H · Verfasst: Di 26.08.14 11:28

Hallo,

Der maximale Wert (3,3e9) ist doch nur knapp über dem Vierfachem ( 0.8e9 )des Startwertes.

Gruß Horst.

Horst_H · Verfasst: Di 26.08.14 12:54

Hallo,

MIst , zurück marsch marsch
Alles falsch verstanden 50*K funktioniert nicht weil das > Cardinal wird.
mit 50.0*K klappt es..
Dennoch habe ich sieben geändert. Zuvor wurden alle Primzahlen bis zur Obergrenze gesucht,damit findet man alle Primzahlen bis zum Quadrat dieser Obergrenze, also reichlich zuviele, was nicht Not tut und bei 1e9 auch etwas dauert...

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:

procedure TGesell.Sieben;
var
  an, Bis, k, z, Anzahl, i, Wurzel, Index, SiebMax: cardinal;
begin
{
   Bis := StrToInt(EditBis.Text);
   Primlistenlaenge abschaetzen  ohne dem fuer  805984760 Lazarus 1.0.14
   200 , so 16 Sekunden 
}
  // Maximale Obergrenze bei 28er Kette ist 43.9
  Bis := Round(sqrt(44.0 * (StrToInt(EditBis.Text))));
  //Primlistenlaenge abschaetzen
  setlength(Primliste,Trunc(Bis/(ln(Bis*1.0)-1.08366))+200);
  an := (Bis - 1) div 2;
  SiebMax := an div 32 + 1;
  SetLength(Sieb, SiebMax + 1);
  for i := 1 to SiebMax do
    Sieb[i] := $FFFFFFFF;
  Wurzel := trunc(sqrt(an / 2 + 0.25) - 0.5); // (2w+1)²=2an+1 !!! nur ungerade Zahlen
  for i := 1 to Wurzel do
  begin
    index := (i - 1) div 32 + 1;
    if TestBit(sieb[index], i mod 32) then
    begin
      z := 2 * i + 1;
      k := i * (1 + z);
      while k <= an do
      begin
        index := (k - 1) div 32 + 1;
        sieb[index] := ClrBit(sieb[index], k mod 32);
        Inc(k, z);
      end;
    end;
  end;
  Anzahl := 1;
  PrimListe[Anzahl] := 2;
  for i := 1 to an do
  begin
    index := (i - 1) div 32 + 1;
    if TestBit(sieb[index], i mod 32) then
    begin
      Inc(Anzahl);
      PrimListe[Anzahl] := 2 * i + 1;
    end;
  end;
  // Laenge jetzt passsend kuerzen
  SetLength(PrimListe, Anzahl + 1);
end;

Gruß Horst

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme