Comparison QB64 compiled with gcc optimizations and without

[Coolman]

if these options exist, it is well to be used. the gain of speed of the programs is very important with freebasic less with qb64. it is strange. it is necessary to make more test to see if it is interesting with qb64. if you have intensive algorithms of calculations, that interests me. here is another example.

after converting the copy of the old forum unfortunately not complete in pdf format. then in text format. i assembled all the files in a single text file of about 200 mo. i have access to a lot of data on qb64. i found a fast sorting algorithm created by Phlashlite.

I lightened the code a bit.

sorting algorithm QuickSort.

4.4x seconds : program compiled with qb64 -Ofast
4.6x seconds : program compiled with original qb64

the gain is negligible. it is not very efficient for sorting algorithms as in the previous example.

Code: (Select All)

_Title "QuickSort Demo"
Screen _NewImage(640, 480, 32)
Const ARRAY_SIZE = 10000000
Dim Shared DemoArray(ARRAY_SIZE)
Print: Print
start = Timer
For i = 0 To ARRAY_SIZE
    DemoArray(i) = Int(Rnd * 255)
Next
Print " generation of an array of"; Str$(ARRAY_SIZE); " elements performed in"; Timer - start; "seconds":
'==============================================================================
Dim Shared ArrayLength
ArrayLength = ARRAY_LENGTH(DemoArray())
Print: Print
start = Timer
QuickSort 0, ARRAY_SIZE - 1, DemoArray()
Print Str$(ARRAY_SIZE); " members sorted in"; Timer - start; "seconds"
'==============================================================================
Sleep

Sub QuickSort (start, finish, array())
    'Straight from the QB64 wiki
    Dim Hi, Lo, Middle
    Hi = finish
    Lo = start
    Middle = array((Lo + Hi) / 2) 'find middle of array
    Do
        Do While array(Lo) < Middle
            Lo = Lo + 1
        Loop
        Do While array(Hi) > Middle
            Hi = Hi - 1
        Loop
        If Lo <= Hi Then
            Swap array(Lo), array(Hi)
            Lo = Lo + 1
            Hi = Hi - 1
        End If
    Loop Until Lo > Hi
    If Hi > start Then Call QuickSort(start, Hi, array())
    If Lo < finish Then Call QuickSort(Lo, finish, array())
End Sub
'______________________________________________________________________________
'Find the largest member of an array set.
Function ARRAY_MAX_VALUE (Array(), ArrayLength)
    For i = 0 To ArrayLength
        If Array(i) > tmx% Then tmx% = Array(i)
    Next
    ARRAY_MAX_VALUE = tmx%
End Function
'______________________________________________________________________________
'Calculate the size of an array.
Function ARRAY_LENGTH (Array())
    ARRAY_LENGTH = UBound(Array) - LBound(Array)
End Function

[Coolman]

I took this code from bplus. I lightened it up a bit.

Alien Trees Reflection - Plasma Mod.

5.9x seconds : program compiled with qb64 -Ofast
9.2x seconds : program compiled with original qb64

the speed gain is confirmed for the graphic controls.

Code: (Select All)

_Title "Alien Trees Reflection - Plasma Mod" 'b+ trans from SB 2022-05-06
Rem trees reflection.bas 2016-02-22 SmallBASIC 0.12.2 [B+=MGA]
'lakeshore demo repurposed with new and improved trees reflected in lake
Randomize Timer
Const xmax = 1024, ymax = 600

Screen _NewImage(xmax, ymax, 32)
_Delay 0.2
_ScreenMove _Middle
Dim Shared As _Unsigned Long qb(15)
Dim Shared pR, pG, pB, cN, dcN

qb(0) = &HFF000000
qb(1) = &HFF000088
qb(2) = &HFF008800
qb(3) = &HFF008888
qb(4) = &HFF880000
qb(5) = &HFF880088
qb(6) = &HFF888800
qb(7) = &HFFCCCCCC
qb(8) = &HFF888888
qb(9) = &HFF0000FF
qb(10) = &HFF00FF00
qb(11) = &HFF00FFFF
qb(12) = &HFFFF0000
qb(13) = &HFFFF00FF
qb(14) = &HFFFFFF00
qb(15) = &HFFFFFFFF

start = Timer
For boucle% = 1 To 500
    For i = 0 To ymax
        Line (0, i)-(xmax, i), _RGB32(70, 60, i / ymax * 160)
    Next
    stars = xmax * ymax * 10 ^ -4
    horizon = .67 * ymax
    For i = 1 To stars 'stars in sky
        PSet (Rnd * xmax, Rnd * horizon), qb(11)
    Next
    stars = stars / 2
    For i = 1 To stars
        fcirc Rnd * xmax, Rnd * horizon, 1, qb(11)
    Next
    stars = stars / 2
    For i = 1 To stars
        fcirc Rnd * xmax, Rnd * horizon, 2, qb(11)
    Next
    For i = .67 * ymax To .8 * ymax
        gc = max(0, 100 - (i - .67 * ymax) * .5)
        Line (0, i)-(xmax, i), _RGB32(gc, gc, gc)
    Next
    resetPlasma
    branch xmax * .6 + Rnd * .3 * xmax, ymax * .75 - .07 * ymax, 6, 90, xmax / 20, 0
    resetPlasma
    branch Rnd * .3 * xmax, ymax * .75 - .05 * ymax, 7, 90, xmax / 18, 0
    resetPlasma
    branch xmax / 2, ymax * .77, 8, 90, xmax / 16, 0

    Line (0, .8 * ymax)-(xmax, .8 * ymax + 1), _RGB32(70, 70, 70), BF
    For y = .8 * ymax To ymax
        For x = 0 To xmax
            yy = .8 * ymax - (y - .8 * ymax) * 4
            PSet (x, y), Point(x, yy)
        Next
    Next
    _Display
Next boucle%
Print Timer - start
End

Sub branch (x, y, startr, angD, lngth, lev)
    ' local x2,y2,dx,dy,bc,i
    Dim bc As _Unsigned Long
    x2 = x + Cos(_D2R(angD)) * lngth
    y2 = y - Sin(_D2R(angD)) * lngth
    dx = (x2 - x) / lngth
    dy = (y2 - y) / lngth
    'bc = _RGB32(10 + lev, 15 + lev, 10)
    For i = 0 To lngth
        fcirc x + dx * i, y + dy * i, startr, changePlasma~&
    Next
    If startr - 1 < 0 Or lev > 11 Or lngth < 5 Then Exit Sub
    lev2 = lev + 1
    branch x2, y2, startr - 1, angD + 10 + 30 * Rnd, .8 * lngth + .2 * Rnd * lngth, lev2
    branch x2, y2, startr - 1, angD - 10 - 30 * Rnd, .8 * lngth + .2 * Rnd * lngth, lev2
End Sub

Sub fcirc (CX As Long, CY As Long, R As Long, C As _Unsigned Long)
    Dim Radius As Long, RadiusError As Long
    Dim X As Long, Y As Long
    Radius = Abs(R): RadiusError = -Radius: X = Radius: Y = 0
    If Radius = 0 Then PSet (CX, CY), C: Exit Sub
    Line (CX - X, CY)-(CX + X, CY), C, BF
    While X > Y
        RadiusError = RadiusError + Y * 2 + 1
        If RadiusError >= 0 Then
            If X <> Y + 1 Then
                Line (CX - Y, CY - X)-(CX + Y, CY - X), C, BF
                Line (CX - Y, CY + X)-(CX + Y, CY + X), C, BF
            End If
            X = X - 1
            RadiusError = RadiusError - X * 2
        End If
        Y = Y + 1
        Line (CX - X, CY - Y)-(CX + X, CY - Y), C, BF
        Line (CX - X, CY + Y)-(CX + X, CY + Y), C, BF
    Wend
End Sub

Function max (x, y)
    If x > y Then max = x Else max = y
End Function

Function changePlasma~& ()
    cN = cN + dcN 'dim shared cN as _Integer64, pR as long, pG as long, pB as long
    changePlasma~& = _RGB32(127 + 127 * Sin(pR * cN), 127 + 127 * Sin(pG * cN), 127 + 127 * Sin(pB * cN))
End Function

Sub resetPlasma ()
    pR = Rnd ^ 2: pG = Rnd ^ 2: pB = Rnd ^ 2: dcN = Rnd
End Sub

[Coolman]

-ffast-math

    Sets the options -fno-math-errno, -funsafe-math-optimizations, -ffinite-math-only, -fno-rounding-math, -fno-signaling-nans, -fcx-limited-range and -fexcess-precision=fast.

    This option causes the preprocessor macro __FAST_MATH__ to be defined.

    This option is not turned on by any -O option besides -Ofast since it can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions. It may, however, yield faster code for programs that do not require the guarantees of these specifications.

https://gcc.gnu.org/onlinedocs/gcc/Optim...tions.html

@Jack. interesting. thanks for the link.

[Coolman]

@Jack

Are you the same person who posted :

suggest adding -O2 to makeline

in the old forum. i can't restore the whole nbody source code.

if it's you, could you post it here so I can test it...

thanks

[Jack]

hello Coolman, yes that's me, here's the code
note about -Ofast, for programs not requiring strict adherence to the IEEE standard especially graphic programs it may work without a problem and possibly much faster

Code: (Select All)

'https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-gcc-1.html
' The Computer Language Benchmarks Game
' https://benchmarksgame-team.pages.debian.net/benchmarksgame/

' contributed by Christoph Bauer
'
' translated to QB64 by Jack
' modified a bit to avoid using pointer to array

declare FUNCTION main# (n AS LONG)

Print Str$(main(50000000)) + " seconds"
End

Const pi = 3.141592653589793
Const solar_mass = (4 * pi) * pi
Const days_per_year = 365.24

Type planet
    x As Double
    y As Double
    z As Double
    vx As Double
    vy As Double
    vz As Double
    mass As Double
End Type

Sub advance (nbodies As Long, bodies() As planet, dt As Double)
    Dim i As Long, j As Long
    Dim dx As Double, dy As Double, dz As Double, distance As Double, mag As Double
    For i = 0 To nbodies - 1
        For j = i + 1 To nbodies - 1
            dx = bodies(i).x - bodies(j).x
            dy = bodies(i).y - bodies(j).y
            dz = bodies(i).z - bodies(j).z
            distance = Sqr(((dx * dx) + (dy * dy)) + (dz * dz))
            mag = dt / ((distance * distance) * distance)
            bodies(i).vx = bodies(i).vx - (dx * bodies(j).mass) * mag
            bodies(i).vy = bodies(i).vy - (dy * bodies(j).mass) * mag
            bodies(i).vz = bodies(i).vz - (dz * bodies(j).mass) * mag
            bodies(j).vx = bodies(j).vx + (dx * bodies(i).mass) * mag
            bodies(j).vy = bodies(j).vy + (dy * bodies(i).mass) * mag
            bodies(j).vz = bodies(j).vz + (dz * bodies(i).mass) * mag
        Next
    Next
    For i = 0 To nbodies - 1
        bodies(i).x = bodies(i).x + dt * bodies(i).vx
        bodies(i).y = bodies(i).y + dt * bodies(i).vy
        bodies(i).z = bodies(i).z + dt * bodies(i).vz
    Next
End Sub

Function energy# (nbodies As Long, bodies() As planet)
    Dim i As Long, j As Long
    Dim e As Double, dx As Double, dy As Double, dz As Double, distance As Double, mag As Double
    e = 0.0
    For i = 0 To nbodies - 1
        e = e + (0.5 * bodies(i).mass) * (((bodies(i).vx * bodies(i).vx) + (bodies(i).vy * bodies(i).vy)) + (bodies(i).vz * bodies(i).vz))
        For j = i + 1 To nbodies - 1
            dx = bodies(i).x - bodies(j).x
            dy = bodies(i).y - bodies(j).y
            dz = bodies(i).z - bodies(j).z
            distance = Sqr(((dx * dx) + (dy * dy)) + (dz * dz))
            e = e - (bodies(i).mass * bodies(j).mass) / distance
        Next
    Next
    energy# = e
End Function

Sub offset_momentum (nbody As Long, bodies() As planet)
    Dim px As Double, py As Double, pz As Double
    Dim I As Long
    For I = 0 To nbody - 1
        px = px + bodies(I).vx * bodies(I).mass
        py = py + bodies(I).vy * bodies(I).mass
        pz = pz + bodies(I).vz * bodies(I).mass
    Next
    bodies(0).vx = (-px) / ((4 * pi) * pi)
    bodies(0).vy = (-py) / ((4 * pi) * pi)
    bodies(0).vz = (-pz) / ((4 * pi) * pi)
End Sub


Function main# (n As Long)
    Const NBODIES = 5
    Dim bodies(0 To 4) As planet
    bodies(0).x = 0: bodies(0).y = 0: bodies(0).z = 0: bodies(0).vx = 0: bodies(0).vy = 0: bodies(0).vz = 0: bodies(0).mass = (4 * pi) * pi
    bodies(1).x = 4.84143144246472090D+00: bodies(1).y = -1.16032004402742839D+00: bodies(1).z = -1.03622044471123109D-01: bodies(1).vx = 1.66007664274403694D-03 * days_per_year: bodies(1).vy = 7.69901118419740425D-03 * days_per_year: bodies(1).vz = (-6.90460016972063023D-05) * days_per_year: bodies(1).mass = 9.54791938424326609D-04 * ((4 * pi) * pi)
    bodies(2).x = 8.34336671824457987D+00: bodies(2).y = 4.12479856412430479D+00: bodies(2).z = -4.03523417114321381D-01: bodies(2).vx = (-2.76742510726862411D-03) * days_per_year: bodies(2).vy = 4.99852801234917238D-03 * days_per_year: bodies(2).vz = 2.30417297573763929D-05 * days_per_year: bodies(2).mass = 2.85885980666130812D-04 * ((4 * pi) * pi)
    bodies(3).x = 1.28943695621391310D+01: bodies(3).y = -1.51111514016986312D+01: bodies(3).z = -2.23307578892655734D-01: bodies(3).vx = 2.96460137564761618D-03 * days_per_year: bodies(3).vy = 2.37847173959480950D-03 * days_per_year: bodies(3).vz = (-2.96589568540237556D-05) * days_per_year: bodies(3).mass = 4.36624404335156298D-05 * ((4 * pi) * pi)
    bodies(4).x = 1.53796971148509165D+01: bodies(4).y = -2.59193146099879641D+01: bodies(4).z = 1.79258772950371181D-01: bodies(4).vx = 2.68067772490389322D-03 * days_per_year: bodies(4).vy = 1.62824170038242295D-03 * days_per_year: bodies(4).vz = (-9.51592254519715870D-05) * days_per_year: bodies(4).mass = 5.15138902046611451D-05 * ((4 * pi) * pi)

    Dim i As Long
    Dim t As Double
    t = Timer
    Call offset_momentum(NBODIES, bodies())
    Print Using "##.#########"; energy#(NBODIES, bodies())
    For i = 1 To n
        Call advance(NBODIES, bodies(), 0.01)
    Next
    Print Using "##.#########"; energy#(NBODIES, bodies())
    main# = Timer - t
End Function

[Coolman]

66.1x seconds : program compiled with qb64 -Ofast
68.5x seconds : program compiled with original qb64

2 seconds apart. that's not good and surprising. you said in your post in the old forum :

time before adding -O2 to makeline_osx.txt 138.84375 seconds, after adding -O2, 42.6875 seconds, that's 3.25 times faster.

if this is what you said. it means that the -O2 option is more efficient than -Ofast for qb64 because with freebasic it is the opposite.

the tests performed on the sorting algorithms are equally surprising and unexpected.

I don't have time today but I'll have to modify my script to compile qb64 totally with the -O2 option to do this test again.

@Jack. Thank you for posting your code.

[Jack]

my times
48.02 seconds : program compiled with original qb64
18.13 seconds : program compiled with -O2
16.60 seconds : program compiled with -O3
17.25 seconds : program compiled with -Ofast
BTW, I make it habit to change the program source by adding and deleting a space to ensure that the source will be recompiled

[ChiaPet]

Code: (Select All)

' Speed Test for g++ optimization parameters
' modifies .\internal\c\makeline_win.txt (restores changes on exit)

o$ = "c_compiler\bin\g++ -s -Wfatal-errors -w -Wall qbx.cpp -lws2_32 -lwinspool parts\core\os\win\src.a "
o$ = o$ + "-lopengl32 -lglu32 -lwinmm -lgdi32 -mwindows -static-libgcc -static-libstdc++ -D GLEW_STATIC "
o$ = o$ + "-D FREEGLUT_STATIC -lksguid -lole32 -lwinmm -ldxguid -o ..\..\"

Shell "copy .\internal\c\makeline_win.txt .\internal\c\makeline_win.bak"

Open "sptest_t1.bas" For Output As #1
Do
    Read c$
    If c$ = "" Then Exit Do
    Print #1, c$
Loop
Close #1
Shell _Hide "copy sptest_t1.bas sptest_t2.bas"
Shell _Hide "copy sptest_t1.bas sptest_t3.bas"
Shell _Hide "copy sptest_t1.bas sptest_t4.bas"

m$(1) = "" '          no optimization
m$(2) = "-O2" '       oxygen
m$(3) = "-O3" '       ozone
m$(4) = "-Ofast" '    Lamborghini?

start_at = 2 '        skip NO optimization because it's a turkey

For pass = start_at To 4
    c$ = Left$(o$, 22) + m$(pass) + " " + Right$(o$, Len(o$) - 22)
    Open ".\internal\c\makeline_win.txt" For Output As #1
    Print #1, c$
    Close #1
    c$ = "qb64 -x sptest_t" + Chr$(48 + pass)
    Shell _Hide c$
    Shell _Hide "copy .\internal\c\makeline_win.bak .\internal\c\makeline_win.txt" '    restore file
    GoSub keycheck '                                                                    exit?
Next pass

Do
    best# = 99999
    For pass = start_at To 4
        start# = Timer '
        c$ = "sptest_t" + Chr$(48 + pass)
        Shell c$ '                                                                      run the test
        elapse# = Timer - start#
        If elapse# < best# Then best# = elapse#: bestpass = pass
        Print Using "####.##  "; elapse#;: Print m$(pass) '                             show the timing
        GoSub keycheck '
    Next pass
    count(bestpass) = count(bestpass) + 1
    For pass = start_at To 4
        Print Using "####  "; count(pass);
    Next pass
    Print Using "####  "; count(1) + count(2) + count(3) + count(4);
    Print
Loop

keycheck:
If _Exit Or Len(InKey$) Then
    Shell _Hide "del .\internal\c\makeline_win.bak"
    Shell _Hide "del sptest_t*.*"
    System
End If
Return

Data "Dim c(256) as _unsigned long"
Data "xm = 100: ym = 100"
Data "Screen _NewImage(xm, ym, 32)"
Data "_screenmove _desktopwidth-xm,0"
Data "x1 = -2: x2 = .6: y1 = -1.25: y2 = 1.25"
Data "xd = (x2 - x1) / xm: yd = (y2 - y1) / ym"
Data "c(0) = _RGB32(0, 0, 0)"
Data "For i = 1 To 31"
Data "c(i) = _RGB32(rnd*255,rnd*255,rnd*255)"
Data "Next i"
Data "For y = 0 To ym - 1"
Data "mandely = y1 + y * yd"
Data "For x = 0 To xm - 1"
Data "mandelx = x1 + x * xd"
Data "Itera = 5000: real = 0: imag = 0: size = -1"
Data "While (Itera > 0) And (size < 0)"
Data "Itera = Itera - 1"
Data "hold = imag"
Data "imag = (real * imag) * 2 + mandely"
Data "real = real * real - hold * hold + mandelx"
Data "size = (real * real + imag * imag) - 99"
Data "Wend"
Data "PSet (x, y), c(Itera Mod 32)"
Data "Next x"
Data "Next y"
Data "System"
Data ""

[ChiaPet]

A summary of the above program, since I don't know how to make it all one message:

1. Creates copies of a test program (a Mandelbrot).
2. Compiles each copy with a different optimization parameter.
3. Endless loop of running the copies, incrementing a counter for each winner, showing the timing and counters.

In 100 runs, -O2 won 44 times, -O3 37 times, and -Ofast 19. 

I suppose the results will vary according to processor, O/S, and what test program is used.

[Coolman]

@Jack. thank you for posting your results.

@ChiaPet. thank you for posting your code. while reviewing it, i realized that i forgot in my script to modify the files makeline_lnx.txt and makeline_lnx_nogui.txt...

after reading some information about gcc, i decided to use the -O3 compiler option. it's a good compromise between security and speed.

I have redone the tests. here are the results :

simple code using pset.

2.3x seconds : program compiled with qb64 -O3
3.3x seconds : program compiled with original qb64

Fractal Tree : I got this code from the site rosettacode.

2.9x seconds : program compiled with qb64 -O3
5.0x seconds : program compiled with original qb64

Bucketsort : I got this code from the site rosettacode.

1.1x seconds : program compiled with qb64 -O3
2.8x seconds : program compiled with original qb64

sorting algorithm QuickSort.

3.0x seconds : program compiled with qb64 -O3
4.6x seconds : program compiled with original qb64

Alien Trees Reflection - Plasma Mod.

5.3x seconds : program compiled with qb64 -O3
9.2x seconds : program compiled with original qb64

Nbody

16.9x seconds : program compiled with qb64 -O3
68.0x seconds : program compiled with original qb64

it becomes very interesting...

Code: (Select All)

'https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-gcc-1.html
' The Computer Language Benchmarks Game
' https://benchmarksgame-team.pages.debian.net/benchmarksgame/

' contributed by Christoph Bauer
'
' translated to QB64 by Jack
' modified a bit to avoid using pointer to array

declare FUNCTION main# (n AS LONG)

Print Str$(main(50000000)) + " seconds"
End

Const pi = 3.141592653589793
Const solar_mass = (4 * pi) * pi
Const days_per_year = 365.24

Type planet
    x As Double
    y As Double
    z As Double
    vx As Double
    vy As Double
    vz As Double
    mass As Double
End Type

Sub advance (nbodies As Long, bodies() As planet, dt As Double)
    Dim i As Long, j As Long
    Dim dx As Double, dy As Double, dz As Double, distance As Double, mag As Double
    For i = 0 To nbodies - 1
        For j = i + 1 To nbodies - 1
            dx = bodies(i).x - bodies(j).x
            dy = bodies(i).y - bodies(j).y
            dz = bodies(i).z - bodies(j).z
            distance = Sqr(((dx * dx) + (dy * dy)) + (dz * dz))
            mag = dt / ((distance * distance) * distance)
            bodies(i).vx = bodies(i).vx - (dx * bodies(j).mass) * mag
            bodies(i).vy = bodies(i).vy - (dy * bodies(j).mass) * mag
            bodies(i).vz = bodies(i).vz - (dz * bodies(j).mass) * mag
            bodies(j).vx = bodies(j).vx + (dx * bodies(i).mass) * mag
            bodies(j).vy = bodies(j).vy + (dy * bodies(i).mass) * mag
            bodies(j).vz = bodies(j).vz + (dz * bodies(i).mass) * mag
        Next
    Next
    For i = 0 To nbodies - 1
        bodies(i).x = bodies(i).x + dt * bodies(i).vx
        bodies(i).y = bodies(i).y + dt * bodies(i).vy
        bodies(i).z = bodies(i).z + dt * bodies(i).vz
    Next
End Sub

Function energy# (nbodies As Long, bodies() As planet)
    Dim i As Long, j As Long
    Dim e As Double, dx As Double, dy As Double, dz As Double, distance As Double, mag As Double
    e = 0.0
    For i = 0 To nbodies - 1
        e = e + (0.5 * bodies(i).mass) * (((bodies(i).vx * bodies(i).vx) + (bodies(i).vy * bodies(i).vy)) + (bodies(i).vz * bodies(i).vz))
        For j = i + 1 To nbodies - 1
            dx = bodies(i).x - bodies(j).x
            dy = bodies(i).y - bodies(j).y
            dz = bodies(i).z - bodies(j).z
            distance = Sqr(((dx * dx) + (dy * dy)) + (dz * dz))
            e = e - (bodies(i).mass * bodies(j).mass) / distance
        Next
    Next
    energy# = e
End Function

Sub offset_momentum (nbody As Long, bodies() As planet)
    Dim px As Double, py As Double, pz As Double
    Dim I As Long
    For I = 0 To nbody - 1
        px = px + bodies(I).vx * bodies(I).mass
        py = py + bodies(I).vy * bodies(I).mass
        pz = pz + bodies(I).vz * bodies(I).mass
    Next
    bodies(0).vx = (-px) / ((4 * pi) * pi)
    bodies(0).vy = (-py) / ((4 * pi) * pi)
    bodies(0).vz = (-pz) / ((4 * pi) * pi)
End Sub


Function main# (n As Long)
    Const NBODIES = 5
    Dim bodies(0 To 4) As planet
    bodies(0).x = 0: bodies(0).y = 0: bodies(0).z = 0: bodies(0).vx = 0: bodies(0).vy = 0: bodies(0).vz = 0: bodies(0).mass = (4 * pi) * pi
    bodies(1).x = 4.84143144246472090D+00: bodies(1).y = -1.16032004402742839D+00: bodies(1).z = -1.03622044471123109D-01: bodies(1).vx = 1.66007664274403694D-03 * days_per_year: bodies(1).vy = 7.69901118419740425D-03 * days_per_year: bodies(1).vz = (-6.90460016972063023D-05) * days_per_year: bodies(1).mass = 9.54791938424326609D-04 * ((4 * pi) * pi)
    bodies(2).x = 8.34336671824457987D+00: bodies(2).y = 4.12479856412430479D+00: bodies(2).z = -4.03523417114321381D-01: bodies(2).vx = (-2.76742510726862411D-03) * days_per_year: bodies(2).vy = 4.99852801234917238D-03 * days_per_year: bodies(2).vz = 2.30417297573763929D-05 * days_per_year: bodies(2).mass = 2.85885980666130812D-04 * ((4 * pi) * pi)
    bodies(3).x = 1.28943695621391310D+01: bodies(3).y = -1.51111514016986312D+01: bodies(3).z = -2.23307578892655734D-01: bodies(3).vx = 2.96460137564761618D-03 * days_per_year: bodies(3).vy = 2.37847173959480950D-03 * days_per_year: bodies(3).vz = (-2.96589568540237556D-05) * days_per_year: bodies(3).mass = 4.36624404335156298D-05 * ((4 * pi) * pi)
    bodies(4).x = 1.53796971148509165D+01: bodies(4).y = -2.59193146099879641D+01: bodies(4).z = 1.79258772950371181D-01: bodies(4).vx = 2.68067772490389322D-03 * days_per_year: bodies(4).vy = 1.62824170038242295D-03 * days_per_year: bodies(4).vz = (-9.51592254519715870D-05) * days_per_year: bodies(4).mass = 5.15138902046611451D-05 * ((4 * pi) * pi)

    Dim i As Long
    Dim t As Double
    t = Timer
    Call offset_momentum(NBODIES, bodies())
    Print Using "##.#########"; energy#(NBODIES, bodies())
    For i = 1 To n
        Call advance(NBODIES, bodies(), 0.01)
    Next
    Print Using "##.#########"; energy#(NBODIES, bodies())
    main# = Timer - t
End Function