Autor |
Beitrag |
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
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
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.
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
Zuletzt bearbeitet von Horst_H am Di 12.08.14 19:42, insgesamt 1-mal bearbeitet
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Fiete
Beiträge: 601
Erhaltene Danke: 339
W7
Delphi 6 pro
|
Verfasst: 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
Einloggen, um Attachments anzusehen!
_________________ Fietes Gesetz: use your brain (THINK)
Für diesen Beitrag haben gedankt: Horst_H, Mathematiker
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: 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:
djm.cc/sociable.txt da kann man mit simplen Brute Force nicht gegen anstinken
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
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. 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
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
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: 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.
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:
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
Für diesen Beitrag haben gedankt: Fiete, Mathematiker
|
|
Fiete
Beiträge: 601
Erhaltene Danke: 339
W7
Delphi 6 pro
|
Verfasst: 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
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
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)
Zuletzt bearbeitet von Fiete am Fr 29.08.14 17:15, insgesamt 1-mal bearbeitet
Für diesen Beitrag haben gedankt: Horst_H, Mathematiker
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
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
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
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:
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
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
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: 601
Erhaltene Danke: 339
W7
Delphi 6 pro
|
Verfasst: 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
Gruß Fiete
_________________ Fietes Gesetz: use your brain (THINK)
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
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
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
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...
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
Beiträge: 601
Erhaltene Danke: 339
W7
Delphi 6 pro
|
Verfasst: Sa 30.08.14 13:21
Moin,
die Version vom Sonntag liegt jetzt in der finalen Version vor:
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
_________________ Fietes Gesetz: use your brain (THINK)
|
|
Horst_H
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: 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.
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
Beiträge: 1652
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Mo 01.09.14 20:01
|
|
|