Autor |
Beitrag |
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: So 27.04.14 10:38
Hallo,
ich möchte Mandelbrotmengen höheren Grades zeichnen, konkret die 5.Grades, d.h.
(a+bî)^5 + c+îd = a^5 - 10·a³·b² + 5·a·b^4 + c + î·(5·a^4·b - 10·a²·b³+ b^5 + d)
Berechne ich 200 Iterationen je komplexe Zahl z0 = c+îd klassisch mit
Delphi-Quelltext 1: 2: 3: 4:
| quadrata:=a*a; quadratb:=b*b; a:=a*sqr(quadrata)-10.0*a*quadrata*quadratb+5.0*a*sqr(quadratb)+c; b:=5.0*sqr(quadrata)*b-10.0*quadrata*b*quadratb+b*sqr(quadratb)+d; |
so wird ein 660x660-Bild für -2 < c < 2 und -2 < d < 2 in 220 ms erstellt.
Vor Jahren habe ich eine ASM-Routine gefunden (Urheber weiß ich nicht mehr), die ich bisher nicht verwendet habe, da ich nichts verstehe.
Der wesentliche Teil ist:
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:
| CONST Drei : double = 3.0; VAR radius,a,b : double; n,i,it : word; ... begin n:=5; radius:=9; it:=200; ... ASM FLD Drei FLD a FLD b FLD st FLD st(2) XOR cx,cx @itloop: INC cx CMP cx,it JE @nopaint FLD st(1) FMUL st(2),st FLD st(1) FMUL st(2),st FPATAN FIMUL n FSINCOS FXCH st(3) FADDP st(2),st FXCH FSQRT FST radius MOV ax, n DEC ax @mult: FMUL radius DEC ax JNZ @mult FMUL st(2),st FMUL FADD st,st(2) FXCH FADD st,st(3) FLD radius FCOMP st(5) FSTSW ax AND ah,1 JNZ @itloop @nopaint: MOV i, cx FINIT END; ... |
Die Darstellung ist korrekt. Jedoch erhöht sich die Berechnungszeit auf knapp 1500 ms, d.h. gut das Sechsfache.
D.h., irgendwo in der ASM-Routine muss ein Fehler sein. Und ich habe keine Ahnung, da ich auf den ASM-Code schaue, wie die berühmte "Kuh ins Uhrwerk".
Assembler ohne Kommentare ist wahrscheinlich eine Zumutung.
Ich wäre Euch daher sehr dankbar, wenn doch jemand einmal einen Blick darauf werfen könnte. Irgendwo muss ein Fehler sein.
Danke und beste Grüße
Mathematiker
|
|
mandras
Beiträge: 430
Erhaltene Danke: 107
Win 10
Delphi 6 Prof, Delphi 10.4 Prof
|
Verfasst: So 27.04.14 11:46
Ich drösel den ASM-Code nun einmal nicht auf,
die höhere Zeiterfordernis ergibt sich aus der Anwendung von Arctan und sincos.
Hier hat wohl jemand gehofft, daß dies schneller sei als Potenzieren
(Multiplikation komplexer Zahlen alternativ als Multiplikation der Vektorlänge und "Weiterdrehen" des Winkels)
ArcTan und SinCos sind jedoch auch in der FPU quählend langsam da über Polynomtabellen realisiert.
Hinweise für den Zeitbedarf der einzelnen FPU-Befehle: www2.deec.uc.pt/~jlobo/tc/opcode_f.html
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: So 27.04.14 12:57
Hallo,
mandras hat folgendes geschrieben : | Hier hat wohl jemand gehofft, daß dies schneller sei als Potenzieren
(Multiplikation komplexer Zahlen alternativ als Multiplikation der Vektorlänge und "Weiterdrehen" des Winkels) |
Danke für den Hinweis. Genau das wird es sein.
Also werde ich es lieber der Delphi-Multiplikation überlassen.
Beste Grüße
Mathematiker
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: So 04.05.14 16:28
Hallo,
Delphi optimiert offensichtlich sehr gut.
EDIT: wegen falscher Berechnung
Ein Beispiel, wo unter freepascal 2.6.4 die Geschwindigkeit von Org zu Neu unter 32 Bit und 64 Bit recht genau wechselt
32 Bit 720/ 900 ms und unter
64-Bit ( xmm ) 907/717 ms
( i3-4330 3.5 Ghz )
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:
| program Mandel5; uses sysutils; const UgC =-2.0; UgD =-2.0; OgC = 2.0; OgD = 2.0;
deltac: double = (OgC-UgC)/660; deltad: double = (OgD-UgD)/660;
function IterOrg(c,d:double):integer; var a,a2,b,b2: double; begin result := 0; a := c; b := d; a2:=sqr(a); b2:=sqr(b); repeat inc(result); a:=a*sqr(a2)-10.0*a*a2*b2+5.0*a*sqr(b2)+c; b:=5.0*sqr(a2)*b-10.0*a2*b*b2+b*sqr(b2)+d; a2:=sqr(a); b2:=sqr(b); until (a2+b2>4.0) OR (result > 1000); end;
function IterNeu(c,d:double):integer; var a,a2,b,b2,b4: double; begin a := c; b := d; result := 0; a2:=a*a; b2:=b*b; repeat b4:= b2*b2; b2:= -10.0*b2+a2; a:=((b2 )*a2+5.0*b4)*a+c; b:=((b2+4*a2)*a2+ b4)*b+d; a2:=a*a; b2:=b*b; inc(result); until (a2+b2 >4.0) or (result>1000); end;
var c,d :double; T2,T1,T0 : TDateTime; cnt : Int64;
begin T0:= time; cnt := 0; c := UgC; repeat d := UgD; repeat inc(cnt,IterOrg(c,d)); d := d+deltaD; until d >= OgD; c := c+deltaC; until c >= OgC; writeln(cnt); T1:= time; cnt := 0; c := UgC; repeat d := UgD; repeat inc(cnt,IterNeu(c,d)); d := d+deltaD; until d >= OgD; c := c+deltaC; until c >= OgC; writeln(cnt); T2:= time; Writeln('Original ',(T1-T0)*86400.0*1000.0:6:3,' ms'); Writeln('Neu ',(T2-T1)*86400.0*1000.0:6:3,' ms'); writeln(cnt); end. |
Gruß Horst
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Blup
Beiträge: 174
Erhaltene Danke: 43
|
Verfasst: Di 06.05.14 16:11
Also ich vermute bei der Umformung der Gleichung zur Prozedur IterNeu hat sich ein Fehler eingeschlichen.
Vergleich doch mal die Ergenisse beider Prozeduren in einem Testdurchlauf.
Die Abbruchbedingung "(a2 + b2) > 4" kann man auch schon vor der ersten Iteration testen.
Bei mir sieht die Prozedur dann so aus:
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:
| function IterNeu(c, d: Double): Integer; var a, a2, b, b2, ab: Double; begin a := c; b := d; for n := 0 to 1000 do begin a2 := sqr(a); b2 := sqr(b); if (a2 + b2) > 4.0 then begin Result := n; Exit; end; ab := 10.0 * a2 * b2; a2 := sqr(a2); b2 := sqr(b2); a := (a2 - ab + (5.0 * b2)) * a + c; b := (b2 - ab + (5.0 * a2)) * b + d; end; Result := 1001; end; |
Für diesen Beitrag haben gedankt: Horst_H
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Di 06.05.14 22:35
Hallo,
eine Assemblerversion, die nicht wirklich viel schneller ist.
SSE müßte was bringen, da man zum Beispiel a*a und b*b in einem Schritt berechnen könnte.
EDIT: Neueste Version, wie unten in der Zip
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:
| program mandel5; {$IFDEF FPC} {$MODE Delphi} {$else} {$Apptype Console} {$ENDIF} uses sysutils; const CntLimit = 1001; vier : double = 4.0; fuenf : double = 5.0; minuszehn : double = -10.0;
UgC =-2.0; UgD =-2.0; OgC = 2.0; OgD = 2.0;
deltaX = 768; deltaY = 768;
deltac: double = (OgC-UgC)/deltaX; deltad: double = (OgD-UgD)/deltaY; type tfIter = function(c,d:double):integer;
tFSave = packed record Byte108 : array[0..27] of byte; ZahlFeld : array[0..7] of Extended; end;
function IterOrg(c,d:double):integer; var a,a2,b,b2: double; begin result := 0; a:= c; b:= d; a2:=sqr(a); b2:=sqr(b); repeat inc(result); a:=fuenf*a*sqr(b2)+minuszehn*a*a2*b2+a*sqr(a2)+c; b:=fuenf*b*sqr(a2)+minuszehn*b*a2*b2+b*sqr(b2)+d;
a2:=sqr(a); b2:=sqr(b); until (a2 + b2 >= vier) OR (result >= cntLimit); end;
function IterAsm(c,d:double):integer;assembler; asm FInit
Fld d Fld c Fld St(1) Fld St(1) MOV EAX,CntLimit @Loop: dec EAX FLD ST(1) FMUL ST(0),St(0) FLD St(1) FMUL ST(0),St(0) FLD ST(1) FADD ST(0),ST(1) fld vier fcomip St(0),St(1) fstp St(0) jbe @HabeFertig
FLD ST(0) FMUL ST(0),ST(2) FMUL Minuszehn FXCH ST(1) FMUL ST(0),St(0) FXCH ST(2) FMUL ST(0),St(0) FXCH ST(2) FLD ST(0) FMUL fuenf FADD St(0),St(2) FADD St(0),St(3) FMUL St(0),St(5) FADD St(0),St(7) FSTP ST(5) FXCH ST(2) FMUL fuenf FADDP St(1),St(0) FADDP St(1),St(0) FMULP St(1),St(0) FADD St(0),St(2) test EAX,EAX JNZ @Loop
@HabeFertig: SUB EAX,CntLimit Neg EAX Fcompp Fcompp end;
function IterNeu(c, d: Double): Integer; var a, a2, b, b2, ab: Double; begin a := c; b := d; for result := 1 to CntLimit do begin a2 := sqr(a); b2 := sqr(b); if (a2 + b2) >= vier then Exit; ab := minuszehn * a2 * b2; a2 := sqr(a2); b2 := sqr(b2); a := ((fuenf * b2)+a2 + ab ) * a + c; b := ((fuenf * a2)+b2 + ab ) * b + d; end; Result := CntLimit; end;
procedure TestLauf(Testfunk: tfIter); var c,d :double; T1,T0 : TDateTime; i,j : integer; cnt : LongWord; begin T0:= time; cnt := 0; For j :=0 to deltaY do begin d := UgD+j*deltaD; For i :=0 to deltaX do begin c := UgC+i*deltaC; cnt := TestFunk(c,d); end; end; T1:= time; Writeln((T1-T0)*86400.0*1000.0:6:3,' ms '); writeln; end;
var c,d :double; cnt1,cnt2 : LongWord;
begin Writeln('Abweichungen '); cnt1 := 0; cnt2 := 0; c := UgC; repeat d := UgD; repeat cnt1 := IterAsm(c,d); cnt2 := IterNeu(c,d); IF cnt1<>cnt2 then begin writeln(c:10:7,d:10:7,cnt1:5,cnt2:5); end; d := d+deltaD; until d >= OgD; c := c+deltaC; until c >= OgC; writeln; write('IterOrg ');TestLauf(IterOrg); write('IterAsm ');TestLauf(IterAsm); write('IterNeu ');TestLauf(IterNeu); {$IFNDEF Linux} readln {$ENDIF} end. |
Aber es gibt doch Verfahren, die nur die Grenzen abfahren, bzw in Gebiete > 1000 nicht weit eintauchen, wenn in der Zeile zuvor schon alles bei > 1000 war.
Gruß Horst
Zuletzt bearbeitet von Horst_H am Do 08.05.14 08:11, insgesamt 2-mal bearbeitet
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: Di 06.05.14 23:32
Hallo,
Danke für den neuen ASM-Text.
Leider kann ich ihn bei mir nicht richtig testen, da ich mit
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7:
| type tFSave = packed record Byte108 : array[0..27] of byte; ZahlFeld : array[0..7] of Extended; end; var FSave : tFSave; |
nicht zurecht komme. Wozu ist das da? Wo stehen nach der Rechnung die Werte a und b?
Wie schon gesagt, von ASM habe ich keine Ahnung. Eigentlich müsste ich es mir mal ansehen, aber ...
Beste Grüße
Mathematiker
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Mi 07.05.14 07:52
Hallo,
der Befehl Fsave speichert den gerade aktuellen Zustand der FPU ( in den ersten 28 (32/64 Bit System)) und deren Register St0..St7.
Ich habe das nur benutzt, um durch Auskommentierung die Zwischenschritte der Berechnung zu verfolgen, und auch, um mit Iterneu vergleichen zu können ( // writeln... ).
Es werden ja alle FPU-Register genutzt und es ist etwas eng..
Als Rückgabe dient ja eigentlich Register EAX.Welchen Wert a,b haben interessiert ja nicht mehr.
Es ging mir nur um die Größenordnung, die Assembler durch den Einsatz aller FPU Register bringen kann.
Das ist eben nicht der Bringer.
Wenn man www.mikusite.de/pages/x86.htm anschaut wird einem dort eine FPU Efficiency berechnet.
Efficiency is calculated like = (1000 Iterations/s) / MHz / Core
Ich habe genau 1 Sekunde für 1e5*1000 Berechnungen bei 1 Kern und und 3500 Mhz.
Also Eff-FPU = 1e5/3500/1 = 28,57.Auf der genannten Seite sind es in der letzten Version für i7-2600K Eff-FPU = 168, was schon erheblich schneller ist.Dort werden aber auch 3 Punkte ( real/ imaginär-Teil) in den FPU Register gehalten.
Eine Berechnung in der inneren Schleife braucht ja hier 29 Befehle, das ist ja nicht wenig.
Das sind hier 1.2 Takte pro FPU Befehl.Etwas lahm...ist ja noch nicht optimiert.Man kann sicher das Aufräumen optimieren ( FCOMpp ) und schon zum Ende hin FADDP einsetzen.
Gruß Horst
EDIT:
Die Version mit Test auf vorzeitigen Abbruch bei 1024x1024 Punkten, getestet mit fpc 2.6.4 unter Linux.
Mehr als 20..25% scheint es nicht zu bringen
Edit2:
Ein Testlauf mit Abbruch bei mehr als 200 Iterationen dauert bei mir IterAsm 181 ms, statt IterOrg 248 ms
Einloggen, um Attachments anzusehen!
Für diesen Beitrag haben gedankt: Mathematiker
|
|
|