* Using x86 and x87 assembly code in a GCC/gfortran environment.
@ 2018-03-23 21:32 Phil Seeger via fortran
2018-03-23 22:44 ` Jerry DeLisle
0 siblings, 1 reply; 3+ messages in thread
From: Phil Seeger via fortran @ 2018-03-23 21:32 UTC (permalink / raw)
To: fortran
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
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: Using x86 and x87 assembly code in a GCC/gfortran environment.
2018-03-23 21:32 Using x86 and x87 assembly code in a GCC/gfortran environment Phil Seeger via fortran
@ 2018-03-23 22:44 ` Jerry DeLisle
2018-03-23 23:21 ` Jerry DeLisle
0 siblings, 1 reply; 3+ messages in thread
From: Jerry DeLisle @ 2018-03-23 22:44 UTC (permalink / raw)
To: Phil Seeger, fortran
On 03/23/2018 02:32 PM, Phil Seeger via fortran wrote:
> 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.
>
I don't know if this is useful to you, but you might want to take a look
at gtk-fortran on github which will give you all the GUI stuff you would
ever need. There is even an example Mandelbrot program.
Also, you might try quadmath (kind=16) which gfortran supports.
I do not know the status of these tools in a Windows environment. I do
know that many applications run on windows using GTK
I have asked your questions on comp.lang.fortran?
Regards,
Jerry
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: Using x86 and x87 assembly code in a GCC/gfortran environment.
2018-03-23 22:44 ` Jerry DeLisle
@ 2018-03-23 23:21 ` Jerry DeLisle
0 siblings, 0 replies; 3+ messages in thread
From: Jerry DeLisle @ 2018-03-23 23:21 UTC (permalink / raw)
To: Phil Seeger, fortran
On 03/23/2018 03:44 PM, Jerry DeLisle wrote:
--- snip ---
>
> I have asked your questions on comp.lang.fortran?
Oops. Have you asked on comp.lang.fortran?
>
> Regards,
>
> Jerry
>
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2018-03-23 23:21 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2018-03-23 21:32 Using x86 and x87 assembly code in a GCC/gfortran environment Phil Seeger via fortran
2018-03-23 22:44 ` Jerry DeLisle
2018-03-23 23:21 ` Jerry DeLisle
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).