Hastebin

#!/usr/bin/env gsi-script
;; gsc-colnbody.scm
;; NBODY exercise. Gambit Scheme variant with body's record fields data rearranged into columns
;; compile: $ gsc -exe gsc-colnbody
;;     run: $ ./gsc-colnbody <n-iterations>
;;

(declare (standard-bindings)
         (not safe)
         ;;(debug)
         )

(define (take-energy n)
  (define PI (* 4.0 (atan 1)))
  (define SOLAR-MASS (* 4 PI PI))
  (define DAYS-PER-YEAR 365.24)

  (define-syntax dpe (syntax-rules () ((_ v) (fl* v DAYS-PER-YEAR))))
  (define-syntax slr (syntax-rules () ((_ m) (fl* m SOLAR-MASS))))

  (define bodies
    (f64vector
     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 body-size 7)
  (define bodies-len (f64vector-length bodies))    ;; (2)
  (define n-body (fxquotient bodies-len body-size))

  (define dt 0.01)
  (define-syntax mkcol (syntax-rules () ((_ ix) (do ((i ix (+ i body-size))
                                                     (j 0 (+ j 1))
                                                     (column (make-f64vector n-body)))
                                                    ((<= bodies-len i) column)
                                                  (f64vector-set! column j (f64vector-ref bodies i))))))
  (define mass-column (mkcol 0))
  (define x-column (mkcol 1))
  (define y-column (mkcol 2))
  (define z-column (mkcol 3))
  (define veloce-x (mkcol 4))
  (define veloce-y (mkcol 5))
  (define veloce-z (mkcol 6))

  (define-syntax let-thval
    (syntax-rules () ((_ ((a ix column) ...) body ...)
                      (let ((a (f64vector-ref column ix)) ...)
                        body ...))))

  (define-syntax let-getset
    (syntax-rules () ((_ ((a ix column) ...) body ...)
                      (begin
                        (define-syntax a
                          (syntax-rules ()
                            ((_ op arg) (f64vector-set! column ix
                                                        (op
                                                         (f64vector-ref column ix)
                                                         arg))))) ...
                                                         body ...))))

  (define-syntax dotimes (syntax-rules () ((_ n body ...) (do ((i_ n (fx- i_ 1))) ((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-syntax ^2 (syntax-rules () ((_ x) (flsquare x))))
  (define-syntax sum^2 (syntax-rules () ((_ x y z) (fl+ (^2 x) (^2 y) (^2 z)))))
  (define-syntax dsqt (syntax-rules () ((_ x) (let ((v x)) (fl* v (flsqrt v))))))
  (define E 0.0)

  (dotimes
   n
   (do ((i 0 (fx+ i 1)))
       ((fx= i n-body))
     (let-thval
      ((massi i mass-column) (Xi i x-column) (Yi i y-column) (Zi i z-column))
      (do ((j (fx+ i 1) (fx+ j 1)))
          ((fx= j n-body))
        (let-thval
         ((massj j mass-column) (Xj j x-column) (Yj j y-column) (Zj j z-column))
         (let-getset
          ((VXi= i veloce-x) (VYi= i veloce-y) (VZi= i veloce-z)
           (VXj= j veloce-x) (VYj= j veloce-y) (VZj= j veloce-z))
          (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 massj mag))
            (VYi= fl- (fl* dy massj mag))
            (VZi= fl- (fl* dz massj mag))
            (VXj= fl+ (fl* dx massi mag))
            (VYj= fl+ (fl* dy massi mag))
            (VZj= fl+ (fl* dz massi mag))
           ))))))
   (do ((i 0 (fx+ i 1)))
       ((fx= i n-body))
     (let-thval
      ((VX i veloce-x) (VY i veloce-y) (VZ i veloce-z))
      (let-getset
       ((X= i x-column) (Y= i y-column) (Z= i z-column))
       (X= fl+ (fl* dt VX))
       (Y= fl+ (fl* dt VY))
       (Z= fl+ (fl* dt VZ))))))

  (do ((i 0 (fx+ i 1)) (E 0.0))
      ((fx= i n-body) E)
    (let-thval
     ((massi i mass-column) (VX i veloce-x) (VY i veloce-y) (VZ i veloce-z))
     (addn! E (fl* 0.5 massi (sum^2 VX VY VZ)))
     (do ((j (fx+ i 1) (fx+ j 1)))
         ((fx= j n-body))
       (let-thval
        ((Xi i x-column) (Yi i y-column) (Zi i z-column)
         (Xj j x-column) (Yj j y-column) (Zj j z-column) (massj j mass-column))
        (subn! E (fl/ (fl* massi massj)
                      (flsqrt (sum^2
                               (fl- Xi Xj)
                               (fl- Yi Yj)
                               (fl- Zi Zj)))))))))
  )

(define (main . args)
  (cond ((pair? args)
         (let* ((n (string->number (car args)))
                (st (current-second))
                (e (take-energy n))
                (et (current-second)))
           (print e " per " (- et st) " sec on " n "\n")))
        (else (print "Usage: gsc-colnbody <n-iterations>~%"))))