+- QB64 Phoenix Edition (https://staging.qb64phoenix.com)
+-- Forum: Chatting and Socializing (https://staging.qb64phoenix.com/forumdisplay.php?fid=11)
+--- Forum: General Discussion (https://staging.qb64phoenix.com/forumdisplay.php?fid=2)
+--- Thread: Comparison QB64 compiled with gcc optimizations and without (/showthread.php?tid=366)
RE: Comparison QB64 compiled with Ofast and without - Coolman - 05-08-2022
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 FunctionRE: Comparison QB64 compiled with Ofast and without - Coolman - 05-08-2022
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 SubRE: Comparison QB64 compiled with Ofast and without - Coolman - 05-08-2022
(05-07-2022, 11:04 PM)Jack Wrote: -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/Optimize-Options.html
@Jack. interesting. thanks for the link.
RE: Comparison QB64 compiled with Ofast and without - Coolman - 05-08-2022
@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
RE: Comparison QB64 compiled with Ofast and without - Jack - 05-08-2022
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 FunctionRE: Comparison QB64 compiled with Ofast and without - Coolman - 05-08-2022
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.
RE: Comparison QB64 compiled with Ofast and without - Jack - 05-08-2022
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
RE: Comparison QB64 compiled with Ofast and without - ChiaPet - 05-09-2022
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 ""RE: Comparison QB64 compiled with Ofast and without - ChiaPet - 05-09-2022
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.
RE: Comparison QB64 compiled with Ofast and without - Coolman - 05-09-2022
@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