Difference between revisions of "Programming language shootout: n-body"
From Gambit wiki
(Add n-body alioth benchmark) |
m (specialize two inequalities to fixnums) |
||
Line 173: | Line 173: | ||
;; ------------------------------- | ;; ------------------------------- | ||
(define (main . args) | (define (main . args) | ||
− | (let ((n (if (< (length args) 1) | + | (let ((n (if (fx< (length args) 1) |
1 | 1 | ||
(string->number (car args)))) | (string->number (car args)))) | ||
Line 180: | Line 180: | ||
(display (roundto 9 (energy system))) (newline) | (display (roundto 9 (energy system))) (newline) | ||
(do ((i 1 (fx+ i 1))) | (do ((i 1 (fx+ i 1))) | ||
− | ((< n i)) | + | ((fx< n i)) |
(advance system 0.01)) | (advance system 0.01)) | ||
(display (roundto 9 (energy system))) (newline))) | (display (roundto 9 (energy system))) (newline))) |
Revision as of 20:59, 22 February 2008
This is a Gambit implementation of the n-body benchmark of the Computer Language Benchmarks Game.
The program
#!gsi-script ;; The Computer Language Benchmarks Game ;; http://shootout.alioth.debian.org/ ;; ;; Derived by Bradley Lucier from the Ikarus variant ;; Derived by Michael D. Adams from the MzScheme variant #| Correct output N = 1000 is -0.169075164 -0.169087605 |# (declare (standard-bindings)(extended-bindings)(block)(not safe)) (define-macro (unless test . body) `(if (not ,test) (begin ,@body))) ;;; Stupid boiler-plate for formatting floating point values (define (roundto digits n) (let* ([e (expt 10 digits)] [num (round (abs (* e (inexact->exact n))))] [str (number->string (modulo num e))]) (string-append (if (negative? n) "-" "") (number->string (quotient num e)) "." (make-string (- digits (string-length str)) #\0) str))) ;; ------------------------------ ;; define planetary masses, initial positions & velocity (define *pi* 3.141592653589793) (define *days-per-year* 365.24) (define *solar-mass* (* 4 *pi* *pi*)) (define (make-body x y z vx vy vz mass) (f64vector x y z vx vy vz mass)) (define (body-x body) (f64vector-ref body 0)) (define (body-x-set! body val) (f64vector-set! body 0 val)) (define (body-y body) (f64vector-ref body 1)) (define (body-y-set! body val) (f64vector-set! body 1 val)) (define (body-z body) (f64vector-ref body 2)) (define (body-z-set! body val) (f64vector-set! body 2 val)) (define (body-vx body) (f64vector-ref body 3)) (define (body-vx-set! body val) (f64vector-set! body 3 val)) (define (body-vy body) (f64vector-ref body 4)) (define (body-vy-set! body val) (f64vector-set! body 4 val)) (define (body-vz body) (f64vector-ref body 5)) (define (body-vz-set! body val) (f64vector-set! body 5 val)) ;; mass is immutable (define (body-mass body) (f64vector-ref body 6)) (define *sun* (make-body 0.0 0.0 0.0 0.0 0.0 0.0 *solar-mass*)) (define *jupiter* (make-body 4.84143144246472090 -1.16032004402742839 -1.03622044471123109e-1 (* 1.66007664274403694e-3 *days-per-year*) (* 7.69901118419740425e-3 *days-per-year*) (* -6.90460016972063023e-5 *days-per-year*) (* 9.54791938424326609e-4 *solar-mass*))) (define *saturn* (make-body 8.34336671824457987 4.12479856412430479 -4.03523417114321381e-1 (* -2.76742510726862411e-3 *days-per-year*) (* 4.99852801234917238e-3 *days-per-year*) (* 2.30417297573763929e-5 *days-per-year*) (* 2.85885980666130812e-4 *solar-mass*))) (define *uranus* (make-body 1.28943695621391310e1 -1.51111514016986312e1 -2.23307578892655734e-1 (* 2.96460137564761618e-03 *days-per-year*) (* 2.37847173959480950e-03 *days-per-year*) (* -2.96589568540237556e-05 *days-per-year*) (* 4.36624404335156298e-05 *solar-mass*))) (define *neptune* (make-body 1.53796971148509165e+01 -2.59193146099879641e+01 1.79258772950371181e-01 (* 2.68067772490389322e-03 *days-per-year*) (* 1.62824170038242295e-03 *days-per-year*) (* -9.51592254519715870e-05 *days-per-year*) (* 5.15138902046611451e-05 *solar-mass*))) ;; ------------------------------- (define (offset-momentum system) (let loop-i ((i system) (px 0.0) (py 0.0) (pz 0.0)) (if (null? i) (begin (body-vx-set! (car system) (fl/ (fl- px) *solar-mass*)) (body-vy-set! (car system) (fl/ (fl- py) *solar-mass*)) (body-vz-set! (car system) (fl/ (fl- pz) *solar-mass*))) (loop-i (cdr i) (fl+ px (fl* (body-vx (car i)) (body-mass (car i)))) (fl+ py (fl* (body-vy (car i)) (body-mass (car i)))) (fl+ pz (fl* (body-vz (car i)) (body-mass (car i)))))))) ;; ------------------------------- (define (energy system) (let loop-o ((o system) (e 0.0)) (if (null? o) e (let ([e (fl+ e (fl* 0.5 (body-mass (car o)) (fl+ (fl* (body-vx (car o)) (body-vx (car o))) (fl* (body-vy (car o)) (body-vy (car o))) (fl* (body-vz (car o)) (body-vz (car o))))))]) (let loop-i ((i (cdr o)) (e e)) (if (null? i) (loop-o (cdr o) e) (let* ((dx (fl- (body-x (car o)) (body-x (car i)))) (dy (fl- (body-y (car o)) (body-y (car i)))) (dz (fl- (body-z (car o)) (body-z (car i)))) (distance (flsqrt (fl+ (fl* dx dx) (fl* dy dy) (fl* dz dz))))) (let ([e (fl- e (fl/ (fl* (body-mass (car o)) (body-mass (car i))) distance))]) (loop-i (cdr i) e))))))))) ;; ------------------------------- (define (advance system dt) (let loop-o ((o system)) (unless (null? o) (let loop-i ((i (cdr o))) (unless (null? i) (let* ((o1 (car o)) (i1 (car i)) (dx (fl- (body-x o1) (body-x i1))) (dy (fl- (body-y o1) (body-y i1))) (dz (fl- (body-z o1) (body-z i1))) (distance (flsqrt (fl+ (fl* dx dx) (fl* dy dy) (fl* dz dz)))) (mag (fl/ dt (fl* distance distance distance))) (dxmag (fl* dx mag)) (dymag (fl* dy mag)) (dzmag (fl* dz mag)) (om (body-mass o1)) (im (body-mass i1))) (body-vx-set! o1 (fl- (body-vx o1) (fl* dxmag im))) (body-vy-set! o1 (fl- (body-vy o1) (fl* dymag im))) (body-vz-set! o1 (fl- (body-vz o1) (fl* dzmag im))) (body-vx-set! i1 (fl+ (body-vx i1) (fl* dxmag om))) (body-vy-set! i1 (fl+ (body-vy i1) (fl* dymag om))) (body-vz-set! i1 (fl+ (body-vz i1) (fl* dzmag om))) (loop-i (cdr i))))) (loop-o (cdr o)))) (let loop-o ((o system)) (unless (null? o) (let ([o1 (car o)]) (body-x-set! o1 (fl+ (body-x o1) (fl* dt (body-vx o1)))) (body-y-set! o1 (fl+ (body-y o1) (fl* dt (body-vy o1)))) (body-z-set! o1 (fl+ (body-z o1) (fl* dt (body-vz o1)))) (loop-o (cdr o)))))) ;; ------------------------------- (define (main . args) (let ((n (if (fx< (length args) 1) 1 (string->number (car args)))) (system (list *sun* *jupiter* *saturn* *uranus* *neptune*))) (offset-momentum system) (display (roundto 9 (energy system))) (newline) (do ((i 1 (fx+ i 1))) ((fx< n i)) (advance system 0.01)) (display (roundto 9 (energy system))) (newline)))
Compiling
gsc n-body
Running
gsi n-body 20000000