Hastebin

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

;; racket-nbody.rkt
;; NBODY in Racket programming language. Columned version.
;; compile: $ raco exe racket-nbody.rkt
;;     run: $ ./racket-nbody <n-iterations>
;;
;; unsafe-flvector-* is hot-hot-hot
;;

(module energy racket
  (provide take-energy)
  (require (only-in racket/math pi)
           racket/flonum
           racket/fixnum
           racket/unsafe/ops)

  (define-syntax (define-constant stx)
    (syntax-case stx () [(_ name constvalue)
                         (syntax (define-syntax (name s)
                                   (syntax-case s () (_ (identifier? #'name)  (syntax constvalue)))))]))

  (define-constant SOLAR-MASS (* 4 pi pi))
  (define-constant DAYS-PER-YEAR 365.24)
  (define-syntax-rule (dpe v) (* v DAYS-PER-YEAR))
  (define-syntax-rule (slr m) (* m SOLAR-MASS))
  (define-syntax-rule (fx1+ v) (fx+ v 1))
  (define-syntax-rule (fx1- v) (fx- v 1))
  (define-syntax-rule (fxzero? v) (fx= v 0))
  (define-syntax-rule (^2 x) (fl* x x))
  (define-syntax sum^2 (syntax-rules () ((_ x y z) (fl+ (^2 x) (^2 y) (^2 z)))))

  (define-syntax-rule (dsqt x) (let ((v x))(fl* v (flsqrt v))))
  (define-syntax-rule (dotimes n body ...) (do ((i_ (fx+ n 0) (fx1- i_)))((fxzero? i_)) 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-constant body-size 7)
  (define bodies
    (flvector
     SOLAR-MASS 0.0 0.0 0.0 0.0 0.0 0.0
     (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 n-body (fxquotient (flvector-length bodies) body-size))
  (define-syntax-rule (mkcol n) (for/flvector ((i (in-range 0 n-body) ))
                                  (flvector-ref bodies (+ n (* i body-size)))))

  (define-syntax (let-thval stx)
    (syntax-case stx ()
      [(_ ((a ix col) ...) body ...)
       #'(let-syntax ((a (make-set!-transformer
                          (lambda (stx)
                            (syntax-case stx (set!)
                              [a (identifier? (syntax a)) (syntax (unsafe-flvector-ref col ix))]
                              [(set! a val) (syntax (unsafe-flvector-set! col ix val))])))) ...)
           body ...)]))

  (define-syntax let-getset
    (syntax-rules ()
      [(_ ((a ix col) ...) body ...)
       (let-syntax ((a (syntax-rules ()
                         ((a op delta)
                          (unsafe-flvector-set!
                           col ix (op (unsafe-flvector-ref col ix) delta))))) ...)
         body ...)]))


  (define (take-energy n)
    (define-values (mass-vec x-vec y-vec z-vec vx-vec vy-vec vz-vec)
      (values (mkcol 0) (mkcol 1) (mkcol 2) (mkcol 3) (mkcol 4) (mkcol 5) (mkcol 6)))

    (define-constant dt 0.01)
    (define E 0.0)

    (dotimes
     n
     (do ((i 0 (fx1+ i)))
         ((fx= i n-body))
       (do ((j (fx1+ i) (fx1+ j)))
           ((fx= j n-body))
         (let-thval
          ((mass_i i mass-vec)
           (mass_j j mass-vec)
           (Xi i x-vec) (Yi i y-vec) (Zi i z-vec) (Xj j x-vec) (Yj j y-vec) (Zj j z-vec))
          (let-getset
           ((VXi i vx-vec)
            (VYi i vy-vec)
            (VZi i vz-vec)
            (VXj j vx-vec)
            (VYj j vy-vec)
            (VZj j vz-vec))
           (let* ((dx (fl- Xi Xj))
                  (dy (fl- Yi Yj))
                  (dz (fl- Zi Zj))
                  (mag (fl/ dt (dsqt (sum^2 dx dy dz)))))
             (VXi fl- (fl* dx mass_j mag))
             (VYi fl- (fl* dy mass_j mag))
             (VZi fl- (fl* dz mass_j mag))
             (VXj fl+ (fl* dx mass_i mag))
             (VYj fl+ (fl* dy mass_i mag))
             (VZj fl+ (fl* dz mass_i mag))
             )))))
     (do ((i 0 (fx1+ i)))
         ((fx= i n-body))
       (let-thval
        ((VX i vx-vec) (VY i vy-vec) (VZ i vz-vec))
        (let-getset
         ((X i x-vec)(Y i y-vec) (Z i z-vec))
         (X fl+ (fl* dt VX))
         (Y fl+ (fl* dt VY))
         (Z fl+ (fl* dt VZ))
         )))
     )

    (do ((i 0 (fx1+ i)))
        ((fx= i n-body))
      (let-thval
       ((mass_i i mass-vec)
        (VX i vx-vec)
        (VY i vy-vec)
        (VZ i vz-vec))
       (addn! E (fl* 0.5 mass_i (sum^2 VX VY VZ)))
      (do ((j (fx1+ i) (fx1+ j)))
          ((fx= j n-body))
        (let-thval
         ((Xi i x-vec) (Yi i y-vec) (Zi i z-vec)
          (Xj j x-vec) (Yj j y-vec) (Zj j z-vec) (mass_j j mass-vec))
         (subn! E (fl/ (fl* mass_i mass_j)
                       (flsqrt (sum^2 (fl- Xi Xj)
                                      (fl- Yi Yj)
                                      (fl- Zi Zj)))))))))
    E
    )
) ;; end-of-module energy

(require 'energy)
(provide (all-from-out 'energy))

;; (require (only-in 'energy take-energy))

(module* main #f
  (define-syntax-rule (->sec from to) (* 1e-3 (- from to)))
  (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 "racket-nbody"
   #:args args
   (cond ((pair? args) (test (string->number (car args))))
         #;(else (take-energy 100000))
         (else (printf "Usage: ~a <n-iterations>~%" "nbody"))))
)