Mandelbrotmengen mit ASM

Mathematiker · 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

markieren

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:

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:

CONST Drei : double = 3.0;
VAR   radius,a,b : double;
      n,i,it : word;
...
begin
      n:=5;     // z^5+z0
      radius:=9;
      it:=200;  //Iterationen
      //Koordinaten des Punktes stehen hier in a,b
...
      ASM
               FLD     Drei
               FLD     a
               FLD     b
               FLD     st        //hier wird wohl c aus a entnommen ?
               FLD     st(2)     //hier wohl d aus b ?
               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;
...  //in i steht die Anzahl der Iterationen bis zur Divergenz bzw. der Wert von it

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 · 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

Mathematiker · 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 · 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 )

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:

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

Blup · Beiträge: 173 Erhaltene Danke: 43

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:

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:

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;

Horst_H · 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

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:

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 //ohne packed wird das nichts
             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

  //a = c und b = d
  Fld d            //  d
  Fld c            //  c   d

  Fld St(1)        //  b   c   d
  Fld St(1)        //  a   b   c   d
  MOV  EAX,CntLimit
@Loop:
  dec  EAX
  FLD ST(1)        //  b   a   b   c   d
  FMUL ST(0),St(0) //b^2   a   b   c   d
  FLD  St(1)       //  a b^2   a   b   c   d
  FMUL ST(0),St(0) //a^2 b^2   a   b   c   d
//Test auf vorzeitigen Abbruch
  FLD  ST(1)       //b^2 a^2 b^2   a   b   c   d
  FADD ST(0),ST(1) //b^2+a^2 a^2 b^2   a   b   c   d
  fld   vier       //4  b^2+a^2 a^2 b^2   a   b   c   d
  fcomip St(0),St(1)     //b^2+a^2 a^2 b^2   a   b   c   d
  fstp  St(0)      //a^2 b^2   a   b   c   d
  jbe   @HabeFertig

//ab=  minuszehn * a2 * b2
  FLD  ST(0)       //a^2 a^2 b^2   a   b   c   d
  FMUL ST(0),ST(2)
  FMUL Minuszehn   // ab a^2 b^2   a   b   c   d

  FXCH ST(1)       //a^2  ab b^2   a   b   c   d

  FMUL ST(0),St(0) //a^4  ab b^2   a   b   c   d
  FXCH ST(2)       //b^2  ab a^4   a   b   c   d
  FMUL ST(0),St(0) //b^4  ab a^4   a   b   c   d

  //    b := ((fuenf * a2)+b2 + ab ) * b + d
  FXCH ST(2)       //a^4  ab b^4   a   b   c   d

  FLD   ST(0)      //a^4 a^4  ab b^4   a   b   c   d

  FMUL  fuenf       //5*a^4 a^4  ab b^4   a   b   c   d
  FADD  St(0),St(2) //+ ab
  FADD  St(0),St(3) //+b^4
  FMUL  St(0),St(5) //*  b
  FADD  St(0),St(7) //+  d
  FSTP  ST(5)       // Das neue b
                    //a^4  ab b^4   a   b   c   d
//a := ((fuenf * b2)+a2 + ab ) * a + c;
  FXCH  ST(2)      //b^4 ab a^4   a   b   c   d
  FMUL  fuenf      //5*b^4 ab a^4   a   b   c   d
  FADDP  St(1),St(0)//..+ab  a^4   a   b   c   d
  FADDP  St(1),St(0)//..+a^4   a   b   c   d
  FMULP  St(1),St(0)//*a   b   c   d
  FADD  St(0),St(2) //+c   b   c   d
                    //a    b   c   d
  test  EAX,EAX
  JNZ  @Loop

@HabeFertig:
  SUB  EAX,CntLimit
  Neg  EAX
  // Die FPU-Register  0..3 leeren
  Fcompp
  Fcompp
//FNSAVE Fsave;
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
// Vergleichslauf ASM und Neu berechnen gleichartig
// Aber in ASM wird intern nur extended gerechnet
  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.
{Abweichungen
-0.8020833 0.6093750  205  209
-0.7708333 0.3125000  360  359
-0.6927083-0.2135417  467  609
-0.6927083 0.2135417  430  513
-0.6875000-0.4218750  451  450
-0.6875000-0.1666667  729  762
-0.6822917-0.4218750  275  226
-0.4218750-0.6875000  451  450
-0.4218750-0.6822917  275  226
-0.2239583 0.7708333  224  238
-0.2135417-0.6927083  467  609
-0.2135417 0.6927083  762  689
-0.1666667-0.6875000  729  762
-0.1666667 0.6875000  770  768
 0.1666667 0.6875000  776  766
 0.2135417-0.6927083  430  513
 0.3125000-0.7708333  360  359
 0.4218750 0.6875000  383  384
 0.6093750-0.8020833  205  209
 0.6875000-0.1666667  770  768
 0.6875000 0.1666667  776  766
 0.6875000 0.4218750  383  384
 0.6927083-0.2135417  762  689
 0.7708333-0.2239583  224  238

IterOrg 1301.000 ms

IterAsm 803.000 ms

IterNeu 1251.000 ms}

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

Mathematiker · 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

markieren

Delphi-Quelltext

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

type
  tFSave = packed record //ohne packed wird das nichts
             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 · 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!

	Mitgliederliste
	Gruppen
	Das Team
	Richtlinien
	Synonyme