/lib/math.arc

  1(def zeros (x)
  2  (n-of x 0))
  3
  4(def ones (x)
  5  (n-of x 1))
  6
  7
  8;calculus fns
  9
 10(def deriv (f)
 11  "provides differential of a function of a single vairable"
 12  (fn (x)
 13    (let dx (max 1d-9 (abs:* x 1d-9))
 14       (/ (- (f (+ x dx))
 15	     (f x))
 16          dx))))
 17
 18(mac partial-diff (f n arity)
 19  "returns deriv function of f (w/ num args ARITY) wrt Nth variable(from 0) "
 20  (withs (arglis (n-of arity (uniq))
 21	  a-lis (firstn n arglis)
 22	  x arglis.n
 23	  b-lis (nthcdr (+ n 1) arglis))
 24    `(fn ,arglis
 25      (let dx (max 1d-9 (abs:* ,x 1d-9)
 26	(/ (- (,f ,@a-lis (+ ,x dx) ,@b-lis)
 27	      (,f ,@a-lis    ,x     ,@b-lis))
 28	   dx)))))
 29
 30(mac partial-diff-vec (f n arity)
 31  "returns vector with each element differentiated with respect to Nth argument"
 32  (withs (arglis (n-of arity (uniq))
 33	  a-lis (firstn n arglis)
 34	  x arglis.n
 35	  b-lis (nthcdr (+ n 1) arglis))
 36    `(fn ,arglis
 37	(let dx (max 1d-9 (abs:* ,x 1d-9)
 38	   (vec-scale (vec- (,f ,@a-lis	(+ ,x dx) ,@b-lis)
 39			    (,f ,@a-lis    ,x     ,@b-lis))
 40		      (/ 1 dx))))))
 41
 42(def grad (f)
 43  "gradient of 3D scalar field given by F"
 44  (fn (x y z)
 45    (list ((partial-diff f 0 3) x y z)
 46	  ((partial-diff f 1 3) x y z)
 47	  ((partial-diff f 2 3) x y z))))
 48
 49(def div (f)
 50  "divergence of 3D vector field given by F"
 51  (fn (x y z)
 52      (+ (((partial-diff-vec f 0 3) x y z) 0)
 53	 (((partial-diff-vec f 1 3) x y z) 1)
 54	 (((partial-diff-vec f 2 3) x y z) 2))))
 55
 56(def curl (f)
 57  "curl of 3D vector field given by F"
 58  (with (d/dx (partial-diff-vec f 0 3)
 59	 d/dy (partial-diff-vec f 1 3)
 60	 d/dz (partial-diff-vec f 2 3))
 61  (fn (x y z)
 62    (list (- ((d/dy x y z) 2)
 63	     ((d/dz x y z) 1))
 64	  (- ((d/dz x y z) 0)
 65	     ((d/dx x y z) 2))
 66	  (- ((d/dx x y z) 1)
 67	     ((d/dy x y z) 0))))))
 68
 69(def integral (f)
 70  "returns the integral of a single argument function"
 71  (fn (lower upper (o its 10000))
 72    (withs (dx      (/ (- upper lower) its)
 73	    x       lower
 74	    current (f x)
 75	    next    (f (+ x dx))
 76	    accum   0) 
 77      (while (<= x (- upper dx))
 78	     (++ accum (* (/ (+ current next) 2) dx))
 79	     (++ x dx)
 80	     (= current next)
 81	     (= next (f:+ x dx)))
 82      accum)))
 83
 84
 85
 86; vector fns
 87
 88(def vec-dot (v1 v2)
 89  (apply + (map (fn (x y) (* x y)) v1 v2)))
 90
 91(def vec-cross (v1 v2 . args)
 92  (if (car args)
 93      (vec-cross (vec-cross v1 v2) (car args) (cdr args))
 94      (list (- (* v1.1 v2.2) (* v2.1 v1.2))
 95	    (- (* v1.2 v2.0) (* v2.2 v1.0))
 96	    (- (* v1.0 v2.1) (* v2.0 v1.1)))))
 97
 98(with (v+ (fn (v1 v2)
 99	     (map (fn (x y) (+ x y)) v1 v2))
100       v- (fn (v1 v2)
101	     (map (fn (x y) (- x y)) v1 v2)))
102  
103  (def vec+ (v1 . args)
104    (if no.args v1
105	(reduce v+ (cons v1 args))))
106  
107  (def vec- (v1 . args)
108    (if no.args (map [- _] v1)
109	(v- v1 (apply vec+ args))))
110
111)
112
113
114(def vec-scale (vec . scalars)
115  (let c (apply * scalars)
116    (map [* _ c] vec)))
117
118(def quad-add args
119  ((afn (tot xs)
120     (if no.xs (sqrt tot)
121	 (self (+ tot (expt (car xs) 2)) (cdr xs))))
122   0 args))
123    
124(def vec-norm (vec)
125  (vec-scale vec (/ (apply quad-add vec))))
126
127
128;others
129
130(def fact (num)
131  "factorial of num (num must be a fixnum)"
132  (if (or (< num 0)(no:isa num 'int)) (err "num must be a positive integer")
133      ((afn (n x)
134	 (if (<= n 1)
135	     x
136	     (self (- n 1) (* n x))))
137       num 1)))
138
139(def n-bessel (n (o terms 100))
140  "gives a fn for the nth bessel function of the first kind evaluated at x"
141  (fn (x)
142    (with (i 0 
143	   tot 0)
144      (while (< i terms)
145	     (++ tot (/ (* (expt -1 i) (expt (/ x 2) (+ n i i)))
146			(* (fact i) (fact (+ n i)))))
147	     (++ i 1))
148      tot)))
149
150(def mean lis
151  ((afn (tot i x . xs)
152     (if no.xs (/ (+ tot x) (+ i 1))
153	 (self (+ tot x) (+ i 1) car.xs cdr.xs)))
154   0 0 car.lis cdr.lis))
155
156(def std-dev lis
157  (let m mean.lis
158    (sqrt:/ (apply + (map [expt (- _ m) 2] lis))
159	    len.lis)))
160
161(def choose (n x)
162  "number of ways of picking n elements from x (ie no. of ways of mixing 2 different sets of identical objects of size n and (- x n))"
163  (if (> n x)
164      (do (prn "")(pr n)(pr " is greater than ")(prn x))
165      (/ (fact x) (* (fact n) (fact (- x n))))))
166
167(def gauss-random (sigma (o middle 0))
168  "aproximation to a gausian distributed random with width sigma around middle, max possible deviation +/- 1000sigma"
169  (let tot 0
170    (for i 0 999
171      (++ tot  (- (* (rand) 2.0) 1)))
172    (+ (* tot (/ sigma 25)) middle)))
173
174(def quad-roots (a b c)
175  "returns roots of the equation ax²+bx+c=0"
176  (let sqroot (sqrt (- (* b b) (* 4 a c)))
177    (rem-dups
178     (list (/ (- sqroot b) 2 a)
179	   (/ (- 0 sqroot b) 2 a)))))