1
;; crunch-nbody.scm
;; NBODY exercise. CHICKEN6 git + Crunch
;; compile: $ csc6 -C "-DNDEBUG -O2 -march=native" crunch-colnbody.scm
;; run: $ ./crunch-colnbody <n-iterations>
;;
(declare (disable-interrupts))
(import (scheme base)
(scheme write)
(scheme time)
(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-syntax mkcol (syntax-rules () ((_) (make-f64vector n-body))))
(define mass-column (mkcol))
(define x-column (mkcol))
(define y-column (mkcol))
(define z-column (mkcol))
(define veloce-x (mkcol))
(define veloce-y (mkcol))
(define veloce-z (mkcol))
(define-syntax set-planet (syntax-rules () ((_ offset mass x y z vx vy vz)
(begin
(f64vector-set! mass-column offset mass)
(f64vector-set! x-column offset x)
(f64vector-set! y-column offset y)
(f64vector-set! z-column offset z)
(f64vector-set! veloce-x offset vx)
(f64vector-set! veloce-y offset vy)
(f64vector-set! veloce-z offset vz)
(add1 offset)))))
(define lp (set-planet 0 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-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 dt 0.01)
(define-syntax let-thval
(syntax-rules () ((_ ((a ix column) ...) body ...)
(let-syntax ((a (syntax-rules () ((a) (f64vector-ref column ix)))) ...)
body ...))))
(define-syntax let-getset
(syntax-rules () ((_ ((a ix column) ...) body ...)
(let-syntax ((a (syntax-rules () ((a op arg)
(f64vector-set! column ix
(op (f64vector-ref column ix)
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 ^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)
;;(define st (current-jiffy))
(dotimes
n
(do ((i 0 (add1 i)))
((= i n-body))
(do ((j (add1 i) (add1 j)))
((= j n-body))
(let-thval
((mass_i i mass-column) (Xi i x-column) (Yi i y-column) (Zi i z-column)
(mass_j j mass-column) (Xj j x-column) (Yj j y-column) (Zj j z-column))
(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 veloce-x) (VYi= i veloce-y) (VZi= i veloce-z)
(VXj= j veloce-x) (VYj= j veloce-y) (VZj= j veloce-z))
(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 (add1 i)))
((= 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= + (* dt (VX)))
(Y= + (* dt (VY)))
(Z= + (* dt (VZ)))))))
(do ((i 0 (add1 i)))
((= i n-body))
(let-thval
((mass_i i mass-column)
(VX i veloce-x)
(VY i veloce-y)
(VZ i veloce-z))
(addn! E (* 0.5 (* (mass_i) (sum^2 (VX) (VY) (VZ)))))
(do ((j (add1 i) (add1 j)))
((= 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) (mass_j j mass-column))
(subn! E (/ (* (mass_i) (mass_j))
(sqrt (sum^2
(- (Xi) (Xj))
(- (Yi) (Yj))
(- (Zi) (Zj))))))
)
)))
;;(define et (current-jiffy))
;;(display (- et st))
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-colnbody <n-iterations>"))))
(cond-expand
(compiling
(main (argv)))
(else))
For immediate assistance, please email our customer support: [email protected]