source code
;; The Computer Language Benchmarks Game
;; https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
;;;
;;; contributed by Patrick Frankenberger
;;; modified by Juho Snellman 2005-11-18
;;; * About 40% speedup on SBCL, 90% speedup on CMUCL
;;; * Represent a body as a DEFSTRUCT with (:TYPE VECTOR DOUBLE-FLOAT), a
;;; not as a structure that contains vectors
;;; * Inline APPLYFORCES
;;; * Replace (/ DT DISTANCE DISTANCE DISTANCE) with
;;; (/ DT (* DISTANCE DISTANCE DISTANCE)), as is done in the other
;;; implementations of this test.
;;; * Add a couple of declarations
;;; * Heavily rewritten for style (represent system as a list instead of
;;; an array to make the nested iterations over it less clumsy, use
;;; INCF/DECF where appropriate, break very long lines, etc)
;;; modified by Marko Kocic
;;; * add optimization declarations
(declaim (optimize (speed 3)(safety 0)(space 0)(debug 0)))
(defconstant +days-per-year+ 365.24d0)
(defconstant +solar-mass+ (* 4d0 pi pi))
(defstruct (body (:type (vector double-float))
(:conc-name nil)
(:constructor make-body (x y z vx vy vz mass)))
x y z
vx vy vz
mass)
(deftype body () '(vector double-float 7))
(defparameter *jupiter*
(make-body 4.84143144246472090d0
-1.16032004402742839d0
-1.03622044471123109d-1
(* 1.66007664274403694d-3 +days-per-year+)
(* 7.69901118419740425d-3 +days-per-year+)
(* -6.90460016972063023d-5 +days-per-year+)
(* 9.54791938424326609d-4 +solar-mass+)))
(defparameter *saturn*
(make-body 8.34336671824457987d0
4.12479856412430479d0
-4.03523417114321381d-1
(* -2.76742510726862411d-3 +days-per-year+)
(* 4.99852801234917238d-3 +days-per-year+)
(* 2.30417297573763929d-5 +days-per-year+)
(* 2.85885980666130812d-4 +solar-mass+)))
(defparameter *uranus*
(make-body 1.28943695621391310d1
-1.51111514016986312d1
-2.23307578892655734d-1
(* 2.96460137564761618d-03 +days-per-year+)
(* 2.37847173959480950d-03 +days-per-year+)
(* -2.96589568540237556d-05 +days-per-year+)
(* 4.36624404335156298d-05 +solar-mass+)))
(defparameter *neptune*
(make-body 1.53796971148509165d+01
-2.59193146099879641d+01
1.79258772950371181d-01
(* 2.68067772490389322d-03 +days-per-year+)
(* 1.62824170038242295d-03 +days-per-year+)
(* -9.51592254519715870d-05 +days-per-year+)
(* 5.15138902046611451d-05 +solar-mass+)))
(defparameter *sun*
(make-body 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 +solar-mass+))
(declaim (inline applyforces))
(defun applyforces (a b dt)
(declare (type body a b) (type double-float dt))
(let* ((dx (- (x a) (x b)))
(dy (- (y a) (y b)))
(dz (- (z a) (z b)))
(distance (sqrt (+ (* dx dx) (* dy dy) (* dz dz))))
(mag (/ dt (* distance distance distance)))
(dxmag (* dx mag))
(dymag (* dy mag))
(dzmag (* dz mag)))
(decf (vx a) (* dxmag (mass b)))
(decf (vy a) (* dymag (mass b)))
(decf (vz a) (* dzmag (mass b)))
(incf (vx b) (* dxmag (mass a)))
(incf (vy b) (* dymag (mass a)))
(incf (vz b) (* dzmag (mass a))))
nil)
(defun advance (system dt)
(declare (double-float dt))
(loop for (a . rest) on system do
(dolist (b rest)
(applyforces a b dt)))
(dolist (b system)
(incf (x b) (* dt (vx b)))
(incf (y b) (* dt (vy b)))
(incf (z b) (* dt (vz b))))
nil)
(defun energy (system)
(let ((e 0.0d0))
(declare (double-float e))
(loop for (a . rest) on system do
(incf e (* 0.5d0
(mass a)
(+ (* (vx a) (vx a))
(* (vy a) (vy a))
(* (vz a) (vz a)))))
(dolist (b rest)
(let* ((dx (- (x a) (x b)))
(dy (- (y a) (y b)))
(dz (- (z a) (z b)))
(dist (sqrt (+ (* dx dx) (* dy dy) (* dz dz)))))
(decf e (/ (* (mass a) (mass b)) dist)))))
e))
(defun offset-momentum (system)
(let ((px 0.0d0)
(py 0.0d0)
(pz 0.0d0))
(dolist (p system)
(incf px (* (vx p) (mass p)))
(incf py (* (vy p) (mass p)))
(incf pz (* (vz p) (mass p))))
(setf (vx (car system)) (/ (- px) +solar-mass+)
(vy (car system)) (/ (- py) +solar-mass+)
(vz (car system)) (/ (- pz) +solar-mass+))
nil))
(defun nbody (n)
(declare (fixnum n))
(let ((system (list *sun* *jupiter* *saturn* *uranus* *neptune*)))
(offset-momentum system)
(format t "~,9F~%" (energy system))
(dotimes (i n)
(advance system 0.01d0))
(format t "~,9F~%" (energy system))))
(defun main ()
(let ((n (parse-integer (or (car (last #+sbcl sb-ext:*posix-argv*
#+cmu extensions:*command-line-strings*
#+gcl si::*command-args*)) "1"))))
(nbody n)))
notes, command-line, and program output
NOTES:
64-bit Ubuntu quad core
SBCL 2.4.2
Mon, 04 Mar 2024 23:22:34 GMT
MAKE:
cp: 'nbody.sbcl-2.sbcl' and './nbody.sbcl-2.sbcl' are the same file
SBCL built with: /opt/src/sbcl-2.4.2/bin/sbcl --userinit /dev/null --batch --eval '(load "nbody.sbcl-2.sbcl_compile")'
### START nbody.sbcl-2.sbcl_compile
(handler-bind ((sb-ext:defconstant-uneql (lambda (c) (abort c)))) (require :sb-concurrency) (load (compile-file "nbody.sbcl-2.sbcl" ))) (save-lisp-and-die "sbcl.core" :purify t)
### END nbody.sbcl-2.sbcl_compile
; compiling file "/home/dunham/all-benchmarksgame/benchmarksgame_i53330/nbody/tmp/nbody.sbcl-2.sbcl" (written 26 APR 2018 12:49:01 PM):
; file: /home/dunham/all-benchmarksgame/benchmarksgame_i53330/nbody/tmp/nbody.sbcl-2.sbcl
; in: DEFSTRUCT BODY
; (DEFSTRUCT
; (BODY (:TYPE (VECTOR DOUBLE-FLOAT)) (:CONC-NAME NIL)
; (:CONSTRUCTOR MAKE-BODY (X Y Z VX VY VZ MASS)))
; X
; Y
; Z
; VX
; VY
; VZ
; MASS)
; --> DEFUN PROGN SB-IMPL::%DEFUN SB-IMPL::%DEFUN SB-INT:NAMED-LAMBDA
; ==>
; #'(SB-INT:NAMED-LAMBDA X
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK X (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 0))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; ==>
; #'(SB-INT:NAMED-LAMBDA Y
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK Y (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 1))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; ==>
; #'(SB-INT:NAMED-LAMBDA Z
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK Z (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 2))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; ==>
; #'(SB-INT:NAMED-LAMBDA VX
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK VX (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 3))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; ==>
; #'(SB-INT:NAMED-LAMBDA VY
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK VY (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 4))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; ==>
; #'(SB-INT:NAMED-LAMBDA VZ
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK VZ (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 5))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; ==>
; #'(SB-INT:NAMED-LAMBDA MASS
; (STRUCTURE)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (DECLARE (TYPE (SIMPLE-ARRAY DOUBLE-FLOAT (*)) STRUCTURE))
; (BLOCK MASS (THE (AND T DOUBLE-FLOAT) (ELT STRUCTURE 6))))
;
; note: doing float to pointer coercion (cost 13) to "<return value>"
; in: DEFUN ENERGY
; (DEFUN ENERGY (SYSTEM)
; (LET ((E 0.0d0))
; (DECLARE (DOUBLE-FLOAT E))
; (LOOP FOR (A . REST) ON SYSTEM
; DO (INCF E (* 0.5d0 # #)) (DOLIST (B REST)
; (LET* #
; #)))
; E))
; --> SB-IMPL::%DEFUN SB-IMPL::%DEFUN SB-INT:NAMED-LAMBDA
; ==>
; #'(SB-INT:NAMED-LAMBDA ENERGY
; (SYSTEM)
; (DECLARE (SB-C::TOP-LEVEL-FORM))
; (BLOCK ENERGY
; (LET ((E 0.0d0))
; (DECLARE (DOUBLE-FLOAT E))
; (LOOP FOR (A . REST) ON SYSTEM
; DO (INCF E #) (DOLIST # #))
; E)))
;
; note: doing float to pointer coercion (cost 13) from E to "<return value>"
; in: DEFUN OFFSET-MOMENTUM
; (/ (- PX) +SOLAR-MASS+)
;
; note: unable to
; convert to multiplication by reciprocal
; due to type uncertainty:
; The first argument is a (OR DOUBLE-FLOAT
; (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; (/ (- PY) +SOLAR-MASS+)
;
; note: unable to
; convert to multiplication by reciprocal
; due to type uncertainty:
; The first argument is a (OR DOUBLE-FLOAT
; (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; (/ (- PZ) +SOLAR-MASS+)
;
; note: unable to
; convert to multiplication by reciprocal
; due to type uncertainty:
; The first argument is a (OR DOUBLE-FLOAT
; (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; (INCF PX (* (VX P) (MASS P)))
; --> THE
; ==>
; (+ (* (VX P) (MASS P)) PX)
;
; note: forced to do GENERIC-+ (cost 10)
; unable to do inline float arithmetic (cost 2) because:
; The second argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES DOUBLE-FLOAT &OPTIONAL).
; unable to do inline float arithmetic (cost 4) because:
; The second argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
; DOUBLE-FLOAT).
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES (COMPLEX DOUBLE-FLOAT)
; &OPTIONAL).
; (INCF PY (* (VY P) (MASS P)))
; --> THE
; ==>
; (+ (* (VY P) (MASS P)) PY)
;
; note: forced to do GENERIC-+ (cost 10)
; unable to do inline float arithmetic (cost 2) because:
; The second argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES DOUBLE-FLOAT &OPTIONAL).
; unable to do inline float arithmetic (cost 4) because:
; The second argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
; DOUBLE-FLOAT).
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES (COMPLEX DOUBLE-FLOAT)
; &OPTIONAL).
; (INCF PZ (* (VZ P) (MASS P)))
; --> THE
; ==>
; (+ (* (VZ P) (MASS P)) PZ)
;
; note: forced to do GENERIC-+ (cost 10)
; unable to do inline float arithmetic (cost 2) because:
; The second argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES DOUBLE-FLOAT &OPTIONAL).
; unable to do inline float arithmetic (cost 4) because:
; The second argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
; DOUBLE-FLOAT).
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES (COMPLEX DOUBLE-FLOAT)
; &OPTIONAL).
; (- PX)
;
; note: forced to do GENERIC-NEGATE (cost 10)
; unable to do inline float arithmetic (cost 1) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES DOUBLE-FLOAT &OPTIONAL).
; unable to do inline float arithmetic (cost 1) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
; DOUBLE-FLOAT).
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES (COMPLEX DOUBLE-FLOAT)
; &OPTIONAL).
; (/ (- PX) +SOLAR-MASS+)
;
; note: forced to do full call
; unable to do inline float arithmetic (cost 19) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; (- PY)
;
; note: forced to do GENERIC-NEGATE (cost 10)
; unable to do inline float arithmetic (cost 1) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES DOUBLE-FLOAT &OPTIONAL).
; unable to do inline float arithmetic (cost 1) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
; DOUBLE-FLOAT).
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES (COMPLEX DOUBLE-FLOAT)
; &OPTIONAL).
; (/ (- PY) +SOLAR-MASS+)
;
; note: forced to do full call
; unable to do inline float arithmetic (cost 19) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; (- PZ)
;
; note: forced to do GENERIC-NEGATE (cost 10)
; unable to do inline float arithmetic (cost 1) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES DOUBLE-FLOAT &OPTIONAL).
; unable to do inline float arithmetic (cost 1) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
; DOUBLE-FLOAT).
; The result is a (VALUES (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT))
; &OPTIONAL), not a (VALUES (COMPLEX DOUBLE-FLOAT)
; &OPTIONAL).
; (/ (- PZ) +SOLAR-MASS+)
;
; note: forced to do full call
; unable to do inline float arithmetic (cost 19) because:
; The first argument is a (OR DOUBLE-FLOAT (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
; (INCF PX (* (VX P) (MASS P)))
; --> THE
; ==>
; (+ (* (VX P) (MASS P)) PX)
;
; note: doing float to pointer coercion (cost 13), for:
; the first argument of GENERIC-+
; (INCF PY (* (VY P) (MASS P)))
; --> THE
; ==>
; (+ (* (VY P) (MASS P)) PY)
;
; note: doing float to pointer coercion (cost 13), for:
; the first argument of GENERIC-+
; (INCF PZ (* (VZ P) (MASS P)))
; --> THE
; ==>
; (+ (* (VZ P) (MASS P)) PZ)
;
; note: doing float to pointer coercion (cost 13), for:
; the first argument of GENERIC-+
;
; compilation unit finished
; printed 23 notes
; wrote /home/dunham/all-benchmarksgame/benchmarksgame_i53330/nbody/tmp/nbody.sbcl-2.fasl
; compilation finished in 0:00:00.172
### START nbody.sbcl-2.sbcl_run
(main) (quit)
### END nbody.sbcl-2.sbcl_run
3.49s to complete and log all make actions
COMMAND LINE:
/opt/src/sbcl-2.4.2/bin/sbcl --dynamic-space-size 500 --noinform --core sbcl.core --userinit /dev/null --load nbody.sbcl-2.sbcl_run 50000000
PROGRAM OUTPUT:
-0.169075164
-0.169059907