[MLton] NBody benchmark

Christoph Bauer ich@christoph-bauer.net
Wed, 16 Feb 2005 09:25:08 +0100


--=-=-=

Hi, 

here is the nbody-benchmark for mlton. It is a 1:1 translation of the
ocaml nbody code written by Christophe Troestler. So maybe
the code could be improved. 

MLTons result is (slightly) faster than ocamls and slower than the gcc
version. BTW in opposite to the two other compilers mlton doesn't use
strict IEEE floation point as default.

Christoph Bauer


--=-=-=
Content-Type: application/octet-stream
Content-Disposition: attachment; filename=nbody.sml
Content-Transfer-Encoding: base64

dmFsIHBpID0gMy4xNDE1OTI2NTM1ODk3OTMKdmFsIHNvbGFyX21hc3MgPSA0LjAgKiBwaSAqIHBp
CnZhbCBkYXlzX3Blcl95ZWFyID0gMzY1LjI0Cgp0eXBlIHBsYW5ldCA9IHsKICB4IDogcmVhbCBy
ZWYsCiAgeSA6IHJlYWwgcmVmLAogIHogOiByZWFsIHJlZiwKICB2eDogcmVhbCByZWYsICAKICB2
eTogcmVhbCByZWYsICAKICB2ejogcmVhbCByZWYsCiAgbWFzcyA6IHJlYWwKfQoKCmZ1biBmb3Ig
KHN0YXJ0LCBzdG9wLCBmKSA9CiAgIGxldAogICAgICBmdW4gbG9vcCBpID0KICAgICAgICAgaWYg
aSA+IHN0b3AKICAgICAgICAgICAgdGhlbiAoKQogICAgICAgICBlbHNlIChmIGk7IGxvb3AgKGkg
KyAxKSkKICAgaW4KICAgICAgbG9vcCBzdGFydAogICBlbmQKCmZ1biBhZHZhbmNlIChib2RpZXMg
OiBwbGFuZXQgYXJyYXkpIGR0ID0KICBsZXQgCiAgICAgdmFsIG4gPSBBcnJheS5sZW5ndGggYm9k
aWVzIC0gMQogIGluIAogICBmb3IgKDAsIG4sIGZuIGkgPT4KICAgbGV0IAogICAgICB2YWwgYiA9
IEFycmF5LnN1YiAoYm9kaWVzLCBpKSAKICAgaW4KICAgICBmb3IgKGkrMSwgbiwgZm4gaiA9Pgog
ICAgICAgICBsZXQgCiAgICAgICAgICAgdmFsIGInID0gQXJyYXkuc3ViIChib2RpZXMsIGopCiAg
ICAgICAgICAgdmFsIGR4ID0gISgjeCBiKSAtICEoI3ggYicpICAKICAgICAgICAgICB2YWwgZHkg
PSAhKCN5IGIpIC0gISgjeSBiJykKICAgICAgICAgICB2YWwgZHogPSAhKCN6IGIpIC0gISgjeiBi
JykKICAgICAgICAgICB2YWwgZGlzdGFuY2UgPSBNYXRoLnNxcnQoZHggKiBkeCArIGR5ICogZHkg
KyBkeiAqIGR6KSAKICAgICAgICAgICB2YWwgbWFnID0gZHQgLyAoZGlzdGFuY2UgKiBkaXN0YW5j
ZSAqIGRpc3RhbmNlKSAKICAgICAgICAgaW4KICAgICAgICAgICAjdnggYiA6PSAhKCN2eCBiKSAt
IGR4ICogKCNtYXNzIGInKSAqIG1hZzsKICAgICAgICAgICAjdnkgYiA6PSAhKCN2eSBiKSAtIGR5
ICogKCNtYXNzIGInKSAqIG1hZzsKICAgICAgICAgICAjdnogYiA6PSAhKCN2eiBiKSAtIGR6ICog
KCNtYXNzIGInKSAqIG1hZzsKCSAgIAogICAgICAgICAgICN2eCBiJyA6PSAhKCN2eCBiJykgKyBk
eCAqICgjbWFzcyBiKSAqIG1hZzsKICAgICAgICAgICAjdnkgYicgOj0gISgjdnkgYicpICsgZHkg
KiAoI21hc3MgYikgKiBtYWc7CiAgICAgICAgICAgI3Z6IGInIDo9ICEoI3Z6IGInKSArIGR6ICog
KCNtYXNzIGIpICogbWFnCiAgICAgICAgZW5kKQogIGVuZCk7CgogIGZvciggMCwgbiwgZm4gaSA9
PgogICAgbGV0IAogICAgICAgdmFsIGIgPSBBcnJheS5zdWIgKGJvZGllcywgaSkgCiAgICBpbgoJ
I3ggYiA6PSAhKCN4IGIpICsgZHQgKiAhKCN2eCBiKTsKCSN5IGIgOj0gISgjeSBiKSArIGR0ICog
ISgjdnkgYik7CgkjeiBiIDo9ICEoI3ogYikgKyBkdCAqICEoI3Z6IGIpCiAgICBlbmQpCmVuZAoK
ZnVuIGVuZXJneSAoYm9kaWVzIDogcGxhbmV0IGFycmF5KSA9CiAgbGV0IAogICAgIHZhbCBlID0g
cmVmIDAuMCAKICBpbgogICAgZm9yKDAsIEFycmF5Lmxlbmd0aCBib2RpZXMgLSAxLCBmbiBpID0+
CiAgICBsZXQgCiAgICAgICAgdmFsIGIgPSBBcnJheS5zdWIgKGJvZGllcywgaSkgCiAgICBpbgog
ICAgICAgZSA6PSAhZSArIDAuNSAqICgjbWFzcyBiKSAqIAoJCSAoISgjdnggYikgKiAhKCN2eCBi
KSArICEoI3Z5IGIpIAoJCSAgKiAhKCN2eSBiKSArICEoI3Z6IGIpICogISgjdnogYikpOwogICAg
ICAgZm9yIChpKzEsIEFycmF5Lmxlbmd0aCBib2RpZXMgLSAxLCBmbiBqID0+IAogICAgICAgICBs
ZXQgCiAgICAgICAgICAgIHZhbCBiJyA9IEFycmF5LnN1YiAoYm9kaWVzLCBqKQogICAgICAgICAg
ICB2YWwgZHggPSAhKCN4IGIpIC0gISgjeCBiJykgCiAgICAgICAgICAgIHZhbCBkeSA9ICEoI3kg
YikgLSAhKCN5IGInKSAgCiAgICAgICAgICAgIHZhbCBkeiA9ICEoI3ogYikgLSAhKCN6IGInKQog
ICAgICAgICAgICB2YWwgZGlzdGFuY2UgPSBNYXRoLnNxcnQoZHggKiBkeCArIGR5ICogZHkgKyBk
eiAqIGR6KSAKICAgICAgICAgaW4KICAgICAgICAgICAgIGUgOj0gIWUgLSAoKCNtYXNzIGIpICog
KCNtYXNzIGInKSkgLyBkaXN0YW5jZQogICAgICAgICBlbmQpCiAgIGVuZCk7CiAgICFlCiAgZW5k
CgoKCmZ1biBvZmZzZXRfbW9tZW50dW0gKGJvZGllcyA6IHBsYW5ldCBhcnJheSkgPQogIGxldCAK
ICAgICAgdmFsIHB4ID0gcmVmIDAuMAogICAgICB2YWwgcHkgPSByZWYgMC4wIAogICAgICB2YWwg
cHogPSByZWYgMC4wCiAgIGluCiAgICAgZm9yICgwLCBBcnJheS5sZW5ndGggYm9kaWVzIC0gMSwg
CgkgIGZuIGkgPT4KCSAgICAgKHB4IDo9ICFweCArICEoI3Z4IChBcnJheS5zdWIgKGJvZGllcywg
aSkpKSAqICgjbWFzcyAoQXJyYXkuc3ViIChib2RpZXMsIGkpKSk7CgkgICAgICBweSA6PSAhcHkg
KyAhKCN2eSAoQXJyYXkuc3ViIChib2RpZXMsIGkpKSkgKiAoI21hc3MgKEFycmF5LnN1YiAoYm9k
aWVzLCBpKSkpOwoJICAgICAgcHogOj0gIXB6ICsgISgjdnogKEFycmF5LnN1YiAoYm9kaWVzLCBp
KSkpICogKCNtYXNzIChBcnJheS5zdWIgKGJvZGllcywgaSkpKSkpOwogICAgICN2eCAoQXJyYXku
c3ViIChib2RpZXMsIDApKSA6PSB+ICghcHggLyBzb2xhcl9tYXNzKTsKICAgICAjdnkgKEFycmF5
LnN1YiAoYm9kaWVzLCAwKSkgOj0gfiAoIXB5IC8gc29sYXJfbWFzcyk7CiAgICAgI3Z6IChBcnJh
eS5zdWIgKGJvZGllcywgMCkpIDo9IH4gKCFweiAvIHNvbGFyX21hc3MpCiAgZW5kCgp2YWwganVw
aXRlciA9IHsKICB4ID0gcmVmIDQuODQxNDMxNDQyNDY0NzIwOTAsCiAgeSA9IHJlZiB+MS4xNjAz
MjAwNDQwMjc0MjgzOSwKICB6ID0gcmVmIH4xLjAzNjIyMDQ0NDcxMTIzMTA5ZX4xLAogIHZ4ID0g
cmVmICgxLjY2MDA3NjY0Mjc0NDAzNjk0ZX4zICogZGF5c19wZXJfeWVhciksCiAgdnkgPSByZWYg
KDcuNjk5MDExMTg0MTk3NDA0MjVlfjMgKiBkYXlzX3Blcl95ZWFyKSwKICB2eiA9IHJlZiAofjYu
OTA0NjAwMTY5NzIwNjMwMjNlfjUgKiBkYXlzX3Blcl95ZWFyKSwKICBtYXNzID0gOS41NDc5MTkz
ODQyNDMyNjYwOWV+NCAqIHNvbGFyX21hc3MKfQoKdmFsIHNhdHVybiA9IHsKICB4ID0gcmVmIDgu
MzQzMzY2NzE4MjQ0NTc5ODcsCiAgeSA9IHJlZiA0LjEyNDc5ODU2NDEyNDMwNDc5ZTAwLAogIHog
PSByZWYgfjQuMDM1MjM0MTcxMTQzMjEzODFlfjAxLAogIHZ4ID0gcmVmICh+Mi43Njc0MjUxMDcy
Njg2MjQxMWV+MDMgKiBkYXlzX3Blcl95ZWFyKSwKICB2eSA9IHJlZiAoNC45OTg1MjgwMTIzNDkx
NzIzOGV+MDMgKiBkYXlzX3Blcl95ZWFyKSwKICB2eiA9IHJlZiAoMi4zMDQxNzI5NzU3Mzc2Mzky
OWV+MDUgKiBkYXlzX3Blcl95ZWFyKSwKICBtYXNzID0gMi44NTg4NTk4MDY2NjEzMDgxMmV+MDQg
KiBzb2xhcl9tYXNzCn0KCnZhbCB1cmFudXMgPSB7CiAgeCA9IHJlZiAxLjI4OTQzNjk1NjIxMzkx
MzEwZTAxLAogIHkgPSByZWYgfjEuNTExMTE1MTQwMTY5ODYzMTJlMDEsCiAgeiA9IHJlZiB+Mi4y
MzMwNzU3ODg5MjY1NTczNGV+MDEsCiAgdnggPSByZWYgKDIuOTY0NjAxMzc1NjQ3NjE2MThlfjAz
ICogZGF5c19wZXJfeWVhciksCiAgdnkgPSByZWYgKDIuMzc4NDcxNzM5NTk0ODA5NTBlfjAzICog
ZGF5c19wZXJfeWVhciksCiAgdnogPSByZWYgKH4yLjk2NTg5NTY4NTQwMjM3NTU2ZX4wNSAqIGRh
eXNfcGVyX3llYXIpLAogIG1hc3MgPSA0LjM2NjI0NDA0MzM1MTU2Mjk4ZX4wNSAqIHNvbGFyX21h
c3MKfQoKdmFsIG5lcHR1bmUgPSB7CiAgeCA9IHJlZiAxLjUzNzk2OTcxMTQ4NTA5MTY1ZTAxLAog
IHkgPSByZWYgfjIuNTkxOTMxNDYwOTk4Nzk2NDFlMDEsCiAgeiA9IHJlZiAxLjc5MjU4NzcyOTUw
MzcxMTgxZX4wMSwKICB2eCA9IHJlZiAoMi42ODA2Nzc3MjQ5MDM4OTMyMmV+MDMgKiBkYXlzX3Bl
cl95ZWFyKSwKICB2eSA9IHJlZiAoMS42MjgyNDE3MDAzODI0MjI5NWV+MDMgKiBkYXlzX3Blcl95
ZWFyKSwKICB2eiA9IHJlZiAofjkuNTE1OTIyNTQ1MTk3MTU4NzBlfjA1ICogZGF5c19wZXJfeWVh
ciksCiAgbWFzcyA9IDUuMTUxMzg5MDIwNDY2MTE0NTFlfjA1ICogc29sYXJfbWFzcwp9Cgp2YWwg
c3VuID0gewogIHggPSByZWYgMC4wLCAgeSA9IHJlZiAwLjAsICB6ID0gcmVmIDAuMCwgIHZ4ID0g
cmVmIDAuMCwgIHZ5ID0gcmVmIDAuMCwgdnogPSByZWYgMC4wLAogIG1hc3M9IHNvbGFyX21hc3MK
fQoKCgpmdW4gcHJpbnRyIHIgPSAKICBsZXQgCiAgICB2YWwgcyA9IFJlYWwuZm10IChTdHJpbmdD
dnQuRklYIChTT01FIDkpKSByCiAgaW4gKHByaW50IHM7IHByaW50ICJcbiIpIGVuZAoKdmFsIGJv
ZGllcyA9IEFycmF5LmZyb21MaXN0IFsgc3VuLCBqdXBpdGVyLCBzYXR1cm4sIHVyYW51cywgbmVw
dHVuZSBdCgoKZnVuIG1haW4gYXJncyA9CiAgICBsZXQgCgl2YWwgbiA9IAoJICAgIGNhc2UgSW50
LmZyb21TdHJpbmcgKExpc3QuaGQgYXJncykgb2YKCQlTT01FIGkgPT4gaQoJICAgICAgfCBOT05F
ID0+IDAKCQkJCiAgICBpbgoJb2Zmc2V0X21vbWVudHVtIGJvZGllczsKCXByaW50ciAoZW5lcmd5
IGJvZGllcyk7Cglmb3IgKDEsIG4sIGZuIF8gPT4gYWR2YW5jZSBib2RpZXMgMC4wMSk7Cglwcmlu
dHIgKGVuZXJneSBib2RpZXMpCiAgICBlbmQKCnZhbCBfID0gbWFpbiAoQ29tbWFuZExpbmUuYXJn
dW1lbnRzICgpKQo=
--=-=-=
Content-Type: application/octet-stream
Content-Disposition: attachment; filename=nbody.ml
Content-Transfer-Encoding: base64

bGV0IHBpID0gMy4xNDE1OTI2NTM1ODk3OTMKbGV0IHNvbGFyX21hc3MgPSA0LiAqLiBwaSAqLiBw
aQpsZXQgZGF5c19wZXJfeWVhciA9IDM2NS4yNAoKdHlwZSBwbGFuZXQgPSB7CiAgbXV0YWJsZSB4
IDogZmxvYXQ7ICBtdXRhYmxlIHkgOiBmbG9hdDsgIG11dGFibGUgeiA6IGZsb2F0OwogIG11dGFi
bGUgdng6IGZsb2F0OyAgbXV0YWJsZSB2eTogZmxvYXQ7ICBtdXRhYmxlIHZ6OiBmbG9hdDsKICBt
YXNzIDogZmxvYXQ7Cn0KCgpsZXQgYWR2YW5jZSBib2RpZXMgZHQgPQogIGxldCBuID0gQXJyYXku
bGVuZ3RoIGJvZGllcyAtIDEgaW4KICBmb3IgaSA9IDAgdG8gbiBkbwogICAgbGV0IGIgPSBib2Rp
ZXMuKGkpIGluCiAgICBmb3IgaiA9IGkrMSB0byBuIGRvCiAgICAgIGxldCBiJyA9IGJvZGllcy4o
aikgaW4KICAgICAgbGV0IGR4ID0gYi54IC0uIGInLnggIGFuZCBkeSA9IGIueSAtLiBiJy55ICBh
bmQgZHogPSBiLnogLS4gYicueiBpbgogICAgICBsZXQgZGlzdGFuY2UgPSBzcXJ0KGR4ICouIGR4
ICsuIGR5ICouIGR5ICsuIGR6ICouIGR6KSBpbgogICAgICBsZXQgbWFnID0gZHQgLy4gKGRpc3Rh
bmNlICouIGRpc3RhbmNlICouIGRpc3RhbmNlKSBpbgoKICAgICAgYi52eCA8LSBiLnZ4IC0uIGR4
ICouIGInLm1hc3MgKi4gbWFnOwogICAgICBiLnZ5IDwtIGIudnkgLS4gZHkgKi4gYicubWFzcyAq
LiBtYWc7CiAgICAgIGIudnogPC0gYi52eiAtLiBkeiAqLiBiJy5tYXNzICouIG1hZzsKCiAgICAg
IGInLnZ4IDwtIGInLnZ4ICsuIGR4ICouIGIubWFzcyAqLiBtYWc7CiAgICAgIGInLnZ5IDwtIGIn
LnZ5ICsuIGR5ICouIGIubWFzcyAqLiBtYWc7CiAgICAgIGInLnZ6IDwtIGInLnZ6ICsuIGR6ICou
IGIubWFzcyAqLiBtYWc7CiAgICBkb25lCiAgZG9uZTsKICBmb3IgaSA9IDAgdG8gbiBkbwogICAg
bGV0IGIgPSBib2RpZXMuKGkpIGluCiAgICBiLnggPC0gYi54ICsuIGR0ICouIGIudng7CiAgICBi
LnkgPC0gYi55ICsuIGR0ICouIGIudnk7CiAgICBiLnogPC0gYi56ICsuIGR0ICouIGIudno7CiAg
ZG9uZQoKCmxldCBlbmVyZ3kgYm9kaWVzID0KICBsZXQgZSA9IHJlZiAwLiBpbgogIGZvciBpID0g
MCB0byBBcnJheS5sZW5ndGggYm9kaWVzIC0gMSBkbwogICAgbGV0IGIgPSBib2RpZXMuKGkpIGlu
CiAgICBlIDo9ICFlICsuIDAuNSAqLiBiLm1hc3MgKi4gKGIudnggKi4gYi52eCArLiBiLnZ5ICou
IGIudnkgKy4gYi52eiAqLiBiLnZ6KTsKICAgIGZvciBqID0gaSsxIHRvIEFycmF5Lmxlbmd0aCBi
b2RpZXMgLSAxIGRvCiAgICAgIGxldCBiJyA9IGJvZGllcy4oaikgaW4KICAgICAgbGV0IGR4ID0g
Yi54IC0uIGInLnggIGFuZCBkeSA9IGIueSAtLiBiJy55ICBhbmQgZHogPSBiLnogLS4gYicueiBp
bgogICAgICBsZXQgZGlzdGFuY2UgPSBzcXJ0KGR4ICouIGR4ICsuIGR5ICouIGR5ICsuIGR6ICou
IGR6KSBpbgogICAgICBlIDo9ICFlIC0uIChiLm1hc3MgKi4gYicubWFzcykgLy4gZGlzdGFuY2UK
ICAgIGRvbmUKICBkb25lOwogICFlCgoKbGV0IG9mZnNldF9tb21lbnR1bSBib2RpZXMgPQogIGxl
dCBweCA9IHJlZiAwLiBhbmQgcHkgPSByZWYgMC4gYW5kIHB6ID0gcmVmIDAuIGluCiAgZm9yIGkg
PSAwIHRvIEFycmF5Lmxlbmd0aCBib2RpZXMgLSAxIGRvCiAgICBweCA6PSAhcHggKy4gYm9kaWVz
LihpKS52eCAqLiBib2RpZXMuKGkpLm1hc3M7CiAgICBweSA6PSAhcHkgKy4gYm9kaWVzLihpKS52
eSAqLiBib2RpZXMuKGkpLm1hc3M7CiAgICBweiA6PSAhcHogKy4gYm9kaWVzLihpKS52eiAqLiBi
b2RpZXMuKGkpLm1hc3M7CiAgZG9uZTsKICBib2RpZXMuKDApLnZ4IDwtIC0uICFweCAvLiBzb2xh
cl9tYXNzOwogIGJvZGllcy4oMCkudnkgPC0gLS4gIXB5IC8uIHNvbGFyX21hc3M7CiAgYm9kaWVz
LigwKS52eiA8LSAtLiAhcHogLy4gc29sYXJfbWFzcwoKCmxldCBqdXBpdGVyID0gewogIHggPSA0
Ljg0MTQzMTQ0MjQ2NDcyMDkwZSswMDsKICB5ID0gLTEuMTYwMzIwMDQ0MDI3NDI4MzllKzAwOwog
IHogPSAtMS4wMzYyMjA0NDQ3MTEyMzEwOWUtMDE7CiAgdnggPSAxLjY2MDA3NjY0Mjc0NDAzNjk0
ZS0wMyAqLiBkYXlzX3Blcl95ZWFyOwogIHZ5ID0gNy42OTkwMTExODQxOTc0MDQyNWUtMDMgKi4g
ZGF5c19wZXJfeWVhcjsKICB2eiA9IC02LjkwNDYwMDE2OTcyMDYzMDIzZS0wNSAqLiBkYXlzX3Bl
cl95ZWFyOwogIG1hc3MgPSA5LjU0NzkxOTM4NDI0MzI2NjA5ZS0wNCAqLiBzb2xhcl9tYXNzOwp9
CgpsZXQgc2F0dXJuID0gewogIHggPSA4LjM0MzM2NjcxODI0NDU3OTg3ZSswMDsKICB5ID0gNC4x
MjQ3OTg1NjQxMjQzMDQ3OWUrMDA7CiAgeiA9IC00LjAzNTIzNDE3MTE0MzIxMzgxZS0wMTsKICB2
eCA9IC0yLjc2NzQyNTEwNzI2ODYyNDExZS0wMyAqLiBkYXlzX3Blcl95ZWFyOwogIHZ5ID0gNC45
OTg1MjgwMTIzNDkxNzIzOGUtMDMgKi4gZGF5c19wZXJfeWVhcjsKICB2eiA9IDIuMzA0MTcyOTc1
NzM3NjM5MjllLTA1ICouIGRheXNfcGVyX3llYXI7CiAgbWFzcyA9IDIuODU4ODU5ODA2NjYxMzA4
MTJlLTA0ICouIHNvbGFyX21hc3M7Cn0KCmxldCB1cmFudXMgPSB7CiAgeCA9IDEuMjg5NDM2OTU2
MjEzOTEzMTBlKzAxOwogIHkgPSAtMS41MTExMTUxNDAxNjk4NjMxMmUrMDE7CiAgeiA9IC0yLjIz
MzA3NTc4ODkyNjU1NzM0ZS0wMTsKICB2eCA9IDIuOTY0NjAxMzc1NjQ3NjE2MThlLTAzICouIGRh
eXNfcGVyX3llYXI7CiAgdnkgPSAyLjM3ODQ3MTczOTU5NDgwOTUwZS0wMyAqLiBkYXlzX3Blcl95
ZWFyOwogIHZ6ID0gLTIuOTY1ODk1Njg1NDAyMzc1NTZlLTA1ICouIGRheXNfcGVyX3llYXI7CiAg
bWFzcyA9IDQuMzY2MjQ0MDQzMzUxNTYyOThlLTA1ICouIHNvbGFyX21hc3M7Cn0KCmxldCBuZXB0
dW5lID0gewogIHggPSAxLjUzNzk2OTcxMTQ4NTA5MTY1ZSswMTsKICB5ID0gLTIuNTkxOTMxNDYw
OTk4Nzk2NDFlKzAxOwogIHogPSAxLjc5MjU4NzcyOTUwMzcxMTgxZS0wMTsKICB2eCA9IDIuNjgw
Njc3NzI0OTAzODkzMjJlLTAzICouIGRheXNfcGVyX3llYXI7CiAgdnkgPSAxLjYyODI0MTcwMDM4
MjQyMjk1ZS0wMyAqLiBkYXlzX3Blcl95ZWFyOwogIHZ6ID0gLTkuNTE1OTIyNTQ1MTk3MTU4NzBl
LTA1ICouIGRheXNfcGVyX3llYXI7CiAgbWFzcyA9IDUuMTUxMzg5MDIwNDY2MTE0NTFlLTA1ICou
IHNvbGFyX21hc3M7Cn0KCmxldCBzdW4gPSB7CiAgeCA9IDAuOyAgeSA9IDAuOyAgeiA9IDAuOyAg
dnggPSAwLjsgIHZ5ID0gMC47IHZ6ID0gMC47CiAgbWFzcz0gc29sYXJfbWFzczsKfQoKbGV0IGJv
ZGllcyA9IFt8IHN1bjsganVwaXRlcjsgc2F0dXJuOyB1cmFudXM7IG5lcHR1bmUgfF0KCmxldCAo
KSA9CiAgbGV0IG4gPSBpbnRfb2Zfc3RyaW5nKFN5cy5hcmd2LigxKSkgaW4KICBvZmZzZXRfbW9t
ZW50dW0gYm9kaWVzOwogIFByaW50Zi5wcmludGYgIiUuOWZcbiIgKGVuZXJneSBib2RpZXMpOwog
IGZvciBpID0gMSB0byBuIGRvCiAgICBhZHZhbmNlIGJvZGllcyAwLjAxCiAgZG9uZTsKICBQcmlu
dGYucHJpbnRmICIlLjlmXG4iIChlbmVyZ3kgYm9kaWVzKQo=
--=-=-=
Content-Type: application/octet-stream
Content-Disposition: attachment; filename=nbody.c

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define pi 3.141592653589793
#define solar_mass (4 * pi * pi)
#define days_per_year 365.24

struct planet {
  double x, y, z;
  double vx, vy, vz;
  double mass;
};

void advance(int nbodies, struct planet * bodies, double dt)
{
  int i, j;

  for (i = 0; i < nbodies; i++) {
    struct planet * b = &(bodies[i]);
    for (j = i + 1; j < nbodies; j++) {
      struct planet * b2 = &(bodies[j]);
      double dx = b->x - b2->x;
      double dy = b->y - b2->y;
      double dz = b->z - b2->z;
      double distance = sqrt(dx * dx + dy * dy + dz * dz);
      double mag = dt / (distance * distance * distance);
      b->vx -= dx * b2->mass * mag;
      b->vy -= dy * b2->mass * mag;
      b->vz -= dz * b2->mass * mag;
      b2->vx += dx * b->mass * mag;
      b2->vy += dy * b->mass * mag;
      b2->vz += dz * b->mass * mag;
    }
  }
  for (i = 0; i < nbodies; i++) {
    struct planet * b = &(bodies[i]);
    b->x += dt * b->vx;
    b->y += dt * b->vy;
    b->z += dt * b->vz;
  }
}

double energy(int nbodies, struct planet * bodies)
{
  double e;
  int i, j;

  e = 0.0;
  for (i = 0; i < nbodies; i++) {
    struct planet * b = &(bodies[i]);
    e += 0.5 * b->mass * (b->vx * b->vx + b->vy * b->vy + b->vz * b->vz);
    for (j = i + 1; j < nbodies; j++) {
      struct planet * b2 = &(bodies[j]);
      double dx = b->x - b2->x;
      double dy = b->y - b2->y;
      double dz = b->z - b2->z;
      double distance = sqrt(dx * dx + dy * dy + dz * dz);
      e -= (b->mass * b2->mass) / distance;
    }
  }
  return e;
}

void offset_momentum(int nbodies, struct planet * bodies)
{
  double px = 0.0, py = 0.0, pz = 0.0;
  int i;
  for (i = 0; i < nbodies; i++) {
    px += bodies[i].vx * bodies[i].mass;
    py += bodies[i].vy * bodies[i].mass;
    pz += bodies[i].vz * bodies[i].mass;
  }
  bodies[0].vx = - px / solar_mass;
  bodies[0].vy = - py / solar_mass;
  bodies[0].vz = - pz / solar_mass;
}

#define NBODIES 5
struct planet bodies[NBODIES] = {
  {                               /* sun */
    0, 0, 0, 0, 0, 0, solar_mass
  },
  {                               /* jupiter */
    4.84143144246472090e+00,
    -1.16032004402742839e+00,
    -1.03622044471123109e-01,
    1.66007664274403694e-03 * days_per_year,
    7.69901118419740425e-03 * days_per_year,
    -6.90460016972063023e-05 * days_per_year,
    9.54791938424326609e-04 * solar_mass
  },
  {                               /* saturn */
    8.34336671824457987e+00,
    4.12479856412430479e+00,
    -4.03523417114321381e-01,
    -2.76742510726862411e-03 * days_per_year,
    4.99852801234917238e-03 * days_per_year,
    2.30417297573763929e-05 * days_per_year,
    2.85885980666130812e-04 * solar_mass
  },
  {                               /* uranus */
    1.28943695621391310e+01,
    -1.51111514016986312e+01,
    -2.23307578892655734e-01,
    2.96460137564761618e-03 * days_per_year,
    2.37847173959480950e-03 * days_per_year,
    -2.96589568540237556e-05 * days_per_year,
    4.36624404335156298e-05 * solar_mass
  },
  {                               /* neptune */
    1.53796971148509165e+01,
    -2.59193146099879641e+01,
    1.79258772950371181e-01,
    2.68067772490389322e-03 * days_per_year,
    1.62824170038242295e-03 * days_per_year,
    -9.51592254519715870e-05 * days_per_year,
    5.15138902046611451e-05 * solar_mass
  }
};

int main(int argc, char ** argv)
{
  int n = atoi(argv[1]);
  int i;

  offset_momentum(NBODIES, bodies);
  printf ("%.9f\n", energy(NBODIES, bodies));
  for (i = 1; i <= n; i++)
    advance(NBODIES, bodies, 0.01);
  printf ("%.9f\n", energy(NBODIES, bodies));
  return 0;
}


--=-=-=


-- 
let () = let rec f a w i j = Printf.printf "%.20f\r" a; let a1 = a *. i /. j in
if w then f a1 false (i +. 2.0) j else f a1 true i (j +. 2.0) in f 2.0 false 2.0 1.0



--=-=-=--