Added Mandelbrot and N-body bechmarking examples.

This commit is contained in:
Thomas Dendale 2017-12-19 17:18:44 +01:00
parent 3656cf274f
commit 51e2eb414c
2 changed files with 241 additions and 0 deletions

View File

@ -0,0 +1,77 @@
;;; Taken from
;;; http://benchmarksgame.alioth.debian.org/u64q/program.php?test=mandelbrot&lang=gcc&id=2
(add-cflag "-Wall -pipe -O3 -fomit-frame-pointer -march=native")
(use Double)
(use Int)
(use String)
(use IO)
(def w 16000.0)
(def h 16000.0)
(system-include "stdio.h")
(register-type FILE)
(register stdout FILE)
(register putc (Fn [Int FILE] ()))
(defdynamic string-append-impl [xs]
(if (= (count xs) 1)
(car xs)
(if (= (count xs) 2)
(list 'append (car xs) (car (cdr xs)))
(list 'append (car xs) (string-append-impl (cdr xs))))))
(defmacro string-append [:rest xs]
(string-append-impl xs))
(defn main []
(do
(print &(string-append @"P4\n" (str w) @" " (str h) @"\n"))
(let [iter 50
limit 2.0
byte_acc 0
bit_num 0]
(for [y 0 (to-int h)]
(for [x 0 (to-int w)]
(let [Zr 0.0
Zi 0.0
Tr 0.0
Ti 0.0
Cr (- (/ (* 2.0 (from-int x)) w) 1.5)
Ci (- (/ (* 2.0 (from-int y)) w) 1.0)]
(do
(let [i 0]
(while (and (< i iter)
(<= (+ Tr Ti) (* limit limit)))
(do
(set! &Zi (+ (* 2.0 (* Zr Zi)) Ci))
(set! &Zr (+ (- Tr Ti) Cr))
(set! &Tr (* Zr Zr))
(set! &Ti (* Zi Zi))
(set! &i (inc i)))))
(set! &byte_acc (* 2 byte_acc))
(when (<= (+ Tr Ti) (* limit limit))
(set! &byte_acc (bit-or byte_acc 1)))
(set! &bit_num (inc bit_num))
(if (= bit_num 8)
(do
(putc byte_acc stdout)
(set! &byte_acc 0)
(set! &bit_num 0))
(if (= x (- (to-int w) 1))
(do
(set! &byte_acc
(bit-shift-left byte_acc (- 8 (mod (to-int w) 8))))
(putc byte_acc stdout)
(set! &byte_acc 0)
(set! &bit_num 0))
())))))))))
;;; (build)

View File

@ -0,0 +1,164 @@
;;; Taken from
;;; http://benchmarksgame.alioth.debian.org/u64q/program.php?test=nbody&lang=gcc&id=1
(add-cflag "-Wall -pipe -O3 -fomit-frame-pointer -march=native -mfpmath=sse -msse3")
(use IO)
(use Double)
(use Array)
(use Int)
;;; pow(double, double) is 5 times slower for (pow x 3)
; int Int_pow(int x, int y) {
; int r = 1;
; while (y) {
; if (y & 1)
; r *= x;
; y >>= 1;
; x *= x;
; }
; return r;
; }
(defn ipow [x y]
(let-do [r 1.0]
(while (/= y 0)
(do
(when (/= (bit-and y 1) 0)
(set! &r (* r x)))
(set! &y (/ y 2))
(set! &x (* x x))))
r))
(def n 50000000)
(defn solar_mass [] (* (* 4.0 pi) pi)) ;; I cannot #DEFINE, so solar_mass cannot be a `def`
(def days_per_year 365.24)
(deftype Planet
[ x Double y Double z Double
vx Double vy Double vz Double
mass Double ])
; for eg. 1.03e-03
(defn e [d exp]
(* d (pow 10.0 (from-int exp))))
; currying would make this a single line
; + lambdas would make this unneeded
(defn reduce-px [t b] (+ @t (* (Planet.vx b) (Planet.mass b))))
(defn reduce-py [t b] (+ @t (* (Planet.vy b) (Planet.mass b))))
(defn reduce-pz [t b] (+ @t (* (Planet.vz b) (Planet.mass b))))
(defn offset_momentum [bodies]
(let [b (first bodies)
px (reduce reduce-px 0.0 bodies)
py (reduce reduce-py 0.0 bodies)
pz (reduce reduce-pz 0.0 bodies)]
(aset! bodies 0
(=> b
(Planet.set-vx (/ (neg px) (solar_mass)))
(Planet.set-vy (/ (neg py) (solar_mass)))
(Planet.set-vz (/ (neg pz) (solar_mass)))))))
(defn energy [bodies]
(let-do [e 0.0]
(for [i 0 (count bodies)]
(let-do [b (nth bodies i)]
(set! &e (+ e (* 0.5 (* (Planet.mass b)
(+ (ipow (Planet.vx b) 2)
(+ (ipow (Planet.vy b) 2)
(ipow (Planet.vz b) 2)))))))
(for [j (+ i 1) (count bodies)]
(let [b2 (nth bodies j)
dx (- (Planet.x b) (Planet.x b2))
dy (- (Planet.y b) (Planet.y b2))
dz (- (Planet.z b) (Planet.z b2))
dist (sqrt (+ (ipow dx 2) (+ (ipow dy 2) (ipow dz 2))))]
(set! &e (- e (/ (* (Planet.mass b) (Planet.mass b2)) dist)))))))
e))
(defn update-planet [i bodies dt]
(let [b (nth bodies i)]
(aset! bodies i
(=> @b
(Planet.set-x (+ (Planet.x b) (* dt (Planet.vx b))))
(Planet.set-y (+ (Planet.y b) (* dt (Planet.vy b))))
(Planet.set-z (+ (Planet.z b) (* dt (Planet.vz b))))))))
(defn advance [bodies dt]
(do
(for [i 0 (count bodies)]
(let [b (nth bodies i)]
(for [j (+ i 1) (count bodies)]
(let-do [b2 (nth bodies j)
dx (- (Planet.x b) (Planet.x b2))
dy (- (Planet.y b) (Planet.y b2))
dz (- (Planet.z b) (Planet.z b2))
dist (sqrt (+ (ipow dx 2) (+ (ipow dy 2) (ipow dz 2))))
mag (/ dt (ipow dist 3))]
(aset! bodies i
(=> @b
(Planet.set-vx (- (Planet.vx b) (* dx (* (Planet.mass b2) mag))))
(Planet.set-vy (- (Planet.vy b) (* dy (* (Planet.mass b2) mag))))
(Planet.set-vz (- (Planet.vz b) (* dz (* (Planet.mass b2) mag))))))
(aset! bodies j
(=> @b2
(Planet.set-vx (+ (Planet.vx b2) (* dx (* (Planet.mass b) mag))))
(Planet.set-vy (+ (Planet.vy b2) (* dy (* (Planet.mass b) mag))))
(Planet.set-vz (+ (Planet.vz b2) (* dz (* (Planet.mass b) mag))))))))))
(for [i 0 (count bodies)]
(update-planet i bodies dt))))
; defining bodies outside main (`def`) results in c compile error:
; - solar mass is a function
; - i cannot do eg. 1.01e-03
(defn main []
(let-do [bodies [
(Planet.init 0.0 0.0 0.0 0.0 0.0 0.0 (solar_mass))
(Planet.init ; jupiter
(e 4.84143144246472090 0)
(e -1.16032004402742839 0)
(e -1.03622044471123109 -1)
(* (e 1.66007664274403694 -3) days_per_year)
(* (e 7.69901118419740425 -3) days_per_year)
(* (e -6.90460016972063023 -5) days_per_year)
(* (e 9.54791938424326609 -4) (solar_mass)))
(Planet.init ; saturn
(e 8.34336671824457987 0)
(e 4.12479856412430479 0)
(e -4.03523417114321381 -1)
(* (e -2.76742510726862411 -3) days_per_year)
(* (e 4.99852801234917238 -3) days_per_year)
(* (e 2.30417297573763929 -5) days_per_year)
(* (e 2.85885980666130812 -4) (solar_mass)))
(Planet.init ; uranus
(e 1.28943695621391310 1)
(e -1.51111514016986312 1)
(e -2.23307578892655734 -1)
(* (e 2.96460137564761618 -3) days_per_year)
(* (e 2.37847173959480950 -3) days_per_year)
(* (e -2.96589568540237556 -5) days_per_year)
(* (e 4.36624404335156298 -5) (solar_mass)))
(Planet.init ; neptune
(e 1.53796971148509165 1)
(e -2.59193146099879641 1)
(e 1.79258772950371181 -1)
(* (e 2.68067772490389322 -3) days_per_year)
(* (e 1.62824170038242295 -3) days_per_year)
(* (e -9.51592254519715870 -5) days_per_year)
(* (e 5.15138902046611451 -5) (solar_mass)))
]]
(offset_momentum &bodies)
(println &(str (energy &bodies))) ; printf %.9f
(for [i 0 n]
(do
(advance &bodies 0.01)))
(println &(str (energy &bodies))))) ;