Solitaire
16th March 2009, 07:51 PM
With all the computer programming talk on this forum I found a site with dozens of different languages, called the Computer Language Shootout. Searching though the listings there I chose one that looked the friendliest and downloaded Free Basic. It’s a procedural programming language based upon QuickBasic for the old IBM machines.
The compiler sits in the Program Files area of the C drive and doesn’t do much except bring up a command prompt.
The IDE works like a word processor with line numbers on the left hand side. You first have to tell it where to find the compiler on the C drive, and then you can begin programming. I personally like the Nightshade Theme for the program. After writing the program you have to save it, compile it, and then run it by clicking on the icons on top.
The user manual on the other hand does not seem easy to use. Keywords don’t always appear in the index, you have to dig around. There’s not much written about them, just a few words followed by “WRITE ME” at the end. Write who I wonder, they didn’t leave an address.
According to the CLS website, it’s supposedly fun to take an algorithm and try to improve it. I took the following program and started tinkering, trying all sorts of different things from moving lines around to inserting different commands. Nothing really worked.
' The Computer Language Shootout
' http://shootout.alioth.debian.org/
' contributed by Antoni Gual
#include "crt.bi"
#define iter 50
dim shared As Integer w,h
const As Double limit=4.0
function calcmandel(byval x As Integer,byval y As Integer) As Integer
dim as double zr=0.0, zi=0.0, cr, ci, tr, ti
dim As Integer ii
cr = (2.0*x/w - 1.5) : ci=(2.0*y/h - 1.0)
for ii = 0 to iter-1
zi = 2.0*zr*zi + ci
zr = tr - ti + cr
tr=zr*zr : ti=zi*zi
if tr+ti > limit then return 0
next
return -1
end function
dim As Integer x, y, i,w1,bb,i1
dim b as uinteger
w = val(command$)
if w < 1 then w = 300
h = w
w1=w-1
printf("P4%c%d %d%c",10, w, h, 10)
'
for y = 0 to h-1
for x = 0 to w1 step 8
b=0:bb=&h80
i1=iif (x+7>w1,w1,x+7)
for i =x to i1
b or= bb and calcmandel(i,y)
bb shr=1
next
putchar b
next
next
end
Then I found the ASM command in the manual so I rewrote the algorithm in terms of ML, FPU, and XMM. Lots of false starts but I did get the program to actually beep faster. I wrote the following code adding just a few short comments.
Rem :: The Computer Language Shootout
Rem :: http://shootout.alioth.debian.org
Rem :: Contributed by John Lockard
Rem :: Date: 2009 March 15
Rem // Compiler Header
#Include "crt.bi"
Const As Double Escape = 4
Const As Integer Dwell = 50
Declare Sub Assemble_PBF(ByVal W As Integer, ByVal H As Integer)
Declare Sub Display_Test(ByVal W As Integer, ByVal H As Integer)
Declare Function Calculate_Mandelbrot(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Declare Function Calculate_Mandelbrot_FPU(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Declare Function Calculate_Mandelbrot_XMMS(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Declare Function Calculate_Mandelbrot_XMMD(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Rem // Main Program
Dim As Integer N
N = Val(Command$) ' [00]
If (N < 1) Then
N = 300
End If
Assemble_PBF(N, N)
'Display_Test(N, N)
Rem -- Program Note
Rem
Rem This program generates and returns a Portable Bitmap Format
Rem file containing a black and white an image of the Mandelbrot
Rem Set to the calling program for selected size N. For an N less
Rem than one, the default width and height becomes 300 pixels.
Rem
Rem -- Program Variable
Rem
Rem N - Size
Rem
Rem -- Program Footnote
Rem
Rem [00] User input for width and height of image file.
Rem // BMP File Assembler
Sub Assemble_PBF(ByVal W As Integer, ByVal H As Integer)
Dim As Integer B, C, I, J, K, L
PrintF("P4%c%d %d%c", 10, W, H, 10) ' [01]
For J = 0 To (H - 1) Step 1
For K = 0 To (W - 1) Step 8
L = IIf((K + 7) > (W - 1), (W - 1), (K + 7)) ' [02]
C = 0
B = 128
For I = K To L Step 1
C = C Or (B And Calculate_Mandelbrot_FPU(I, J, W, H))' [03]
B = B Shr 1 ' [04]
Next
PutChar C ' [05]
Next
Next
Rem -- Procedural Notes
Rem
Rem This procedure produces a Portable Bitmap Format file
Rem containing the selected portion of the Manderbolt Set.
Rem Each bit set in the image corresponds to a point in
Rem or near the Mandelbrot Set.
Rem
Rem Since PBF files store the bottom of the image first and
Rem this subroutine starts with the negative imaginary axis
Rem first, the output file will display the Mandelbrot Set
Rem with positive imaginary axis on top.
Rem
Rem -- Procedural Variables
Rem
Rem B - Bit Position
Rem C - Character
Rem
Rem W - Image Width
Rem H - Image Height
Rem
Rem I - Horizontal Array Position
Rem J - Vertical Array Position
Rem
Rem K - First Bit Position In Horizontal Array
Rem L - Last Bit Position In Horizontal Array
Rem
Rem -- Procedural Footnotes
Rem
Rem [01] Fast C code output routine from the "crt.bi" file.
Rem Write the header for the PBF file to the calling program.
Rem [02] A very compact form of:
Rem If ((K + 7) > (W - 1)) Then
Rem L = (W - 1)
Rem Else
Rem L = (K + 7)
Rem End If
Rem [03] Integer bit math, sets bit in character if a member.
Rem [04] A faster version of:
Rem B = B / 2
Rem [05] Fast C code output routine from the "crt.bi" file.
Rem Writes characters to the calling program.
End Sub
Rem // Display Test Set
Sub Display_Test(ByVal W As Integer, ByVal H As Integer)
Dim As Double Ci, Cr
Dim As Integer I, J
For J = (H - 1) To 0 Step -1
For I = 0 To (W - 1) Step 1
If (Calculate_Mandelbrot_FPU(I, J, W, H) = -1) Then
Print "*";
Else
Print " ";
End If
Next
Print
Next
Sleep ' [06]
Rem -- Procedural Note
Rem
Rem This procedure displays the Mandelbrot Set as asterisks
Rem within the Console Window. Not recommended for widths
Rem greater than 79 characters.
Rem
Rem -- Procedural Variables
Rem
Rem I - Horizontal Position
Rem J - Vertical Position
Rem
Rem W - Width
Rem H - Height
Rem
Rem -- Procedural Footnotes
Rem
Rem [06] This command tells the computer to wait for user input.
End Sub
Rem // Mandelbrot Set Membership Calculation
Function Calculate_Mandelbrot(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Double Ci, Cr, Si, Sr, Zi, Zr
Dim As Integer C
Ci = (J / H) * 2.0 - 1.0 ' [07]
Cr = (I / W) * 2.0 - 1.5 ' [08]
For C = 0 To (Dwell - 1) Step 1
Zi = Ci + Zr * Zi * 2.0
Zr = Cr + Sr - Si
Si = Zi * Zi
Sr = Zr * Zr
If ((Sr + Si) > Escape) Then Return 0
Next
Return -1
Rem -- Procedural Note
Rem
Rem This procedure determines if a point exists in or near
Rem the Mandelbrot Set by iterating a complex number for a
Rem predetermined number of steps. The larger the dwell limit
Rem the more accurately the better the returned value, -1,
Rem represents a point within the set.
Rem
Rem -- Procedural Variables
Rem
Rem C - Iteration
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [07] This equation translates vertical positions in the
Rem image grid into imaginary values on the complex plain.
Rem
Rem The equation starts by dividing J in pixels by H in pixels
Rem producing a dimensionless scalar. The equation then multiplies
Rem this dimensionless scalar by an imaginary range value of 2.0
Rem producing an imaginary offset. Finally, the equation adds
Rem the imaginary offset to the imaginary anchor point of -1.0
Rem producing the imaginary position part of the complex number.
Rem
Rem [08] This equation translates horizontal positions in the
Rem image grid into real values on the complex plain.
Rem
Rem The equation starts by dividing I in pixels by W in pixels
Rem producing a dimensionless scalar. The equation then multiplies
Rem this dimensionless scalar by a real range value of 2.0
Rem producing an real offset. Finally, the equation adds
Rem the real offset to the real anchor point of -1.5
Rem producing the real position part of the complex number.
End Function
Rem // Mandelbrot Set Membership Calculation - FPU Version
Function Calculate_Mandelbrot_FPU(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Double Ci = (J / H) * 2.0 - 1.0
Dim As Double Cr = (I / W) * 2.0 - 1.5
Dim As Double E = Escape
Dim As Integer D = Dwell
ASM
REM --------------------------- ' Load CPU Extended Registers
MOV EAX, 0 ' A = False
MOV EBX, 1 ' B = 1
MOV ECX, 0 ' C = 0
MOV EDX, [D] ' D = Dwell
REM --------------------------- ' Load FPU Stack Registers
FLD QWORD PTR [Cr] ' ST[7] = Cr
FLD QWORD PTR [Ci] ' ST[6] = Ci
FLDZ ' ST[5] = Sr
FLDZ ' ST[4] = Si
FLD ST(1) ' ST[3] = Zr
FLD ST(1) ' ST[2] = Zi
FLD QWORD PTR [E] ' ST[1] = Escape
FLDZ ' ST[0] = Mu
REM --------------------------- ' Sr = (Zr * Zr) ' [09]
L1: FXCH ST(3) ' ST[0] = Zr
FST ST(3) ' ST[3] = Zr
FMUL ST(0), ST(0) ' ST[0] = Zr * Zr
FST ST(5) ' ST[5] = Zr * Zr
REM --------------------------- ' Si = (Zi * Zi)
FXCH ST(2) ' ST[0] = Zi
FST ST(2) ' ST[2] = Zi
FMUL ST(0), ST(0) ' ST[0] = Zi * Zi
FST ST(4) ' ST[4] = Zi * Zi
REM --------------------------- ' Mu = (Sr + Si)
FADD ST(5) ' ST[0] = (Si + Sr)
REM --------------------------- ' If ((Mu > Escape) Then L2
FCOMI ST(0), ST(1) '
JA L2 '
REM --------------------------- ' Zi = 2 * (Zi * Zr) + Ci
FXCH ST(2) ' ST[0] = Zi
FMUL ST(3) ' ST[0] = Zi * Zr
FADD ST(0) ' ST[0] = 2 * Zi * Zr
FADD ST(6) ' ST[0] = 2 * Zi * Zr + Ci
FXCH ST(2) ' ST[2] = 2 * Zi * Zr + Ci
REM --------------------------- ' Zr = (Sr - Si) + Cr
FXCH ST(5) ' ST[0] = Sr
FSUB ST(4) ' ST[0] = Sr - Si
FADD ST(7) ' ST[0] = Sr - Si + Cr
FXCH ST(3) ' ST[2] = Sr - Si + Cr
REM --------------------------- ' C = C + 1
ADD ECX, EBX '
REM --------------------------- ' If (C <> Dwell) Then L1
CMP ECX, EDX '
JNE L1 '
REM --------------------------- ' A = True
SUB EAX, EBX '
REM --------------------------- ' Empty FPU Stack Registers ' [10]
L2: FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
REM --------------------------- ' Return A
MOV [Function], EAX '
END ASM
Rem -- Procedural Note
Rem
Rem This procedure, using the Floating Point Unit, determines
Rem if a point exists in or near the Mandelbrot Set by iterating
Rem a complex number for a predetermined number of steps.
Rem
Rem -- Procedural Variables
Rem
Rem A - Accumulator
Rem B - One
Rem C - Iteration
Rem D - Dwell
Rem E - Escape
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [09] The FLD ST(I) instruction pushes new items onto the stack
Rem instead of copying, as one would expect scrambling the
Rem stack. Using FXCH ST(I) and FST ST(I) prevents the stack
Rem from being overloaded, but takes extra time for the copying
Rem back of stored values.
Rem
Rem [10] Floating-point register must be emptied after use.
Rem Otherwise, the stack will become misaligned producing
Rem numerous errors.
End Function
Rem // Mandelbrot Set Membership Calculation - XMM Single Version
Function Calculate_Mandelbrot_XMMS(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Single Zero = 0
Dim As Single Ci = (J / H) * 2.0 - 1.0
Dim As Single Cr = (I / W) * 2.0 - 1.5
Dim As Single E = Escape
Dim As Integer D = Dwell
ASM
REM --------------------------- ' Load FPU Stack Registers ' [11]
FLDZ ' ST[3] A = False
FLD1 ' ST[2] B = 1
FILD DWORD PTR [D] ' ST[1] D = Dwell
FLDZ ' ST[0] C = 0
REM --------------------------- ' Load XMM Registers ' [12]
MOVSS XMM7, DWORD [Cr] ' XMM[7] = Cr
MOVSS XMM6, DWORD [Ci] ' XMM[6] = Ci
MOVSS XMM5, DWORD [Zero] ' XMM[5] = Sr
MOVSS XMM4, DWORD [Zero] ' XMM[4] = Si
MOVSS XMM3, DWORD [Cr] ' XMM[3] = Zr
MOVSS XMM2, DWORD [Ci] ' XMM[2] = Zi
MOVSS XMM1, DWORD [E] ' XMM[1] = Escape
MOVSS XMM0, DWORD [Zero] ' XMM[0] = Mu
REM --------------------------- ' Sr = (Zr * Zr)
L3: MOVSS XMM5, XMM3 ' XMM[5] = Zr
MULSS XMM5, XMM5 ' XMM[5] = Zr * Zr
REM --------------------------- ' Si = (Zi * Zi)
MOVSS XMM4, XMM2 ' XMM[4] = Zi
MULSS XMM4, XMM4 ' XMM[4] = Zi * Zi
REM --------------------------- ' M = (Sr + Si)
MOVSS XMM0, XMM5 ' XMM[0] = Sr
ADDSS XMM0, XMM4 ' XMM[0] = (Si + Sr)
REM --------------------------- ' If ((Mu > Escape) Then L4
COMISS XMM0, XMM1 '
JA L4 '
REM --------------------------- ' Zi = 2 * (Zi * Zr) + Ci
MULSS XMM2, XMM3 ' XMM[2] = Zi * Zr
ADDSS XMM2, XMM2 ' XMM[2] = 2 * Zi * Zr
ADDSS XMM2, XMM6 ' XMM[2] = 2 * Zi * Zr + Ci
REM --------------------------- ' Zr = Cr + (Zr - Zi)
MOVSS XMM3, XMM5 ' XMM[3] = Sr
SUBSS XMM3, XMM4 ' XMM[3] = Sr - Si
ADDSS XMM3, XMM7 ' XMM[3] = Sr - Si + Cr
REM --------------------------- ' C = C + 1
FADD ST(2) '
REM --------------------------- ' If (C <> Dwell) Then L3
FCOMI ST, ST(1) '
JNE L3 '
REM --------------------------- ' A = True
FXCH ST(3) ' ST[0] = A
FSUB ST(2) ' ST[0] = False
FXCH ST(3) ' ST[3] = A
REM --------------------------- ' Empty FPU Stack Registers
L4: FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
Rem --------------------------- ' Return A
FISTP DWORD PTR [Function] '
END ASM
Rem // Procedural Note
Rem
Rem This procedure, using the first Supplemental Streaming Extension
Rem set, determines if a point exists in or near the Mandelbrot Set
Rem by iterating a complex number for a predetermined number of steps.
Rem
Rem -- Procedural Variables
Rem
Rem A - Accumulator
Rem B - One
Rem C - Iteration
Rem D - Dwell
Rem E - Escape
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [11] In SIMD mode the general purpose registers EAX, EBX, ECX,
Rem and EDX store instructions for SSE1 operations, requiring
Rem the use of the floating point unit for both iteration and
Rem membership state.
Rem
Rem [12] Use of single floating points in this procedure limits
Rem the resolution of the set to one in eight million parts over
Rem the range of -2.0 to +2.0 in both imaginary and real dimensions.
End Function
Rem // Mandelbrot Set Membership Calculation - XMM Double Version
Function Calculate_Mandelbrot_XMMD(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Double Zero = 0
Dim As Double Ci = (J / H) * 2.0 - 1.0
Dim As Double Cr = (I / W) * 2.0 - 1.5
Dim As Double E = Escape
Dim As Integer D = Dwell
ASM
REM --------------------------- ' Load FPU Stack Registers ' [13]
FLDZ ' ST[3] A = 0
FLD1 ' ST[2] B = 1
FILD DWORD PTR [D] ' ST[1] D = Dwell
FLDZ ' ST[0] C = 0
REM --------------------------- ' Load FPU Stack Registers ' [14]
MOVSD XMM7, QWORD [Cr] ' XMM[7] = Cr
MOVSD XMM6, QWORD [Ci] ' XMM[6] = Ci
MOVSD XMM5, QWORD [Zero] ' XMM[5] = Sr
MOVSD XMM4, QWORD [Zero] ' XMM[4] = Si
MOVSD XMM3, QWORD [Cr] ' XMM[3] = Zr
MOVSD XMM2, QWORD [Ci] ' XMM[2] = Zi
MOVSD XMM1, QWORD [E] ' XMM[1] = Escape
MOVSD XMM0, QWORD [Zero] ' XMM[0] = Mu
REM --------------------------- ' Sr = (Zr * Zr)
L5: MOVSD XMM5, XMM3 ' XMM[5] = Zr
MULSD XMM5, XMM5 ' XMM[5] = Zr * Zr
REM --------------------------- ' Si = (Zi * Zi)
MOVSD XMM4, XMM2 ' XMM[4] = Zi
MULSD XMM4, XMM4 ' XMM[4] = Zi * Zi
REM --------------------------- ' M = (Sr + Si)
MOVSD XMM0, XMM5 ' XMM[0] = Sr
ADDSD XMM0, XMM4 ' XMM[0] = (Si + Sr)
REM --------------------------- ' If ((Mu > Escape) Then L5
COMISD XMM0, XMM1 '
JA L6 '
REM --------------------------- ' Zi = 2 * (Zi * Zr) + Ci
MULSD XMM2, XMM3 ' XMM[2] = Zi * Zr
ADDSD XMM2, XMM2 ' XMM[2] = 2 * Zi * Zr
ADDSD XMM2, XMM6 ' XMM[2] = 2 * Zi * Zr + Ci
REM --------------------------- ' Zr = Cr + (Zr - Zi)
MOVSD XMM3, XMM5 ' XMM[3] = Sr
SUBSD XMM3, XMM4 ' XMM[3] = Sr - Si
ADDSD XMM3, XMM7 ' XMM[3] = Sr - Si + Cr
REM --------------------------- ' C = C + 1
FADD ST(2) '
REM --------------------------- ' If (C <> Dwell) Then L5
FCOMI ST, ST(1) '
JNE L5 '
REM --------------------------- ' A = True
FXCH ST(3) ' ST[0] = A
FSUB ST(2) ' ST[0] = False
FXCH ST(3) ' ST[3] = A
REM ----------------------------' Empty FPU Stack Registers
L6: FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
Rem --------------------------- ' Return A
FISTP DWORD PTR [Function] '
END ASM
Rem // Procedural Note
Rem
Rem This procedure, using the second Supplemental Streaming Extension
Rem set, determines if a point exists in or near the Mandelbrot Set
Rem by iterating a complex number for a predetermined number of steps.
Rem
Rem -- Procedural Variables
Rem
Rem A - Accumulator
Rem B - One
Rem C - Iteration
Rem D - Dwell
Rem E - Escape
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [13] In SIMD mode the general purpose registers EAX, EBX, ECX,
Rem and EDX store instructions for SSE2 operations requiring
Rem the use of the floating point unit for both iteration and
Rem membership state.
Rem
Rem [14] Use of double floating points in this procedure limits
Rem the resolution of the set to one in four quadrillion parts
Rem over the range of -2.0 to +2.0 in both imaginary and real
Rem dimensions.
End Function
I’m not sure an actual computer programmer would take this approach to the problem, but it was interesting to work with the most basic level of computer programming short of writing everything in ones and zeros.
P.S. Here’s a list of the websites.
The Computer Language Shootout (http://shootout.alioth.debian.org/)
Free Basic (http://www.freebasic.net/)
Intel Software Development Manuals (http://www.intel.com/products/processor/manuals/)
The compiler sits in the Program Files area of the C drive and doesn’t do much except bring up a command prompt.
The IDE works like a word processor with line numbers on the left hand side. You first have to tell it where to find the compiler on the C drive, and then you can begin programming. I personally like the Nightshade Theme for the program. After writing the program you have to save it, compile it, and then run it by clicking on the icons on top.
The user manual on the other hand does not seem easy to use. Keywords don’t always appear in the index, you have to dig around. There’s not much written about them, just a few words followed by “WRITE ME” at the end. Write who I wonder, they didn’t leave an address.
According to the CLS website, it’s supposedly fun to take an algorithm and try to improve it. I took the following program and started tinkering, trying all sorts of different things from moving lines around to inserting different commands. Nothing really worked.
' The Computer Language Shootout
' http://shootout.alioth.debian.org/
' contributed by Antoni Gual
#include "crt.bi"
#define iter 50
dim shared As Integer w,h
const As Double limit=4.0
function calcmandel(byval x As Integer,byval y As Integer) As Integer
dim as double zr=0.0, zi=0.0, cr, ci, tr, ti
dim As Integer ii
cr = (2.0*x/w - 1.5) : ci=(2.0*y/h - 1.0)
for ii = 0 to iter-1
zi = 2.0*zr*zi + ci
zr = tr - ti + cr
tr=zr*zr : ti=zi*zi
if tr+ti > limit then return 0
next
return -1
end function
dim As Integer x, y, i,w1,bb,i1
dim b as uinteger
w = val(command$)
if w < 1 then w = 300
h = w
w1=w-1
printf("P4%c%d %d%c",10, w, h, 10)
'
for y = 0 to h-1
for x = 0 to w1 step 8
b=0:bb=&h80
i1=iif (x+7>w1,w1,x+7)
for i =x to i1
b or= bb and calcmandel(i,y)
bb shr=1
next
putchar b
next
next
end
Then I found the ASM command in the manual so I rewrote the algorithm in terms of ML, FPU, and XMM. Lots of false starts but I did get the program to actually beep faster. I wrote the following code adding just a few short comments.
Rem :: The Computer Language Shootout
Rem :: http://shootout.alioth.debian.org
Rem :: Contributed by John Lockard
Rem :: Date: 2009 March 15
Rem // Compiler Header
#Include "crt.bi"
Const As Double Escape = 4
Const As Integer Dwell = 50
Declare Sub Assemble_PBF(ByVal W As Integer, ByVal H As Integer)
Declare Sub Display_Test(ByVal W As Integer, ByVal H As Integer)
Declare Function Calculate_Mandelbrot(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Declare Function Calculate_Mandelbrot_FPU(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Declare Function Calculate_Mandelbrot_XMMS(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Declare Function Calculate_Mandelbrot_XMMD(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Rem // Main Program
Dim As Integer N
N = Val(Command$) ' [00]
If (N < 1) Then
N = 300
End If
Assemble_PBF(N, N)
'Display_Test(N, N)
Rem -- Program Note
Rem
Rem This program generates and returns a Portable Bitmap Format
Rem file containing a black and white an image of the Mandelbrot
Rem Set to the calling program for selected size N. For an N less
Rem than one, the default width and height becomes 300 pixels.
Rem
Rem -- Program Variable
Rem
Rem N - Size
Rem
Rem -- Program Footnote
Rem
Rem [00] User input for width and height of image file.
Rem // BMP File Assembler
Sub Assemble_PBF(ByVal W As Integer, ByVal H As Integer)
Dim As Integer B, C, I, J, K, L
PrintF("P4%c%d %d%c", 10, W, H, 10) ' [01]
For J = 0 To (H - 1) Step 1
For K = 0 To (W - 1) Step 8
L = IIf((K + 7) > (W - 1), (W - 1), (K + 7)) ' [02]
C = 0
B = 128
For I = K To L Step 1
C = C Or (B And Calculate_Mandelbrot_FPU(I, J, W, H))' [03]
B = B Shr 1 ' [04]
Next
PutChar C ' [05]
Next
Next
Rem -- Procedural Notes
Rem
Rem This procedure produces a Portable Bitmap Format file
Rem containing the selected portion of the Manderbolt Set.
Rem Each bit set in the image corresponds to a point in
Rem or near the Mandelbrot Set.
Rem
Rem Since PBF files store the bottom of the image first and
Rem this subroutine starts with the negative imaginary axis
Rem first, the output file will display the Mandelbrot Set
Rem with positive imaginary axis on top.
Rem
Rem -- Procedural Variables
Rem
Rem B - Bit Position
Rem C - Character
Rem
Rem W - Image Width
Rem H - Image Height
Rem
Rem I - Horizontal Array Position
Rem J - Vertical Array Position
Rem
Rem K - First Bit Position In Horizontal Array
Rem L - Last Bit Position In Horizontal Array
Rem
Rem -- Procedural Footnotes
Rem
Rem [01] Fast C code output routine from the "crt.bi" file.
Rem Write the header for the PBF file to the calling program.
Rem [02] A very compact form of:
Rem If ((K + 7) > (W - 1)) Then
Rem L = (W - 1)
Rem Else
Rem L = (K + 7)
Rem End If
Rem [03] Integer bit math, sets bit in character if a member.
Rem [04] A faster version of:
Rem B = B / 2
Rem [05] Fast C code output routine from the "crt.bi" file.
Rem Writes characters to the calling program.
End Sub
Rem // Display Test Set
Sub Display_Test(ByVal W As Integer, ByVal H As Integer)
Dim As Double Ci, Cr
Dim As Integer I, J
For J = (H - 1) To 0 Step -1
For I = 0 To (W - 1) Step 1
If (Calculate_Mandelbrot_FPU(I, J, W, H) = -1) Then
Print "*";
Else
Print " ";
End If
Next
Next
Sleep ' [06]
Rem -- Procedural Note
Rem
Rem This procedure displays the Mandelbrot Set as asterisks
Rem within the Console Window. Not recommended for widths
Rem greater than 79 characters.
Rem
Rem -- Procedural Variables
Rem
Rem I - Horizontal Position
Rem J - Vertical Position
Rem
Rem W - Width
Rem H - Height
Rem
Rem -- Procedural Footnotes
Rem
Rem [06] This command tells the computer to wait for user input.
End Sub
Rem // Mandelbrot Set Membership Calculation
Function Calculate_Mandelbrot(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Double Ci, Cr, Si, Sr, Zi, Zr
Dim As Integer C
Ci = (J / H) * 2.0 - 1.0 ' [07]
Cr = (I / W) * 2.0 - 1.5 ' [08]
For C = 0 To (Dwell - 1) Step 1
Zi = Ci + Zr * Zi * 2.0
Zr = Cr + Sr - Si
Si = Zi * Zi
Sr = Zr * Zr
If ((Sr + Si) > Escape) Then Return 0
Next
Return -1
Rem -- Procedural Note
Rem
Rem This procedure determines if a point exists in or near
Rem the Mandelbrot Set by iterating a complex number for a
Rem predetermined number of steps. The larger the dwell limit
Rem the more accurately the better the returned value, -1,
Rem represents a point within the set.
Rem
Rem -- Procedural Variables
Rem
Rem C - Iteration
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [07] This equation translates vertical positions in the
Rem image grid into imaginary values on the complex plain.
Rem
Rem The equation starts by dividing J in pixels by H in pixels
Rem producing a dimensionless scalar. The equation then multiplies
Rem this dimensionless scalar by an imaginary range value of 2.0
Rem producing an imaginary offset. Finally, the equation adds
Rem the imaginary offset to the imaginary anchor point of -1.0
Rem producing the imaginary position part of the complex number.
Rem
Rem [08] This equation translates horizontal positions in the
Rem image grid into real values on the complex plain.
Rem
Rem The equation starts by dividing I in pixels by W in pixels
Rem producing a dimensionless scalar. The equation then multiplies
Rem this dimensionless scalar by a real range value of 2.0
Rem producing an real offset. Finally, the equation adds
Rem the real offset to the real anchor point of -1.5
Rem producing the real position part of the complex number.
End Function
Rem // Mandelbrot Set Membership Calculation - FPU Version
Function Calculate_Mandelbrot_FPU(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Double Ci = (J / H) * 2.0 - 1.0
Dim As Double Cr = (I / W) * 2.0 - 1.5
Dim As Double E = Escape
Dim As Integer D = Dwell
ASM
REM --------------------------- ' Load CPU Extended Registers
MOV EAX, 0 ' A = False
MOV EBX, 1 ' B = 1
MOV ECX, 0 ' C = 0
MOV EDX, [D] ' D = Dwell
REM --------------------------- ' Load FPU Stack Registers
FLD QWORD PTR [Cr] ' ST[7] = Cr
FLD QWORD PTR [Ci] ' ST[6] = Ci
FLDZ ' ST[5] = Sr
FLDZ ' ST[4] = Si
FLD ST(1) ' ST[3] = Zr
FLD ST(1) ' ST[2] = Zi
FLD QWORD PTR [E] ' ST[1] = Escape
FLDZ ' ST[0] = Mu
REM --------------------------- ' Sr = (Zr * Zr) ' [09]
L1: FXCH ST(3) ' ST[0] = Zr
FST ST(3) ' ST[3] = Zr
FMUL ST(0), ST(0) ' ST[0] = Zr * Zr
FST ST(5) ' ST[5] = Zr * Zr
REM --------------------------- ' Si = (Zi * Zi)
FXCH ST(2) ' ST[0] = Zi
FST ST(2) ' ST[2] = Zi
FMUL ST(0), ST(0) ' ST[0] = Zi * Zi
FST ST(4) ' ST[4] = Zi * Zi
REM --------------------------- ' Mu = (Sr + Si)
FADD ST(5) ' ST[0] = (Si + Sr)
REM --------------------------- ' If ((Mu > Escape) Then L2
FCOMI ST(0), ST(1) '
JA L2 '
REM --------------------------- ' Zi = 2 * (Zi * Zr) + Ci
FXCH ST(2) ' ST[0] = Zi
FMUL ST(3) ' ST[0] = Zi * Zr
FADD ST(0) ' ST[0] = 2 * Zi * Zr
FADD ST(6) ' ST[0] = 2 * Zi * Zr + Ci
FXCH ST(2) ' ST[2] = 2 * Zi * Zr + Ci
REM --------------------------- ' Zr = (Sr - Si) + Cr
FXCH ST(5) ' ST[0] = Sr
FSUB ST(4) ' ST[0] = Sr - Si
FADD ST(7) ' ST[0] = Sr - Si + Cr
FXCH ST(3) ' ST[2] = Sr - Si + Cr
REM --------------------------- ' C = C + 1
ADD ECX, EBX '
REM --------------------------- ' If (C <> Dwell) Then L1
CMP ECX, EDX '
JNE L1 '
REM --------------------------- ' A = True
SUB EAX, EBX '
REM --------------------------- ' Empty FPU Stack Registers ' [10]
L2: FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
REM --------------------------- ' Return A
MOV [Function], EAX '
END ASM
Rem -- Procedural Note
Rem
Rem This procedure, using the Floating Point Unit, determines
Rem if a point exists in or near the Mandelbrot Set by iterating
Rem a complex number for a predetermined number of steps.
Rem
Rem -- Procedural Variables
Rem
Rem A - Accumulator
Rem B - One
Rem C - Iteration
Rem D - Dwell
Rem E - Escape
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [09] The FLD ST(I) instruction pushes new items onto the stack
Rem instead of copying, as one would expect scrambling the
Rem stack. Using FXCH ST(I) and FST ST(I) prevents the stack
Rem from being overloaded, but takes extra time for the copying
Rem back of stored values.
Rem
Rem [10] Floating-point register must be emptied after use.
Rem Otherwise, the stack will become misaligned producing
Rem numerous errors.
End Function
Rem // Mandelbrot Set Membership Calculation - XMM Single Version
Function Calculate_Mandelbrot_XMMS(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Single Zero = 0
Dim As Single Ci = (J / H) * 2.0 - 1.0
Dim As Single Cr = (I / W) * 2.0 - 1.5
Dim As Single E = Escape
Dim As Integer D = Dwell
ASM
REM --------------------------- ' Load FPU Stack Registers ' [11]
FLDZ ' ST[3] A = False
FLD1 ' ST[2] B = 1
FILD DWORD PTR [D] ' ST[1] D = Dwell
FLDZ ' ST[0] C = 0
REM --------------------------- ' Load XMM Registers ' [12]
MOVSS XMM7, DWORD [Cr] ' XMM[7] = Cr
MOVSS XMM6, DWORD [Ci] ' XMM[6] = Ci
MOVSS XMM5, DWORD [Zero] ' XMM[5] = Sr
MOVSS XMM4, DWORD [Zero] ' XMM[4] = Si
MOVSS XMM3, DWORD [Cr] ' XMM[3] = Zr
MOVSS XMM2, DWORD [Ci] ' XMM[2] = Zi
MOVSS XMM1, DWORD [E] ' XMM[1] = Escape
MOVSS XMM0, DWORD [Zero] ' XMM[0] = Mu
REM --------------------------- ' Sr = (Zr * Zr)
L3: MOVSS XMM5, XMM3 ' XMM[5] = Zr
MULSS XMM5, XMM5 ' XMM[5] = Zr * Zr
REM --------------------------- ' Si = (Zi * Zi)
MOVSS XMM4, XMM2 ' XMM[4] = Zi
MULSS XMM4, XMM4 ' XMM[4] = Zi * Zi
REM --------------------------- ' M = (Sr + Si)
MOVSS XMM0, XMM5 ' XMM[0] = Sr
ADDSS XMM0, XMM4 ' XMM[0] = (Si + Sr)
REM --------------------------- ' If ((Mu > Escape) Then L4
COMISS XMM0, XMM1 '
JA L4 '
REM --------------------------- ' Zi = 2 * (Zi * Zr) + Ci
MULSS XMM2, XMM3 ' XMM[2] = Zi * Zr
ADDSS XMM2, XMM2 ' XMM[2] = 2 * Zi * Zr
ADDSS XMM2, XMM6 ' XMM[2] = 2 * Zi * Zr + Ci
REM --------------------------- ' Zr = Cr + (Zr - Zi)
MOVSS XMM3, XMM5 ' XMM[3] = Sr
SUBSS XMM3, XMM4 ' XMM[3] = Sr - Si
ADDSS XMM3, XMM7 ' XMM[3] = Sr - Si + Cr
REM --------------------------- ' C = C + 1
FADD ST(2) '
REM --------------------------- ' If (C <> Dwell) Then L3
FCOMI ST, ST(1) '
JNE L3 '
REM --------------------------- ' A = True
FXCH ST(3) ' ST[0] = A
FSUB ST(2) ' ST[0] = False
FXCH ST(3) ' ST[3] = A
REM --------------------------- ' Empty FPU Stack Registers
L4: FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
Rem --------------------------- ' Return A
FISTP DWORD PTR [Function] '
END ASM
Rem // Procedural Note
Rem
Rem This procedure, using the first Supplemental Streaming Extension
Rem set, determines if a point exists in or near the Mandelbrot Set
Rem by iterating a complex number for a predetermined number of steps.
Rem
Rem -- Procedural Variables
Rem
Rem A - Accumulator
Rem B - One
Rem C - Iteration
Rem D - Dwell
Rem E - Escape
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [11] In SIMD mode the general purpose registers EAX, EBX, ECX,
Rem and EDX store instructions for SSE1 operations, requiring
Rem the use of the floating point unit for both iteration and
Rem membership state.
Rem
Rem [12] Use of single floating points in this procedure limits
Rem the resolution of the set to one in eight million parts over
Rem the range of -2.0 to +2.0 in both imaginary and real dimensions.
End Function
Rem // Mandelbrot Set Membership Calculation - XMM Double Version
Function Calculate_Mandelbrot_XMMD(ByVal I As Integer, ByVal J As Integer, ByVal W As Integer, ByVal H As Integer) As Integer
Dim As Double Zero = 0
Dim As Double Ci = (J / H) * 2.0 - 1.0
Dim As Double Cr = (I / W) * 2.0 - 1.5
Dim As Double E = Escape
Dim As Integer D = Dwell
ASM
REM --------------------------- ' Load FPU Stack Registers ' [13]
FLDZ ' ST[3] A = 0
FLD1 ' ST[2] B = 1
FILD DWORD PTR [D] ' ST[1] D = Dwell
FLDZ ' ST[0] C = 0
REM --------------------------- ' Load FPU Stack Registers ' [14]
MOVSD XMM7, QWORD [Cr] ' XMM[7] = Cr
MOVSD XMM6, QWORD [Ci] ' XMM[6] = Ci
MOVSD XMM5, QWORD [Zero] ' XMM[5] = Sr
MOVSD XMM4, QWORD [Zero] ' XMM[4] = Si
MOVSD XMM3, QWORD [Cr] ' XMM[3] = Zr
MOVSD XMM2, QWORD [Ci] ' XMM[2] = Zi
MOVSD XMM1, QWORD [E] ' XMM[1] = Escape
MOVSD XMM0, QWORD [Zero] ' XMM[0] = Mu
REM --------------------------- ' Sr = (Zr * Zr)
L5: MOVSD XMM5, XMM3 ' XMM[5] = Zr
MULSD XMM5, XMM5 ' XMM[5] = Zr * Zr
REM --------------------------- ' Si = (Zi * Zi)
MOVSD XMM4, XMM2 ' XMM[4] = Zi
MULSD XMM4, XMM4 ' XMM[4] = Zi * Zi
REM --------------------------- ' M = (Sr + Si)
MOVSD XMM0, XMM5 ' XMM[0] = Sr
ADDSD XMM0, XMM4 ' XMM[0] = (Si + Sr)
REM --------------------------- ' If ((Mu > Escape) Then L5
COMISD XMM0, XMM1 '
JA L6 '
REM --------------------------- ' Zi = 2 * (Zi * Zr) + Ci
MULSD XMM2, XMM3 ' XMM[2] = Zi * Zr
ADDSD XMM2, XMM2 ' XMM[2] = 2 * Zi * Zr
ADDSD XMM2, XMM6 ' XMM[2] = 2 * Zi * Zr + Ci
REM --------------------------- ' Zr = Cr + (Zr - Zi)
MOVSD XMM3, XMM5 ' XMM[3] = Sr
SUBSD XMM3, XMM4 ' XMM[3] = Sr - Si
ADDSD XMM3, XMM7 ' XMM[3] = Sr - Si + Cr
REM --------------------------- ' C = C + 1
FADD ST(2) '
REM --------------------------- ' If (C <> Dwell) Then L5
FCOMI ST, ST(1) '
JNE L5 '
REM --------------------------- ' A = True
FXCH ST(3) ' ST[0] = A
FSUB ST(2) ' ST[0] = False
FXCH ST(3) ' ST[3] = A
REM ----------------------------' Empty FPU Stack Registers
L6: FSTP ST(0) '
FSTP ST(0) '
FSTP ST(0) '
Rem --------------------------- ' Return A
FISTP DWORD PTR [Function] '
END ASM
Rem // Procedural Note
Rem
Rem This procedure, using the second Supplemental Streaming Extension
Rem set, determines if a point exists in or near the Mandelbrot Set
Rem by iterating a complex number for a predetermined number of steps.
Rem
Rem -- Procedural Variables
Rem
Rem A - Accumulator
Rem B - One
Rem C - Iteration
Rem D - Dwell
Rem E - Escape
Rem
Rem Dwell - Dwell Limit
Rem Escape - Escape Value
Rem
Rem Ci - Coordinate Imaginary
Rem Cr - Coordinate Real
Rem
Rem Si - Squared Imaginary
Rem Sr - Squared Real
Rem
Rem Zi - Iterated Imaginary
Rem Zr - Iterated Real
Rem
Rem -- Procedural Footnotes
Rem
Rem [13] In SIMD mode the general purpose registers EAX, EBX, ECX,
Rem and EDX store instructions for SSE2 operations requiring
Rem the use of the floating point unit for both iteration and
Rem membership state.
Rem
Rem [14] Use of double floating points in this procedure limits
Rem the resolution of the set to one in four quadrillion parts
Rem over the range of -2.0 to +2.0 in both imaginary and real
Rem dimensions.
End Function
I’m not sure an actual computer programmer would take this approach to the problem, but it was interesting to work with the most basic level of computer programming short of writing everything in ones and zeros.
P.S. Here’s a list of the websites.
The Computer Language Shootout (http://shootout.alioth.debian.org/)
Free Basic (http://www.freebasic.net/)
Intel Software Development Manuals (http://www.intel.com/products/processor/manuals/)