Using x86 and x87 assembly code in a GCC/gfortran environment.
Phil Seeger via fortran
fortran@gcc.gnu.org
Fri Mar 23 21:32:00 GMT 2018
Background. I have been using Fortran for more than 50 years. One of my
favorite projects was to compute the Mandelbrot set at magnifications up
to 10^10, maintaining precision by performing continuing iterations
totally within the x87 stack, with no memory references and no
round-offs, with a resolution of 19 significant digits. If the entire
Mandelbrot set were plotted at that magnification and at 150 dpi, the
width would be 1.6 million miles. A subroutine with the iteration loop
is appended.
The Problem. Because the program used GUI to get parameters for each
plot, and especially because I used the mouse to start and stop, the
program fails catastrophically under Win7 and later. The Microsoft API I
used to detect a mouse click no longer works! I have recently
decommissioned my last XP system, so I need to recompile. I have
installed MinGW (and gfortran), and have the program running with the
iteration loop written in REAL*8, but it looks a bit fuzzy.
Questions. Is there a (free) assembler? What does a fortran calling
sequence look like (I assume call-by-reference)? How to create a child
window for the plot? (Currently I write a BMP file and look at it
later.) How to make a GUI, with stuff like Dialog Boxes, Radio buttons,
etc.?
Thank you for any advice!
TITLE SUBROUTINE MANDEL(a, b, NMAX, NIT, ESCAPE)
;
;P. A. Seeger, Los Alamos National Laboratory, March 28, 1989
; (derived from MANDEL, March 28, 1988)
; 14 Nov 1991: keep x2 and y2 in stack (20% faster)
; 13 May 1992: fix error in ESCAPE for NIT > 32767
; 19 Apr 1996: inserted into template produced by MS PowerStation
Fortran 4.0
; 21 Mar 2018: always return ESCAPE time (delete AZIMUTH and OUTFILE)
;
;Iterate z := z*z + c, starting from z = 0 + 0*i, where c = a + b*i, for a
; maximum of NMAX iterations, quitting whenever abs(z*z) > 1000000.
Returns
; the number of iterations, NIT, and the ESCAPE time function defined by
; H.-O. Pietgen in Chapter 4.2 of "The Science of Fractal Images" [H.-O.
; Pietgen and D. Saupe, eds., Springer-Verlag, 1988]. Points which do not
; diverge are included in the Mandelbrot set, flagged by ESCAPE=0.
;
; a = Re(c), REAL*8, input
; b = Im(c), REAL*8, input
; NMAX = maximum number of iterations, INTEGER*2, input
; NIT = number of iterations done, INTEGER*2, output
; ESCAPE = "Escape Time" = -log2[(ln|z|)/2**NIT], REAL*4, output
;
; Iteration:
; Z := z*z + c
; X = Re(Z) := Re(z)*Re(z) - Im(z)*Im(z) + Re(c) = x*x - y*y + a
; Y = Im(Z) := 2*Re(z)*Im(z) + Im(c) = 2*x*y + b
;
.386P
.model FLAT
EXTRN __FFljj:NEAR
_DATA SEGMENT
_Mandel DB 'Mandel.asm', 00H
_DATA ENDS
PUBLIC _MANDEL@28
EXTRN __fltused:NEAR
_CONST SEGMENT
LARGE DD 1000000.0 ;value to compare with abs(z*z)
_CONST ENDS
;
_TEXT SEGMENT
_a$ = 8
_b$ = 12
_NMAX$ = 16
_NIT$ = 20
_ESCAPE$ = 24
_MANDEL@28 PROC NEAR
;
mov ecx, DWORD PTR _a$[esp-4] ;address of Re(c)
sub esp, 4
push ebx
fld QWORD PTR [ecx] ;load a=Re(c) in processor stack
push esi
mov esi, DWORD PTR _IMAG$[esp+8] ;address of Im(c)
push edi
push ebp
fld QWORD PTR [esi] ;load b=Im(c) in processor stack
;
FLD ST(1) ;initial Re(z) = x := Re(c) = a
FLD ST(1) ;initial Im(z) = y := Im(c) = b
;
FLD ST(1) ;x
FMUL ST(0),ST(0) ;x2 := x*x
mov ebp, DWORD PTR _NMAX$[esp+16] ;address of NMAX
FLD ST(1) ;y
FMUL ST(0),ST(0) ;y2 := y*y, stack is y2 x2 y x b a
;
mov ECX, DWORD PTR [ebp] ;NMAX to loop count register
XOR EDX,EDX ;set iteration count to zero
mov ebp, DWORD PTR _NIT$[esp+16] ;address of NIT
AGAIN: INC DX ;iteration counter
;Stack: y2 x2 y x b
a - -
FSUBP ST(1),ST(0) ; x2-y2 y x b a
- - -
FADD ST(0),ST(4) ; X y x b a
- - -
;
FXCH ST(2) ; x y X b a
- - -
FMULP ST(1),ST(0) ; x*y X b a -
- - -
FADD ST(0),ST(0) ; 2*x*y X b a -
- - -
FADD ST(0),ST(2) ; Y X b a -
- - -
;
FLD ST(1) ; X Y X b a
- - -
FMUL ST(0),ST(0) ; X2 Y X b a
- - -
FLD ST(1) ; Y X2 Y X b
a - -
FMUL ST(0),ST(0) ; Y2 X2 Y X b
a - -
;
FLD ST(1) ; X2 Y2 X2 Y X
b a -
FADD ST(0),ST(1) ; |Z*Z| Y2 X2 Y X
b a -
fcomp DWORD PTR LARGE ; Y2 X2 Y X b
a - -
; Stack is now ready for next loop
fnstsw ax ;get floating-point status word
test ah, 65 ; 00000041H
je SHORT DIVERGED ;abs(z*z) > 1000000.0
LOOPW AGAIN ; no -- try again
; completed NMAX loops without
diverging
mov DWORD PTR [ebp], edx ;return NIT to calling program
FLDZ ;set ESCAPE := 0 as flag for no
divergence
JMP SHORT DONE
;
DIVERGED: mov DWORD PTR [ebp], edx ;return NIT to calling program
;use intrinsic log functions to find "ESCAPE time"
FADDP ST(1),ST(0) ; |Z*Z| Y X b a
- - -
FLD1 ; 1 |Z*Z| Y X
b a - -
FLDLN2 ; ln(2) 1 |Z*Z| Y
X b a -
FLD1 ; 1 ln(2) 1 |Z*Z|
Y X b a
FXCH ST(3) ; |Z*Z| ln(2) 1 1
Y X b a
FYL2X ; 2*ln|Z| 1 1 Y X
b a -
;ln(2)*log2(z**2) = ln(z**2) = 2*ln|z|
FYL2X ; log2(2*ln|Z|) 1 Y X b
a - -
FSUBP ST(1),ST ; -log2(ln|Z|) Y X b a
- - -
;1.0-log2(2*ln|z|) = -log2(ln|z|)
FIADD DWORD PTR [ebp] ; ESCAPE on top of stack
;N-log2(ln|z|) = -log2[(ln|z|)/2**N]
;
DONE: mov DWORD PTR [ebp], edx ;return NIT to calling program
mov ecx, DWORD PTR _ESCAPE$[esp+16] ;address of ESCAPE
fstp DWORD PTR [ecx] ;return ESCAPE to calling program
pop ebp
pop edi
pop esi
pop ebx
finit
add esp, 4
ret 28
_MANDEL@28 ENDP
_TEXT ENDS
END
More information about the Fortran
mailing list