Hastebin

#!/usr/bin/env racket
#lang racket

(require racket/fixnum
         racket/flonum)

(provide take-energy)

(define PI pi)
(define SOLAR-MASS (* 4 PI PI))
(define DAYS-PER-YEAR 365.24)
(define-syntax dpe (syntax-rules () ((_ v) (* v DAYS-PER-YEAR))))
(define-syntax slr (syntax-rules () ((_ m) (* m SOLAR-MASS))))
(define-syntax ^2 (syntax-rules () ((_ x) (fl* x x))))
(define-syntax dotimes (syntax-rules () ((_ n body ...) (do ((i 0 (fx+ 1 i)))((fx= i n)) body ...))))
(define-syntax subn! (syntax-rules () ((_ what withs) (set! what (fl- what withs)))))
(define-syntax addn! (syntax-rules () ((_ what withs) (set! what (fl+ what withs)))))


(define-syntax _mass
  (syntax-rules ()
    ((_ v base) (flvector-ref v base))))

(define-syntax _x
  (syntax-rules ()
    ((_ v base) (flvector-ref v (fx+ 1 base)))))

(define-syntax _y
  (syntax-rules ()
    ((_ v base) (flvector-ref v (fx+ 2 base)))))

(define-syntax _z
  (syntax-rules ()
    ((_ v base) (flvector-ref v (fx+ 3 base)))))

(define-syntax _vx
  (syntax-rules ()
    ((_ v base) (flvector-ref v (fx+ 4 base)))))

(define-syntax _vy
  (syntax-rules ()
    ((_ v base) (flvector-ref v (fx+ 5 base)))))

(define-syntax _vz
  (syntax-rules ()
    ((_ v base) (flvector-ref v (fx+ 6 base)))))

;; --- setters ----

(define-syntax _x!
  (syntax-rules ()
    ((_ v base val) (flvector-set! v (fx+ 1 base) val))))

(define-syntax _y!
  (syntax-rules ()
    ((_ v base val) (flvector-set! v (fx+ 2 base) val))))

(define-syntax _z!
  (syntax-rules ()
    ((_ v base val) (flvector-set! v (fx+ 3 base) val))))

(define-syntax _vx!
  (syntax-rules ()
    ((_ v base val) (flvector-set! v (fx+ 4 base) val))))

(define-syntax _vy!
  (syntax-rules ()
    ((_ v base val) (flvector-set! v (fx+ 5 base) val))))

(define-syntax _vz!
  (syntax-rules ()
    ((_ v base val) (flvector-set! v (fx+ 6 base) val))))


(define-syntax sum^2
  (syntax-rules ()
    ((_ x y z) (fl+ (^2 x) (^2 y) (^2 z)))))

;;(: take-energy (Integer -> Flonum))
(define (take-energy n)
;;  (: bodies (Vectorof Flonum))
  (define bodies
    (flvector
     SOLAR-MASS 0.0 0.0 0.0 0.0 0.0 0.0  ;; SUN
     (slr 9.54791938424326609e-04) 4.84143144246472090e+00 -1.16032004402742839e+00 -1.03622044471123109e-01 (dpe 1.66007664274403694e-03) (dpe 7.69901118419740425e-03) (dpe -6.90460016972063023e-05)
     (slr 2.85885980666130812e-04) 8.34336671824457987e+00 4.12479856412430479e+00 -4.03523417114321381e-01 (dpe -2.76742510726862411e-03) (dpe 4.99852801234917238e-03) (dpe 2.30417297573763929e-05)
     (slr 4.36624404335156298e-05) 1.28943695621391310e+01 -1.51111514016986312e+01 -2.23307578892655734e-01 (dpe 2.96460137564761618e-03) (dpe 2.37847173959480950e-03) (dpe -2.96589568540237556e-05)
     (slr 5.15138902046611451e-05) 1.53796971148509165e+01 -2.59193146099879641e+01 1.79258772950371181e-01 (dpe 2.68067772490389322e-03) (dpe 1.62824170038242295e-03) (dpe -9.51592254519715870e-05)
     ))

  (define dt 0.01)
  (define bodies-len (flvector-length bodies))
  (define body-size 7)
  (dotimes
      n
      (for ((i (in-range 0 bodies-len body-size)))
        (for ((j (in-range (+ i body-size) bodies-len body-size)))
          (let* ((dx (fl- (_x bodies i) (_x bodies j)))
                 (dy (fl- (_y bodies i) (_y bodies j)))
                 (dz (fl- (_z bodies i) (_z bodies j)))
                 (dsq (sum^2 dx dy dz))
                 (mag (fl/ dt (fl* dsq (sqrt dsq)))))
            (_vx! bodies i (fl- (_vx bodies i) (fl* dx (_mass bodies j) mag)))
            (_vy! bodies i (fl- (_vy bodies i) (fl* dy (_mass bodies j) mag)))
            (_vz! bodies i (fl- (_vz bodies i) (fl* dz (_mass bodies j) mag)))
            (_vx! bodies j (fl+ (_vx bodies j) (fl* dx (_mass bodies i) mag)))
            (_vy! bodies j (fl+ (_vy bodies j) (fl* dy (_mass bodies i) mag)))
            (_vz! bodies j (fl+ (_vz bodies j) (fl* dz (_mass bodies i) mag))))
          ))
    (for ((i (in-range 0 bodies-len body-size)))
      (_x! bodies i (fl+ (_x bodies i) (fl* dt (_vx bodies i))))
      (_y! bodies i (fl+ (_y bodies i) (fl* dt (_vy bodies i))))
      (_z! bodies i (fl+ (_z bodies i) (fl* dt (_vz bodies i))))))
  (do ((i 0 (+ body-size i))
       (E 0.0 (fl+ E (* 0.5 (_mass bodies i) (sum^2 (_vx bodies i) (_vy bodies i) (_vz bodies i))))))
      ((fx= i bodies-len) E)
    (do ((j (+ i body-size) (+ body-size j)))
        ((= j bodies-len))
      (let ((dx (- (_x bodies i) (_x bodies j)))
            (dy (- (_y bodies i) (_y bodies j)))
            (dz (- (_z bodies i) (_z bodies j))))
        (subn! E (/ (* (_mass bodies i) (_mass bodies j))
                    (sqrt (sum^2 dx dy dz)))))))
  )

(define-syntax-rule (->sec from to) (* 1e-3 (- from to)))

(module* main #f
  (define (test n)
    (let* ((st (current-inexact-milliseconds))
           (e (take-energy n))
           (et (current-inexact-milliseconds)))
      (printf "~a  per ~a sec ~%" (~r e #:precision 9)
              (~r (->sec et st) #:precision 3))))
  (command-line
   #:program "nbody"
   #:args args
   (cond ((pair? args) (test (string->number (car args))))
         (else (take-energy 100000))
         #;(else (printf "Usage: ~a <n-iterations>~%" "nbody"))))
)