Entwickler-Ecke
Open Source Projekte - befreundete Zahlen
Fiete - Di 05.08.14 16:57
Titel: befreundete Zahlen
Das Programm berechnet in einem Intervall befreundete Zahlen.
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:
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
Mathematiker - 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
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
Fiete - Fr 08.08.14 17:08
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
Horst_H - 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.
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);
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 - So 10.08.14 10:16
Moin Horst_H,
es gibt nichts was man nicht verbessern könnte,
hätte eigentlich selbst drauf kommen müssen(betriebsblind).
Gruß Fiete
Mathematiker - 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 - 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
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..}
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.
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:
http://de.wikipedia.org/wiki/Inhaltskette
Mathematiker - Mo 11.08.14 15:35
Hallo Horst,
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 - 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
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.
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) 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;
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); IF sum > max then
BREAK;
TeilerFeld[sum] := TeilerFeld[sum]*Faktor[IncP(p)]
until false ;
repeat
inc(p);
until (TeilerFeld[p] = 1);
end;
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:
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 - Di 12.08.14 19:15
Hier mein erster Versuch:
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); 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
Horst_H - Di 12.08.14 20:22
Hallo,
das scheint ja schon gut zu funktionieren.
505757684 bis 505757684+1E6 =506757684
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:
Delphi-Quelltext
1:
2:
3:
4:
5:
6:
7:
| tListe = record
Zahl,
Teilersumme : LongWord;
FolgeWertMitGleichemHash : longWord; end;
Liste : array of tList;
HashListe : array of integer |
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...
Delphi-Quelltext
1:
2:
3:
4:
5:
| tListe = record
Zahl : LongWord,
FolgeWertMitGleichemHash : longWord; end; |
Bei jedem finden einer Kette muss Liste und HashListe wieder löschen
Jetzt sehe ich erst die Tabelle, der gefundenen:
http://djm.cc/sociable.txt da kann man mit simplen Brute Force nicht gegen anstinken ;-)
Horst_H - Do 14.08.14 00:37
Hallo,
hat jemand überhaupt eine Vorstellung, wie schnell die Suche ist
Aus
http://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.
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 - Do 14.08.14 13:04
Hallo,
ein kleines Update zu
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.
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); 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.
Delphi-Quelltext
1:
2:
3:
4:
| or (Start > 1000000000) or
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:
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}
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;
UpLim := (Bis-1) SHR 1;
Setlength(sieve,UpLim+1); Setlength(PL,GuessPrimCnt(Bis));
primeCnt := 0;
InsertPrime(2);
sieveprime := 1;
PrimePos := 0; while sqr(sieveprime) <= Bis do
begin
repeat
inc(PrimePos);
inc(SievePrime,2);
until sieve[PrimePos] = 0;
InsertPrime(sievePrime);
actPos := (sqr(SievePrime) -1) shr 1;
IF actPos>UpLim then
BREAK;
repeat
sieve[actPos] := 1;
inc(actPos,SievePrime);
until actPos>UpLim;
end;
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;
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;var
Start, von, bis, K,Max: Uint64;
KLaenge,Nr,cnt,cnt2,cntTS,dcntTs,LowCntTs,MaxCntTs : NativeInt;
TSek: extended;
begin
cnt := 0;
cnt2 := 0;
cntTS := 0;
LowCntTs := 0;
MaxCntTs := 0;
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); 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));
end;
inc(K);
end;
TSek := Time - TSek;
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
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 - So 24.08.14 15:36
Hier meine aktuelle Version, bin gestern erst aus dem Harz zurück.
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
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
Mathematiker untersuchen :wink:
Kettenlänge=4 für 1820741168916950
2051026009623850
2310219512619350
2050935748436650
1820741168916950
Viel Spaß beim Testen
Fiete
Horst_H - 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 - 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:
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 - 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 - Di 26.08.14 11:21
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
Horst_H - 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 - 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...
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 := Round(sqrt(44.0 * (StrToInt(EditBis.Text))));
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); 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;
SetLength(PrimListe, Anzahl + 1);
end; |
Gruß Horst
Fiete - Sa 30.08.14 13:21
Moin,
die Version vom Sonntag liegt jetzt in der finalen Version vor:
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:
| unit Kette;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls;
type
TGesell = class(TForm)
ButtonRechne: TButton;
Ausgabe: TMemo;
LabelRZ: TLabel;
LabelZeit: TLabel;
LabelVon: TLabel;
LabelBis: TLabel;
EditVon: TEdit;
EditBis: TEdit;
LabelKL: TLabel;
EditKL: TEdit;
AusgabeF: TMemo;
AusgabeV: TMemo;
LabelK: TLabel;
LabelF: TLabel;
LabelV: TLabel;
procedure ButtonRechneClick(Sender: TObject);
private
PrimListe:array of Int64;
Sieb:array of Cardinal;
procedure Sieben;
procedure TeilerSumme(N:Int64;var Summe:Int64);
function TestBit(Zahl:Cardinal;BitNr:Byte):Boolean;
function SetBit(Zahl:Cardinal;BitNr:Byte):Cardinal;
function ClrBit(Zahl:Cardinal;BitNr:Byte):Cardinal;
function TimeSekunden:Extended;
public
end;
var
Gesell: TGesell;
implementation
{$R *.dfm}
{$R-,Q-}
function TGesell.TimeSekunden:Extended;
var H, M, S, MS : Word;
begin
DecodeTime(Now,H,M,S,MS);
TimeSekunden:=3600.0*H+60.0*M+S+MS/1000
end;
function TGesell.TestBit(Zahl:Cardinal;BitNr:Byte):Boolean;
begin TestBit:=(((Zahl shr BitNr) and 1)=1) end;
function TGesell.SetBit(Zahl:Cardinal;BitNr:Byte):Cardinal;
begin SetBit:=Zahl or (1 shl BitNr) end;
function TGesell.ClrBit(Zahl:Cardinal;BitNr:Byte):Cardinal;
begin ClrBit:=Zahl and not(1 shl BitNr) end;
procedure TGesell.Sieben;
var an,Bis,k,z,Anzahl,i,Wurzel,Index,SiebMax:Cardinal;
begin
Bis:=Round(sqrt(50.0*(StrToInt64(EditBis.Text))));
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); 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;
SetLength(PrimListe,Anzahl+1);
end;
procedure TGesell.TeilerSumme(N:Int64;var Summe:Int64);
var Quotient,Dividend,PrimPotSum,PrimPot,S,Primzahl:Int64;
Nr:Cardinal;
begin
Dividend:=N;S:=1;Nr:=1;
Primzahl:=PrimListe[Nr];
repeat
PrimPotSum:=1;
PrimPot:=Primzahl;
Quotient:=Dividend div Primzahl;
if Dividend=Quotient*Primzahl then
begin
repeat
Dividend:=Quotient;
Quotient:=Quotient div Primzahl;
PrimPotSum:=PrimPotSum+PrimPot;
PrimPot:=PrimPot*Primzahl;
until Dividend<>Quotient*Primzahl;
S:=S*PrimPotSum;
if Dividend=1 then
begin
Summe:=S-N;
exit
end;
end;
inc(Nr);
Primzahl:=PrimListe[Nr];
if Primzahl*Primzahl>Dividend then
begin
S:=S*(Dividend+1);
Summe:=S-N;
exit
end;
until Dividend=1;
Summe:=S-N
end;
procedure TGesell.ButtonRechneClick(Sender: TObject);
var Start,Summe,von,bis,K:Int64;
Nr,L,KLaenge:Cardinal;
TSek:Extended;
Kette:Array of Int64;
begin
if EditVon.Text='1' then EditVon.Text:='2';
if EditVon.Text='' then EditVon.Text:='2';
if EditBis.Text='' then EditBis.Text:='10000000';
von:=StrToInt64(EditVon.Text);
bis:=StrToInt64(EditBis.Text);
KLaenge:=StrToInt(EditKL.Text);
SetLength(Kette,KLaenge+1);
LabelZeit.Caption:='';
Ausgabe.Clear;AusgabeF.Clear;AusgabeV.Clear;
Screen.Cursor:=crHourGlass;
Sieben;
TSek:=TimeSekunden;
K:=von;
while K<=bis do
begin
Start:=K;
Nr:=0;
repeat
TeilerSumme(Start,Summe);
Start:=Summe;
inc(Nr);
Kette[Nr]:=Summe;
until (Start<=K) or (Start=1) or (Start>50.0*K) or (Nr=KLaenge);
if Start=K then
begin
if Nr=1 then AusgabeV.Lines.Add(IntToStr(K))
else if Nr=2 then AusgabeF.Lines.Add(IntToStr(Kette[2])+' - '+IntToStr(Kette[1]))
else
begin
Ausgabe.Lines.Add('Kettenlänge='+IntToStr(Nr)+' für '+IntToStr(K));
for L:=1 to Nr do Ausgabe.Lines.Add(IntToStr(Kette[L]));
end
end;
inc(K)
end;
TSek:=TimeSekunden-TSek;
Screen.Cursor:=crDefault;
LabelZeit.Caption:=Format('%0.2f',[TSek])+' sek.';
end;
end. |
Viel Spaß beim Testen
Fiete
Horst_H - Sa 30.08.14 15:47
Hallo,
eine kleine Änderung, zur Ausgabe der Anzahl der Kette gewisser Länge.
Weil K>2 reicht Start> K, dann ist K >1.
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:
| implementation
const
cFaktor = 50;
...
procedure TGesell.Sieben;
var an,Bis,k,z,Anzahl,i,Wurzel,Index,SiebMax:Cardinal;
begin
Bis:=Round(sqrt(cFaktor*(StrToInt64(EditBis.Text))));
...
procedure TGesell.ButtonRechneClick(Sender: TObject);
var Start,Summe,von,bis,K,KMax:Int64;
Nr,L,KLaenge:Cardinal;
TSek:Extended;
Kette:Array of Int64;
KettenCnt : array of integer;
begin
if EditVon.Text='1' then EditVon.Text:='2';
if EditVon.Text='' then EditVon.Text:='2';
if EditBis.Text='' then EditBis.Text:='10000000';
von:=StrToInt64(EditVon.Text);
bis:=StrToInt64(EditBis.Text);
KLaenge:=StrToInt(EditKL.Text);
SetLength(Kette,KLaenge+1);
SetLength(KettenCnt,KLaenge+1);
LabelZeit.Caption:='';
Ausgabe.Clear;AusgabeF.Clear;AusgabeV.Clear;
Application.ProcessMessages;
Screen.Cursor:=crHourGlass;
Sieben;
TSek:=TimeSekunden;
K:=von;
while K<=bis do
begin
Start:=K;
KMax := cFaktor*K;
Nr:=0;
repeat
TeilerSumme(Start,Summe);
Start:=Summe;
inc(Nr);
Kette[Nr]:=Summe;
until (Start<=K) or (Start>KMax) or (Nr=KLaenge);
if Start=K then
begin
if Nr=1 then
AusgabeV.Lines.Add(IntToStr(K))
else
if Nr=2 then
AusgabeF.Lines.Add(IntToStr(Kette[2])+' - '+IntToStr(Kette[1]))
else
begin
Ausgabe.Lines.Add('Kettenlänge='+IntToStr(Nr)+' für '+IntToStr(K));
for L:=1 to Nr do
Ausgabe.Lines.Add(IntToStr(Kette[L]));
end;
inc(KettenCnt[Nr]);
end;
inc(K)
end;
TSek:=TimeSekunden-TSek;
Screen.Cursor:=crDefault;
LabelZeit.Caption:=Format('%0.2f',[TSek])+' sek.';
AusgabeV.Lines.Add('');
For Nr := 1 to KLaenge do
begin
L := KettenCnt[Nr];
IF L<> 0 then
begin
AusgabeV.Lines.Add(Format('Laenge %2d',[Nr]));
AusgabeV.Lines.Add(Format('Anzahl %2d',[L]));
end;
end;
end; |
Wenn Mathematiker eine obere Schranke für cFaktor= f(Größenordnung von K )angeben könnte.
Lelf ist doch Spezialist für Faktorisirung, der müsste doch sicher schnellere Verfahren kennen.
Vielleicht die 2 als Spezialfall einführen?Kommt sehr häufig vor und wäre nur Bit-Check.. And 1 = 0
und die Division nur ein Shr 1... Bringt sagenhafte 2% bei 2..1e7 ;-)
Gruß Horst
Horst_H - Mo 01.09.14 20:01
falsche Rubrik...
Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2025 by Christian Stelzmann Alle Rechte vorbehalten.
Alle Beiträge stammen von dritten Personen und dürfen geltendes Recht nicht verletzen.
Entwickler-Ecke und die zugehörigen Webseiten distanzieren sich ausdrücklich von Fremdinhalten jeglicher Art!