1
;; crunch-nbody.scm
;; NBODY exercise. CHICKEN6 git + Crunch
;; preinstall: $ chicken-install6 crunch
;; compile: $ csc6 -C "-O3 -march=native" crunch-nbody.scm
;; run: $ ./crunch-nbody <n-iterations>
;;
(declare (disable-interrupts))
(import (scheme base)
(scheme write)
(chicken base)
(chicken fixnum)
(chicken flonum)
(chicken time)
(chicken process-context)
(chicken number-vector)
(chicken foreign)
(crunch)
(crunch c)
(crunch declarations))
(crunch
(: take-energy (procedure (integer) float))
(define (take-energy n)
(define PI (* 4.0 (fpatan 1.)))
(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 N-BODY 5)
(define body-size 7)
(define bodies-len (* N-BODY body-size))
(define bodies (make-f64vector bodies-len 0.0))
;;(: set-planet (procedure ((numvector f64) integer float float float float float float float) integer))
(define-syntax set-planet (syntax-rules () ((_ offset mass x y z vx vy vz)
(begin
(f64vector-set! bodies (+ offset 0) mass)
(f64vector-set! bodies (+ offset 1) x)
(f64vector-set! bodies (+ offset 2) y)
(f64vector-set! bodies (+ offset 3) z)
(f64vector-set! bodies (+ offset 4) vx)
(f64vector-set! bodies (+ offset 5) vy)
(f64vector-set! bodies (+ offset 6) vz)
(+ offset body-size)))))
(define lp 0)
(set! lp (set-planet lp SOLAR-MASS 0.0 0.0 0.0 0.0 0.0 0.0))
(set! lp (set-planet lp (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)))
(set! lp (set-planet lp (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)))
(set! lp (set-planet lp (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)))
(set! lp (set-planet lp (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-syntax let-attr
(syntax-rules () ((_ ((a ix offset) ...) body ...)
(let-syntax ((a (syntax-rules ()((a) (f64vector-ref bodies (+ ix offset))))) ...)
body ...))))
(define-syntax let-getset
(syntax-rules () ((_ ((a ix offset) ...) body ...)
(let-syntax ((a (syntax-rules () ((a op arg)
(f64vector-set! bodies (+ ix offset)
(op (f64vector-ref bodies (+ ix offset))
arg))))) ...)
body ...))))
(define-syntax subn! (syntax-rules () ((_ what withs) (set! what (- what withs)))))
(define-syntax addn! (syntax-rules () ((_ what withs) (set! what (+ what withs)))))
(define-syntax subnv! (syntax-rules () ((_ ix offset val deltaval) (f64vector-set! bodies (+ ix offset) (- val deltaval)))))
(define-syntax addnv! (syntax-rules () ((_ ix offset val deltaval) (f64vector-set! bodies (+ ix offset) (+ val deltaval)))))
(define-syntax ^2 (syntax-rules () ((_ x) (* x x))))
(define-syntax sum^2 (syntax-rules () ((_ x y z) (+ (^2 x) (+ (^2 y) (^2 z))))))
(define-syntax dsqt (syntax-rules () ((_ v) (let ((v_ v)) (* v_ (sqrt v_))))))
(define-syntax dotimes (syntax-rules () ((_ n body ...) (do ((i_ n (- i_ 1))) ((= i_ 0) ) body ...))))
(define E 0.0)
(dotimes
n
(do ((i 0 (+ i body-size)))
((= i bodies-len))
(do ((j (+ i body-size) (+ j body-size)))
((= j bodies-len))
(let-attr
((mass_i i 0) (Xi i 1) (Yi i 2) (Zi i 3)
(mass_j j 0) (Xj j 1) (Yj j 2) (Zj j 3))
(let* ((dx (- (Xi) (Xj)))
(dy (- (Yi) (Yj)))
(dz (- (Zi) (Zj)))
(mag (/ dt (dsqt (sum^2 dx dy dz))))
(mag_j (* (mass_j) mag))
(mag_i (* (mass_i) mag)))
(let-getset
((VXi= i 4) (VYi= i 5) (VZi= i 6)
(VXj= j 4) (VYj= j 5) (VZj= j 6))
(VXi= - (* dx mag_j))
(VYi= - (* dy mag_j))
(VZi= - (* dz mag_j))
(VXj= + (* dx mag_i))
(VYj= + (* dy mag_i))
(VZj= + (* dz mag_i))
)))))
(do ((i 0 (+ i body-size)))
((= i bodies-len))
(let-attr
((VX i 4) (VY i 5) (VZ i 6))
(let-getset
((X= i 1) (Y= i 2) (Z= i 3))
(X= + (* dt (VX)))
(Y= + (* dt (VY)))
(Z= + (* dt (VZ)))))))
(do ((i 0 (+ i body-size)))
((= i bodies-len))
(let-attr
((mass_i i 0)
(VX i 4)
(VY i 5)
(VZ i 6))
(addn! E (* 0.5 (* (mass_i) (sum^2 (VX) (VY) (VZ)))))
(do ((j (+ i body-size) (+ j body-size)))
((= j bodies-len))
(let-attr
((Xi i 1) (Yi i 2) (Zi i 3)
(Xj j 1) (Yj j 2) (Zj j 3) (mass_j j 0))
(subn! E (/ (* (mass_i) (mass_j))
(sqrt (sum^2
(- (Xi) (Xj))
(- (Yi) (Yj))
(- (Zi) (Zj))))))
)
)))
E
)
);;end-of-crunch
(define (main args)
(cond ((and (pair? args) (pair? (cdr args)))
(let* ((n (string->number (cadr args)))
(st (current-process-milliseconds))
(e (take-energy n))
(et (current-process-milliseconds)))
(print e " per " (/ (- et st) 1000.) " sec on " n "\n"))
(exit 0))
(else (print "Usage: crunch-nbody <n-iterations>"))))
(cond-expand
(compiling
(main (argv)))
(else))
For immediate assistance, please email our customer support: [email protected]