ikarus/benchmarks/r6rs-benchmarks/nbody.ss

1447 lines
52 KiB
Scheme

; This benchmark was obtained from Andrew Wright.
; 970215 / wdc Changed {box, unbox, set-box!} to {list, car, set-car!},
; flushed #% prefixes, defined void,
; and added nbody-benchmark.
; 981116 / wdc Replaced nbody-benchmark by main, added apply:+.
(library (r6rs-benchmarks nbody)
(export main)
(import
(r6rs)
(r6rs arithmetic flonums)
(r6rs mutable-pairs)
(r6rs-benchmarks))
(define void
(let ((invisible (string->symbol "")))
(lambda () invisible)))
(define (apply:+ xs)
(do ((result 0.0 (fl+ result (car xs)))
(xs xs (cdr xs)))
((null? xs) result)))
;; code is slightly broken...
(define vect (vector))
; minimal standard random number generator
; 32 bit integer version
; cacm 31 10, oct 88
;
(define *seed* (list 1))
(define (srand seed)
(set-car! *seed* seed))
(define (rand)
(let* ((hi (quotient (car *seed*) 127773))
(lo (modulo (car *seed*) 127773))
(test (- (* 16807 lo) (* 2836 hi))))
(if (> test 0)
(set-car! *seed* test)
(set-car! *seed* (+ test 2147483647)))
(car *seed*)))
;; return a random number in the interval [0,n)
(define random
(lambda (n)
(modulo (abs (rand)) n)))
(define (array-ref a . indices)
(let loop ((a a) (indices indices))
(if (null? indices)
a
(loop (vector-ref a (car indices)) (cdr indices)))))
(define (atan0 y x)
(if (and (= x y) (= x 0)) 0 (atan y x)))
;change this to desired precision
;measured in order of expansions calculated
(define precision 10)
;;; =========================================================
;;; Algorithm
;;; =========================================================
(define (cartesian-algorithm tree)
(up! tree
cartesian-make-multipole-expansion
cartesian-multipole-shift
cartesian-expansion-sum)
(down! tree
cartesian-multipole-to-local-convert
cartesian-local-shift
cartesian-eval-local-expansion
cartesian-expansion-sum
cartesian-zero-expansion))
(define (spherical-algorithm tree)
(up! tree
spherical-make-multipole-expansion
spherical-multipole-shift
spherical-expansion-sum)
(down! tree
spherical-multipole-to-local-convert
spherical-local-shift
spherical-eval-local-expansion
spherical-expansion-sum
spherical-zero-expansion))
;;; The upward path in the algorithm calculates
;;; multipole expansions at every node.
(define (up! tree
make-multipole-expansion
multipole-shift
expansion-sum)
(let loop ((node (tree-body tree)))
(define center (node-center node))
(if (leaf-node? node)
(let ((multipole-expansion
(expansion-sum
(map (lambda (particle)
(make-multipole-expansion
(pt- (particle-position particle) center)
(particle-strength particle)))
(node-particles node)))))
(set-node-multipole-expansion! node multipole-expansion)
(cons center multipole-expansion))
(let ((multipole-expansion
(expansion-sum
(map (lambda (child)
(define center-and-expansion (loop child))
(multipole-shift (pt- (car center-and-expansion)
center)
(cdr center-and-expansion)))
(node-children node)))))
(set-node-multipole-expansion! node multipole-expansion)
(cons center multipole-expansion)))))
;;; Downward path of the algorithm which calculates local expansionss
;;; at every node and accelerations and potentials at each particle.
(define (down! tree
multipole-to-local-convert
local-shift
eval-local-expansion
expansion-sum
zero-expansion)
(let loop ((node (tree-body tree))
(parent-local-expansion (zero-expansion))
(parent-center (node-center (tree-body tree))))
(let* ((center (node-center node))
(interactive-sum
(expansion-sum
(map (lambda (interactive)
(multipole-to-local-convert
(pt- (node-center interactive) center)
(node-multipole-expansion interactive)))
(node-interactive-field node))))
(local-expansion
(expansion-sum (list (local-shift (pt- center parent-center)
parent-local-expansion)
interactive-sum))))
(if (leaf-node? node)
(eval-potentials-and-accelerations
node
local-expansion
eval-local-expansion)
(for-each (lambda (child) (loop child local-expansion center))
(node-children node))))))
(define (eval-potentials-and-accelerations
node
local-expansion
eval-local-expansion)
(let ((center (node-center node))
(near-field (apply append (map node-particles
(node-near-field node)))))
(for-each (lambda (particle)
(let* ((pos (particle-position particle))
(far-field-accel-and-poten
(eval-local-expansion (pt- pos center)
local-expansion)))
(set-particle-acceleration!
particle
(pt+ (car far-field-accel-and-poten)
(sum-vectors
(map (lambda (near)
(direct-accel (pt- (particle-position near)
pos)
(particle-strength near)))
(nfilter near-field
(lambda (near)
(not (eq? near particle))))))))
(set-particle-potential!
particle
(+ (cdr far-field-accel-and-poten)
(apply:+ (map (lambda (near)
(direct-poten
(pt- (particle-position near) pos)
(particle-strength near)))
(nfilter near-field
(lambda (near)
(not
(eq? near particle))))))))))
(node-particles node))))
;;; ================================================================
;;; Expansion Theorems
;;; ================================================================
(define (cartesian-make-multipole-expansion pt strength)
(let ((min-x (- (pt-x pt)))
(min-y (- (pt-y pt)))
(min-z (- (pt-z pt))))
(make-cartesian-expansion
(lambda (i j k) (* strength
(expt min-x i)
(expt min-y j)
(expt min-z k)
(n1/prod-fac i j k))))))
(define (spherical-make-multipole-expansion pt strength)
(let ((r (pt-r pt))
(theta (pt-theta pt))
(phi (pt-phi pt)))
(make-spherical-expansion
(lambda (l m)
(* strength
(expt r l)
(eval-spherical-harmonic l (- m) theta phi))))))
;;; ==================================================================
;;; Shifting lemmas
;;; ==================================================================
;;; Shift multipole expansion
(define (cartesian-multipole-shift pt multipole-expansion)
(let ((pt-expansion (cartesian-make-multipole-expansion pt 1)))
(make-cartesian-expansion
(lambda (i j k)
(sum-3d i j k
(lambda (l m n)
(* (array-ref multipole-expansion l m n)
(array-ref pt-expansion (- i l) (- j m) (- k n)))))))))
(define (spherical-multipole-shift pt multipole-expansion)
(let ((r (pt-r pt))
(theta (pt-theta pt))
(phi (pt-phi pt)))
(letrec ((foo (lambda (a b)
(if (< (* a b) 0)
(expt -1 (min (abs a) (abs b)))
1)))
(bar (lambda (a b)
(/ (expt -1 a)
(sqrt (* (fac (- a b)) (fac (+ a b))))))))
(let ((pt-expansion (make-spherical-expansion
(lambda (l m)
(* (eval-spherical-harmonic l m theta phi)
(bar l m)
(expt r l))))))
(make-spherical-expansion
(lambda (j k)
(sum-2d j
(lambda (l m)
(if (> (abs (- k m)) (- j l))
0
(* (spherical-ref multipole-expansion (- j l) (- k m))
(foo m (- k m))
(bar (- j l) (- k m))
(spherical-ref pt-expansion l (- m))
(/ (bar j k))))))))))))
;;; Convert multipole to local
(define (cartesian-multipole-to-local-convert pt multipole-expansion)
(define pt-expansion
(let* ((n1/radius (/ (pt-r vect)))
(n2cosines (pt-scalar* (* 2 n1/radius) vect))
(x (pt-x n2cosines))
(y (pt-y n2cosines))
(z (pt-z n2cosines)))
(make-cartesian-expansion
(lambda (i j k)
(define ijk (+ i j k))
(* (expt -1 ijk)
(expt n1/radius (+ 1 ijk))
(prod-fac i j k)
(sum-3d (/ i 2) (/ j 2) (/ k 2)
(lambda (l m n)
(* (fac-1 (- ijk l m n))
(n1/prod-fac l m n)
(n1/prod-fac (- i (* 2 l))
(- j (* 2 m))
(- k (* 2 n)))
(expt x (- i (* 2 l)))
(expt y (- j (* 2 m)))
(expt z (- k (* 2 n)))))))))))
(make-cartesian-expansion
(lambda (i j k)
(sum2-3d (- precision i j k)
(lambda (l m n)
(* (array-ref multipole-expansion l m n)
(array-ref pt-expansion (+ i l) (+ j m) (+ k n))))))))
(define (spherical-multipole-to-local-convert pt multipole-expansion)
(let* ((r (pt-r pt))
(theta (pt-theta pt))
(phi (pt-phi pt))
(foo (lambda (a b c)
(* (expt -1 b)
(if (> (* a c) 0)
(expt -1 (min (abs a) (abs c)))
1))))
(bar (lambda (a b)
(/ (expt -1 a)
(sqrt (* (fac (- a b)) (fac (+ a b)))))))
(pt-expansion (make-spherical-expansion
(lambda (l m)
(/ (eval-spherical-harmonic l m theta phi)
(bar l m)
(expt r (+ 1 l)))))))
(make-spherical-expansion
(lambda (j k)
(* (bar j k)
(sum-2d (- precision j 1)
(lambda (l m)
(* (spherical-ref multipole-expansion l m)
(foo k l m)
(bar l m)
(spherical-ref pt-expansion (+ j l) (- m k))))))))))
;;; Shift local expansion
(define (cartesian-local-shift pt local-expansion)
(let* ((x (pt-x pt))
(y (pt-y pt))
(z (pt-z pt))
(expts (make-cartesian-expansion
(lambda (l m n) (* (expt x l)
(expt y m)
(expt z n))))))
(make-cartesian-expansion
(lambda (i j k)
(sum2-3d (- precision i j k)
(lambda (l m n)
(* (array-ref local-expansion (+ i l) (+ j m) (+ k n))
(array-ref expts l m n)
(n1/prod-fac l m n))))))))
(define (spherical-local-shift pt local-expansion)
(let* ((pt (pt- (make-pt 0 0 0) pt))
(r (pt-r pt))
(theta (pt-theta pt))
(phi (pt-phi pt))
(foo (lambda (a b c)
(* (expt -1 a)
(expt -1 (/ (+ (abs (- b c))
(abs b)
(- (abs c)))
2)))))
(bar (lambda (a b)
(/ (expt -1 a)
(sqrt (* (fac (- a b))
(fac (+ a b)))))))
(stuff (make-spherical-expansion
(lambda (l m)
(* (eval-spherical-harmonic l m theta phi)
(expt r l))))))
(make-spherical-expansion
(lambda (j k)
(sum2-2d j (lambda (l m)
(if (> (abs (- m k)) (- l j))
0
(* (spherical-ref local-expansion l m)
(bar (- l j) (- m k))
(bar j k)
(spherical-ref stuff (- l j) (- m k))
(/ (foo (- l j) (- m k) m) (bar l m))))))))))
;;; Evaluate the resulting local expansion at point pt.
(define (cartesian-eval-local-expansion pt local-expansion)
(let* ((x (pt-x pt))
(y (pt-y pt))
(z (pt-z pt))
(local-expansion
(make-cartesian-expansion
(lambda (i j k)
(* (array-ref local-expansion i j k)
(n1/prod-fac i j k)))))
(expts (make-cartesian-expansion
(lambda (i j k)
(* (expt x i) (expt y j) (expt z k))))))
(cons (make-pt
(sum2-3d (- precision 1)
(lambda (l m n)
(* (array-ref expts l m n)
(+ 1 l)
(array-ref local-expansion (+ 1 l) m n))))
(sum2-3d (- precision 1)
(lambda (l m n)
(* (array-ref expts l m n)
(+ 1 m)
(array-ref local-expansion l (+ 1 m) n))))
(sum2-3d (- precision 1)
(lambda (l m n)
(* (array-ref expts l m n)
(+ 1 n)
(array-ref local-expansion l m (+ 1 n))))))
(sum2-3d precision
(lambda (l m n)
(* (array-ref expts l m n)
(array-ref local-expansion l m n)))))))
(define (spherical-eval-local-expansion pt local-expansion)
(let* ((r (pt-r pt))
(x (pt-x pt))
(y (pt-y pt))
(z (pt-z pt))
(rho-sq (+ (* x x) (* y y)))
(theta (pt-theta pt))
(phi (pt-phi pt))
(r-deriv
(real-part
(sum-2d (- precision 1)
(lambda (l m)
(* (spherical-ref local-expansion l m)
(expt r (- l 1))
l
(eval-spherical-harmonic l m theta phi))))))
(theta-deriv
(real-part
(sum-2d (- precision 1)
(lambda (l m)
(* (spherical-ref local-expansion l m)
(expt r l)
(eval-spher-harm-theta-deriv l m theta phi))))))
(phi-deriv
(real-part
(sum-2d (- precision 1)
(lambda (l m)
(* (spherical-ref local-expansion l m)
(expt r l)
(eval-spher-harm-phi-deriv l m theta phi)))))))
(cons (make-pt
(+ (* r-deriv (/ x r))
(* theta-deriv (/ x (sqrt rho-sq) (+ z (/ rho-sq z))))
(* phi-deriv -1 (/ y rho-sq)))
(+ (* r-deriv (/ y r))
(* theta-deriv (/ y (sqrt rho-sq) (+ z (/ rho-sq z))))
(/ phi-deriv x (+ 1 (/ (* y y) x x))))
(+ (* r-deriv (/ z r))
(* theta-deriv (/ -1 r r) (sqrt rho-sq))))
(real-part (sum-2d (- precision 1)
(lambda (l m)
(* (spherical-ref local-expansion l m)
(eval-spherical-harmonic l m theta phi)
(expt r l))))))))
;;; Direct calculation of acceleration and potential
(define (direct-accel pt strength)
(pt-scalar* (/ strength (expt (pt-r pt) 3))
pt))
(define (direct-poten pt strength)
(/ strength (pt-r pt)))
;;; =================================================================
;;; TREES NODES PARTICLES and POINTS
;;; =================================================================
(begin (begin (begin (define make-raw-tree
(lambda (tree-1 tree-2 tree-3)
(vector '<tree> tree-1 tree-2 tree-3)))
(define tree?
(lambda (obj)
(if (vector? obj)
(if (= (vector-length obj) 4)
(eq? (vector-ref obj 0) '<tree>)
#f)
#f)))
(define tree-1
(lambda (obj)
(if (tree? obj)
(void)
(error
'tree-1
"~s is not a ~s"
obj
'<tree>))
(vector-ref obj 1)))
(define tree-2
(lambda (obj)
(if (tree? obj)
(void)
(error
'tree-2
"~s is not a ~s"
obj
'<tree>))
(vector-ref obj 2)))
(define tree-3
(lambda (obj)
(if (tree? obj)
(void)
(error
'tree-3
"~s is not a ~s"
obj
'<tree>))
(vector-ref obj 3)))
(define set-tree-1!
(lambda (obj newval)
(if (tree? obj)
(void)
(error
'set-tree-1!
"~s is not a ~s"
obj
'<tree>))
(vector-set! obj 1 newval)))
(define set-tree-2!
(lambda (obj newval)
(if (tree? obj)
(void)
(error
'set-tree-2!
"~s is not a ~s"
obj
'<tree>))
(vector-set! obj 2 newval)))
(define set-tree-3!
(lambda (obj newval)
(if (tree? obj)
(void)
(error
'set-tree-3!
"~s is not a ~s"
obj
'<tree>))
(vector-set! obj 3 newval))))
(define make-tree
(lambda (body low-left-front-vertex up-right-back-vertex)
((lambda ()
(make-raw-tree
body
low-left-front-vertex
up-right-back-vertex)))))
(define tree-body tree-1)
(define tree-low-left-front-vertex tree-2)
(define tree-up-right-back-vertex tree-3)
(define set-tree-body! set-tree-1!)
(define set-tree-low-left-front-vertex! set-tree-2!)
(define set-tree-up-right-back-vertex! set-tree-3!))
(begin (begin (define make-raw-node
(lambda (node-1
node-2
node-3
node-4
node-5
node-6
node-7
node-8)
(vector
'<node>
node-1
node-2
node-3
node-4
node-5
node-6
node-7
node-8)))
(define node?
(lambda (obj)
(if (vector? obj)
(if (= (vector-length obj) 9)
(eq? (vector-ref obj 0) '<node>)
#f)
#f)))
(define node-1
(lambda (obj)
(if (node? obj)
(void)
(error
'node-1
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 1)))
(define node-2
(lambda (obj)
(if (node? obj)
(void)
(error
'node-2
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 2)))
(define node-3
(lambda (obj)
(if (node? obj)
(void)
(error
'node-3
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 3)))
(define node-4
(lambda (obj)
(if (node? obj)
(void)
(error
'node-4
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 4)))
(define node-5
(lambda (obj)
(if (node? obj)
(void)
(error
'node-5
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 5)))
(define node-6
(lambda (obj)
(if (node? obj)
(void)
(error
'node-6
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 6)))
(define node-7
(lambda (obj)
(if (node? obj)
(void)
(error
'node-7
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 7)))
(define node-8
(lambda (obj)
(if (node? obj)
(void)
(error
'node-8
"~s is not a ~s"
obj
'<node>))
(vector-ref obj 8)))
(define set-node-1!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-1!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 1 newval)))
(define set-node-2!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-2!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 2 newval)))
(define set-node-3!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-3!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 3 newval)))
(define set-node-4!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-4!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 4 newval)))
(define set-node-5!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-5!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 5 newval)))
(define set-node-6!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-6!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 6 newval)))
(define set-node-7!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-7!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 7 newval)))
(define set-node-8!
(lambda (obj newval)
(if (node? obj)
(void)
(error
'set-node-8!
"~s is not a ~s"
obj
'<node>))
(vector-set! obj 8 newval))))
(define make-node
(lambda (center
low-left-front-vertex
up-right-back-vertex
children
particles
multipole-expansion
near-field
interactive-field)
((lambda ()
(make-raw-node
center
low-left-front-vertex
up-right-back-vertex
children
particles
multipole-expansion
near-field
interactive-field)))))
(define node-center node-1)
(define node-low-left-front-vertex node-2)
(define node-up-right-back-vertex node-3)
(define node-children node-4)
(define node-particles node-5)
(define node-multipole-expansion node-6)
(define node-near-field node-7)
(define node-interactive-field node-8)
(define set-node-center! set-node-1!)
(define set-node-low-left-front-vertex! set-node-2!)
(define set-node-up-right-back-vertex! set-node-3!)
(define set-node-children! set-node-4!)
(define set-node-particles! set-node-5!)
(define set-node-multipole-expansion! set-node-6!)
(define set-node-near-field! set-node-7!)
(define set-node-interactive-field! set-node-8!))
(define leaf-node?
(lambda (node) (null? (node-children node))))
(begin (begin (define make-raw-particle
(lambda (particle-1
particle-2
particle-3
particle-4
particle-5
particle-6)
(vector
'<particle>
particle-1
particle-2
particle-3
particle-4
particle-5
particle-6)))
(define particle?
(lambda (obj)
(if (vector? obj)
(if (= (vector-length obj) 7)
(eq? (vector-ref obj 0) '<particle>)
#f)
#f)))
(define particle-1
(lambda (obj)
(if (particle? obj)
(void)
(error
'particle-1
"~s is not a ~s"
obj
'<particle>))
(vector-ref obj 1)))
(define particle-2
(lambda (obj)
(if (particle? obj)
(void)
(error
'particle-2
"~s is not a ~s"
obj
'<particle>))
(vector-ref obj 2)))
(define particle-3
(lambda (obj)
(if (particle? obj)
(void)
(error
'particle-3
"~s is not a ~s"
obj
'<particle>))
(vector-ref obj 3)))
(define particle-4
(lambda (obj)
(if (particle? obj)
(void)
(error
'particle-4
"~s is not a ~s"
obj
'<particle>))
(vector-ref obj 4)))
(define particle-5
(lambda (obj)
(if (particle? obj)
(void)
(error
'particle-5
"~s is not a ~s"
obj
'<particle>))
(vector-ref obj 5)))
(define particle-6
(lambda (obj)
(if (particle? obj)
(void)
(error
'particle-6
"~s is not a ~s"
obj
'<particle>))
(vector-ref obj 6)))
(define set-particle-1!
(lambda (obj newval)
(if (particle? obj)
(void)
(error
'set-particle-1!
"~s is not a ~s"
obj
'<particle>))
(vector-set! obj 1 newval)))
(define set-particle-2!
(lambda (obj newval)
(if (particle? obj)
(void)
(error
'set-particle-2!
"~s is not a ~s"
obj
'<particle>))
(vector-set! obj 2 newval)))
(define set-particle-3!
(lambda (obj newval)
(if (particle? obj)
(void)
(error
'set-particle-3!
"~s is not a ~s"
obj
'<particle>))
(vector-set! obj 3 newval)))
(define set-particle-4!
(lambda (obj newval)
(if (particle? obj)
(void)
(error
'set-particle-4!
"~s is not a ~s"
obj
'<particle>))
(vector-set! obj 4 newval)))
(define set-particle-5!
(lambda (obj newval)
(if (particle? obj)
(void)
(error
'set-particle-5!
"~s is not a ~s"
obj
'<particle>))
(vector-set! obj 5 newval)))
(define set-particle-6!
(lambda (obj newval)
(if (particle? obj)
(void)
(error
'set-particle-6!
"~s is not a ~s"
obj
'<particle>))
(vector-set! obj 6 newval))))
(define make-particle
(lambda (position
acceleration
d-acceleration
potential
d-potential
strength)
((lambda ()
(make-raw-particle
position
acceleration
d-acceleration
potential
d-potential
strength)))))
(define particle-position particle-1)
(define particle-acceleration particle-2)
(define particle-d-acceleration particle-3)
(define particle-potential particle-4)
(define particle-d-potential particle-5)
(define particle-strength particle-6)
(define set-particle-position! set-particle-1!)
(define set-particle-acceleration! set-particle-2!)
(define set-particle-d-acceleration! set-particle-3!)
(define set-particle-potential! set-particle-4!)
(define set-particle-d-potential! set-particle-5!)
(define set-particle-strength! set-particle-6!))
(begin (begin (define make-raw-pt
(lambda (pt-1 pt-2 pt-3)
(vector '<pt> pt-1 pt-2 pt-3)))
(define pt?
(lambda (obj)
(if (vector? obj)
(if (= (vector-length obj) 4)
(eq? (vector-ref obj 0) '<pt>)
#f)
#f)))
(define pt-1
(lambda (obj)
(if (pt? obj)
(void)
(error 'pt-1 "~s is not a ~s" obj '<pt>))
(vector-ref obj 1)))
(define pt-2
(lambda (obj)
(if (pt? obj)
(void)
(error 'pt-2 "~s is not a ~s" obj '<pt>))
(vector-ref obj 2)))
(define pt-3
(lambda (obj)
(if (pt? obj)
(void)
(error 'pt-3 "~s is not a ~s" obj '<pt>))
(vector-ref obj 3)))
(define set-pt-1!
(lambda (obj newval)
(if (pt? obj)
(void)
(error
'set-pt-1!
"~s is not a ~s"
obj
'<pt>))
(vector-set! obj 1 newval)))
(define set-pt-2!
(lambda (obj newval)
(if (pt? obj)
(void)
(error
'set-pt-2!
"~s is not a ~s"
obj
'<pt>))
(vector-set! obj 2 newval)))
(define set-pt-3!
(lambda (obj newval)
(if (pt? obj)
(void)
(error
'set-pt-3!
"~s is not a ~s"
obj
'<pt>))
(vector-set! obj 3 newval))))
(define make-pt
(lambda (x y z) ((lambda () (make-raw-pt x y z)))))
(define pt-x pt-1)
(define pt-y pt-2)
(define pt-z pt-3)
(define set-pt-x! set-pt-1!)
(define set-pt-y! set-pt-2!)
(define set-pt-z! set-pt-3!)))
;(define-structure (tree
; body
; low-left-front-vertex
; up-right-back-vertex))
;
;(define-structure (node
; center
; low-left-front-vertex
; up-right-back-vertex
; children
; particles
; multipole-expansion
; near-field
; interactive-field))
;
;(define (leaf-node? node)
; (null? (node-children node)))
;
;(define-structure (particle
; position
; acceleration
; d-acceleration
; potential
; d-potential
; strength))
;
;(define-structure (pt x y z))
;
(define (pt-r pt)
(sqrt (+ (* (pt-x pt) (pt-x pt))
(* (pt-y pt) (pt-y pt))
(* (pt-z pt) (pt-z pt)))))
(define (pt-theta pt)
(let ((x (pt-x pt))
(y (pt-y pt))
(z (pt-z pt)))
(atan0 (sqrt (+ (* x x) (* y y))) z)))
(define (pt-phi pt)
(let ((x (pt-x pt)) (y (pt-y pt)))
(atan0 y x)))
(define (pt+ pt1 pt2)
(make-pt (+ (pt-x pt1) (pt-x pt2))
(+ (pt-y pt1) (pt-y pt2))
(+ (pt-z pt1) (pt-z pt2))))
(define (sum-vectors vectors)
(make-pt (apply:+ (map pt-x vectors))
(apply:+ (map pt-y vectors))
(apply:+ (map pt-z vectors))))
(define (pt- pt1 pt2)
(make-pt (- (pt-x pt1) (pt-x pt2))
(- (pt-y pt1) (pt-y pt2))
(- (pt-z pt1) (pt-z pt2))))
(define (pt-average pt1 pt2)
(pt-scalar* .5 (pt+ pt1 pt2)))
(define (pt-scalar* scalar pt)
(make-pt (* scalar (pt-x pt))
(* scalar (pt-y pt))
(* scalar (pt-z pt))))
(define (within-box? pt pt1 pt2)
(and (<= (pt-x pt) (pt-x pt2)) (> (pt-x pt) (pt-x pt1))
(<= (pt-y pt) (pt-y pt2)) (> (pt-y pt) (pt-y pt1))
(<= (pt-z pt) (pt-z pt2)) (> (pt-z pt) (pt-z pt1))))
;;; ==========================================================
;;; Useful Things
;;; ==========================================================
(define (nfilter list predicate)
(let loop ((list list))
(cond ((null? list) '())
((predicate (car list)) (cons (car list) (loop (cdr list))))
(else (loop (cdr list))))))
;;; array in the shape of a pyramid with each
;;; element a function of the indices
(define (make-cartesian-expansion func)
(let ((expansion (make-vector precision 0)))
(let loop1 ((i 0))
(if (= i precision)
expansion
(let ((foo (make-vector (- precision i) 0)))
(vector-set! expansion i foo)
(let loop2 ((j 0))
(if (= j (- precision i))
(loop1 (+ 1 i))
(let ((bar (make-vector (- precision i j) 0)))
(vector-set! foo j bar)
(let loop3 ((k 0))
(if (= k (- precision i j))
(loop2 (+ 1 j))
(begin (vector-set! bar k (func i j k))
(loop3 (+ 1 k)))))))))))))
;;; array in the shape of a triangle with each
;;; element a function of the indices
(define (make-spherical-expansion func)
(let ((expansion (make-vector precision 0)))
(let loop1 ((l 0))
(if (= l precision)
expansion
(let ((foo (make-vector (+ 1 l) 0)))
(vector-set! expansion l foo)
(let loop2 ((m 0))
(if (= m (+ 1 l))
(loop1 (+ 1 l))
(begin (vector-set! foo m (func l m))
(loop2 (+ 1 m))))))))))
(define (spherical-ref expansion l m)
(let ((conj (lambda (z) (make-rectangular (real-part z) (- (imag-part z))))))
(if (negative? m)
(conj (array-ref expansion l (- m)))
(array-ref expansion l m))))
(define (cartesian-expansion-sum expansions)
(make-cartesian-expansion
(lambda (i j k)
(apply:+ (map (lambda (expansion)
(array-ref expansion i j k))
expansions)))))
(define (spherical-expansion-sum expansions)
(make-spherical-expansion
(lambda (l m)
(apply:+ (map (lambda (expansion)
(spherical-ref expansion l m))
expansions)))))
(define (cartesian-zero-expansion)
(make-cartesian-expansion (lambda (i j k) 0)))
(define (spherical-zero-expansion)
(make-spherical-expansion (lambda (l m) 0)))
(define (sum-3d end1 end2 end3 func)
(let loop1 ((l 0) (sum 0))
(if (> l end1)
sum
(loop1
(+ 1 l)
(+ sum
(let loop2 ((m 0) (sum 0))
(if (> m end2)
sum
(loop2
(+ 1 m)
(+ sum
(let loop3 ((n 0) (sum 0))
(if (> n end3)
sum
(loop3 (+ 1 n)
(+ sum (func l m n))))))))))))))
(define (sum2-3d end func)
(let loop1 ((l 0) (sum 0))
(if (= l end)
sum
(loop1
(+ 1 l)
(+ sum
(let loop2 ((m 0) (sum 0))
(if (= (+ l m) end)
sum
(loop2
(+ 1 m)
(+ sum
(let loop3 ((n 0) (sum 0))
(if (= (+ l m n) end)
sum
(loop3 (+ 1 n)
(+ sum (func l m n))))))))))))))
(define (sum-2d end func)
(let loop1 ((l 0) (sum 0))
(if (> l end)
sum
(loop1 (+ 1 l)
(+ sum (let loop2 ((m (- l)) (sum 0))
(if (> m l)
sum
(loop2 (+ 1 m)
(+ sum (func l m))))))))))
(define (sum2-2d init func)
(let loop1 ((l init) (sum 0))
(if (= l precision)
sum
(loop1 (+ 1 l)
(+ sum (let loop2 ((m (- l)) (sum 0))
(if (> m l)
sum
(loop2 (+ 1 m)
(+ sum (func l m))))))))))
;;; Storing factorials in a table
(define fac
(let ((table (make-vector (* 4 precision) 0)))
(vector-set! table 0 1)
(let loop ((n 1))
(if (= n (* 4 precision))
(lambda (x) (vector-ref table x))
(begin (vector-set! table
n
(* n (vector-ref table (- n 1))))
(loop (+ 1 n)))))))
;;; The table for (* (-0.5) (-1.5) (-2.5) ... (+ -0.5 -i 1))
(define fac-1
(let ((table (make-vector precision 0)))
(vector-set! table 0 1)
(let loop ((n 1))
(if (= n precision)
(lambda (x) (vector-ref table x))
(begin (vector-set! table
n
(* (- .5 n)
(vector-ref table (- n 1))))
(loop (+ 1 n)))))))
(define fac-2
(let ((table (make-vector (* 4 precision) 0)))
(vector-set! table 0 1)
(let loop ((n 1))
(if (= n (* 4 precision))
(lambda (n) (if (< n 0)
1
(vector-ref table n)))
(begin (vector-set! table n (* (if (even? n) 1 n)
(vector-ref table (- n 1))))
(loop (+ 1 n)))))))
;;; Storing the products of factorials in a table.
(define prod-fac
(let ((table (make-cartesian-expansion
(lambda (i j k) (* (fac i) (fac j) (fac k))))))
(lambda (i j k) (array-ref table i j k))))
(define n1/prod-fac
(let ((table (make-cartesian-expansion
(lambda (i j k) (/ (prod-fac i j k))))))
(lambda (i j k) (array-ref table i j k))))
(define (assoc-legendre l m x)
(cond ((= l m) (* (expt -1 m)
(fac-2 (- (* 2 m) 1))
(expt (- 1 (* x x)) (/ m 2))))
((= l (+ 1 m)) (* x (+ 1 (* 2 m)) (assoc-legendre m m x)))
(else (/ (- (* x (- (* 2 l) 1) (assoc-legendre (- l 1) m x))
(* (+ l m -1) (assoc-legendre (- l 2) m x)))
(- l m)))))
(define (eval-spherical-harmonic l m theta phi)
(let ((mm (abs m)))
(* (sqrt (/ (fac (- l mm)) (fac (+ l mm))))
(assoc-legendre l mm (cos theta))
(make-polar 1 (* m phi)))))
(define (eval-spher-harm-phi-deriv l m theta phi)
(* (eval-spherical-harmonic l m theta phi)
m
(make-rectangular 0 1)))
(define (eval-spher-harm-theta-deriv l m theta phi)
(let ((mm (abs m)))
(* (sqrt (/ (fac (- l mm)) (fac (+ l mm))))
(make-polar 1 (* m phi))
(- (sin theta))
(assoc-legendre-deriv l mm (cos theta)))))
(define (assoc-legendre-deriv l m x)
(cond ((= l m) (* (expt -1 (+ 1 m))
(fac-2 (- (* 2 m) 1))
m
(expt (- 1 (* x x)) (- (/ m 2) 1))
x))
((= l (+ 1 m)) (* (+ 1 (* 2 m))
(+ (assoc-legendre m m x)
(* x (assoc-legendre-deriv m m x)))))
(else (/ (- (* (- (* 2 l) 1)
(+ (assoc-legendre (- l 1) m x)
(* x (assoc-legendre-deriv (- l 1) m x))))
(* (+ l m -1) (assoc-legendre-deriv (- l 2) m x)))
(- l m)))))
;;; ================================================================
;;; TREE CODE
;;; ================================================================
(define (build-tree height near-size)
(let* ((vertex1 (make-pt -10 -10 -10))
(vertex2 (make-pt 10 10 10))
(tree (make-tree '() vertex1 vertex2)))
(let loop ((level 0) (pt1 vertex1) (pt2 vertex2))
(let* ((half-diagonal (pt-scalar* .5 (pt- pt2 pt1)))
(diag-length/2 (pt-x half-diagonal)))
(insert-node tree level pt1 pt2)
(if (< level height)
(let ((child-pt1s
(map (lambda (offset)
(pt+ pt1 (pt-scalar* diag-length/2 offset)))
(list (make-pt 0 0 0)
(make-pt 0 0 1)
(make-pt 0 1 0)
(make-pt 1 0 0)
(make-pt 0 1 1)
(make-pt 1 0 1)
(make-pt 1 1 0)
(make-pt 1 1 1)))))
(for-each (lambda (child-pt1)
(loop (+ 1 level)
child-pt1
(pt+ child-pt1 half-diagonal)))
child-pt1s)))))
(calc-near-and-interaction tree near-size)
tree))
(define (insert-node tree level pt1 pt2)
(let* ((center (pt-average pt1 pt2))
(new-node (make-node center pt1 pt2 '() '() '() '() '())))
(letrec ((insert-internal (lambda (node depth)
(if (= level depth)
(set-node-children! node
(cons new-node
(node-children node)))
(insert-internal (find-child node center)
(+ 1 depth))))))
(if (= level 0)
(set-tree-body! tree new-node)
(insert-internal (tree-body tree) 1)))))
(define (find-child node pos)
(let loop ((children (node-children node)))
(let ((child (car children)))
(if (within-box? pos
(node-low-left-front-vertex child)
(node-up-right-back-vertex child))
child
(loop (cdr children))))))
(define (insert-particle tree particle)
(let* ((pos (particle-position particle)))
(letrec ((insert-internal (lambda (node)
(if (leaf-node? node)
(set-node-particles!
node
(cons particle (node-particles node)))
(insert-internal (find-child node pos))))))
(if (within-box? pos
(tree-low-left-front-vertex tree)
(tree-up-right-back-vertex tree))
(insert-internal (tree-body tree))
(error 'insert-particle "particle not within boundaries of tree" particle)))))
;;; This function finds the near and
;;; interaction fields for every node in the tree.
(define (calc-near-and-interaction tree near-size)
(set-node-near-field! (tree-body tree) (list (tree-body tree)))
(let loop ((node (tree-body tree)) (parent #f))
(if parent
(let* ((center (node-center node))
(dist (* near-size
(abs (- (pt-x center)
(pt-x (node-center parent)))))))
(for-each
(lambda (parent-near)
(let ((interactives (list '())))
(for-each
(lambda (child)
(if (> (pt-r (pt- center (node-center child))) dist)
(set-car! interactives (cons child (car interactives)))
(set-node-near-field!
node
(cons child (node-near-field node)))))
(node-children parent-near))
(set-node-interactive-field!
node
(append (car interactives) (node-interactive-field node)))))
(node-near-field parent))))
(for-each (lambda (child) (loop child node)) (node-children node))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; GO
(define (initial-particle x y z m)
(make-particle (make-pt x y z)
(make-pt 0 0 0)
(make-pt 0 0 0)
0 0 m))
(define (random-float bot top)
(+ (* (- top bot)
(/ (* (random 1000000) 1.0) 1000000.0))
bot))
(define (random-particle)
(make-particle (make-pt (random-float -10.0 10.0)
(random-float -10.0 10.0)
(random-float -10.0 10.0))
(make-pt 0 0 0)
(make-pt 0 0 0)
0 0 1.0))
(define *particles* (list '()))
(define (go depth precision n-particles)
(let ((tree (build-tree depth 27))
(particles (let next ((i 0) (ps '()))
(if (<= i n-particles)
(next (+ i 1) (cons (random-particle) ps))
ps))))
(for-each (lambda (p) (insert-particle tree p)) particles)
(cartesian-algorithm tree)
(set-car! *particles* particles)))
;;; virtual time for cartesian-algorithm step
;;; (go 1 3 10) 0.31 seconds
;;; (go 3 5 128) 1397.31
;;; (go 3 5 256) 1625.29
;;; (go 3 5 512) 2380.35
;;; (go 2 5 128) 27.44 seconds
(define (main . args)
(run-benchmark
"nbody"
nbody-iters
(lambda () #t)
(lambda (i j k) (go i j k))
2 5 128)))