Hastebin

;;; N-body benchmark - ISLisp version
;;; nbody-islisp.lsp
;;; compile eisl > (compile-file "nbody-islisp.lsp")
;;; ...
;;; T
;;; eisl > (load "nbody-islisp.o")
;;; eisl > (time (take-energy 50000000))
;;;  it takes a Time to calculate... about 2600 seconds.

(defconstant +solar-mass+ (* 4.0 3.141592653589793 3.141592653589793))
(defconstant +days-per-year+ 365.24)
(defconstant +dt+ 0.01)

(defglobal *x*
           (vector 0.0
                   4.84143144246472090e+00
                   8.34336671824457987e+00
                   1.28943695621391310e+01
                   1.53796971148509165e+01))

(defglobal *y*
           (vector 0.0
                   -1.16032004402742839e+00
                   4.12479856412430479e+00
                   -1.51111514016986312e+01
                   -2.59193146099879641e+01))

(defglobal *z*
           (vector 0.0
                   -1.03622044471123109e-01
                   -4.03523417114321381e-01
                   -2.23307578892655734e-01
                   1.79258772950371181e-01))

(defglobal *vx*
           (vector 0.0
                   (* 1.66007664274403694e-03 +days-per-year+)
                   (* -2.76742510726862411e-03 +days-per-year+)
                   (* 2.96460137564761618e-03 +days-per-year+)
                   (* 2.68067772490389322e-03 +days-per-year+)))

(defglobal *vy*
           (vector 0.0
                   (* 7.69901118419740425e-03 +days-per-year+)
                   (* 4.99852801234917238e-03 +days-per-year+)
                   (* 2.37847173959480950e-03 +days-per-year+)
                   (* 1.62824170038242295e-03 +days-per-year+)))

(defglobal *vz*
           (vector 0.0
                   (* -6.90460016972063023e-05 +days-per-year+)
                   (* 2.30417297573763929e-05 +days-per-year+)
                   (* -2.96589568540237556e-05 +days-per-year+)
                   (* -9.51592254519715870e-05 +days-per-year+)))

(defglobal *mass*
           (vector +solar-mass+
                   (* 9.54791938424326609e-04 +solar-mass+)
                   (* 2.85885980666130812e-04 +solar-mass+)
                   (* 4.36624404335156298e-05 +solar-mass+)
                   (* 5.15138902046611451e-05 +solar-mass+)))

(defmacro incf (var addval)
  `(setf ,var (+ ,var ,addval)))

(defmacro decf (var decval)
  `(setf ,var (- ,var ,decval)))

(defun advance ()
  (let ((x *x*) (y *y*) (z *z*)
        (vx *vx*) (vy *vy*) (vz *vz*)
        (mass *mass*))
    (for ((i 0 (+ i 1)))
         ((= i 5))
      (let ((xi (aref x i)) (yi (aref y i)) (zi (aref z i))
            (vxi (aref vx i)) (vyi (aref vy i)) (vzi (aref vz i))
            (mi (aref mass i)))
        (for ((j (+ i 1) (+ j 1)))
              ((= j 5))
          (let* ((dx (- xi (aref x j)))
                 (dy (- yi (aref y j)))
                 (dz (- zi (aref z j)))
                 (dsq (+ (* dx dx) (* dy dy) (* dz dz)))
                 (mag (quotient +dt+ (* dsq (sqrt dsq))))
                 (mj-mag (* (aref mass j) mag))
                 (mi-mag (* mi mag)))
            (decf vxi (* dx mj-mag))
            (decf vyi (* dy mj-mag))
            (decf vzi (* dz mj-mag))
            (incf (aref vx j) (* dx mi-mag))
            (incf (aref vy j) (* dy mi-mag))
            (incf (aref vz j) (* dz mi-mag))))
        (setf (aref vx i) vxi)
        (setf (aref vy i) vyi)
        (setf (aref vz i) vzi)))
    (for ((i 0 (+ i 1)))
         ((= i 5))
         (incf (aref x i) (* +dt+ (aref vx i)))
         (incf (aref y i) (* +dt+ (aref vy i)))
         (incf (aref z i) (* +dt+ (aref vz i))))
  ))

(defun energy ()
  (let ((x *x*) (y *y*) (z *z*)
        (vx *vx*) (vy *vy*) (vz *vz*)
        (mass *mass*)
        (e 0.0))
    (for ((i 0 (+ i 1)))
         ((= i 5))
      (incf e (* 0.5 (aref mass i)
                (+ (* (aref vx i) (aref vx i))
                   (* (aref vy i) (aref vy i))
                   (* (aref vz i) (aref vz i)))))
      (for ((j (+ i 1) (+ j 1)))
           ((= j 5))
        (let* ((dx (- (aref x i) (aref x j)))
               (dy (- (aref y i) (aref y j)))
               (dz (- (aref z i) (aref z j))))
          (decf e (quotient (* (aref mass i) (aref mass j))
                       (sqrt (+ (* dx dx) (* dy dy) (* dz dz))))))))
    e))


(defun take-energy (n)
  ;;(declare (type fixnum n))
  (for ((i n (- i 1)))
       ((= i 0))
       (advance))
  (format (standard-output) "~A~%" (energy)))