numbers.lisp - clic - Clic is an command line interactive client for gopher wri… | |
git clone git://bitreich.org/clic/ git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws65… | |
Log | |
Files | |
Refs | |
Tags | |
README | |
LICENSE | |
--- | |
numbers.lisp (11218B) | |
--- | |
1 (in-package :alexandria) | |
2 | |
3 (declaim (inline clamp)) | |
4 (defun clamp (number min max) | |
5 "Clamps the NUMBER into [min, max] range. Returns MIN if NUMBER is les… | |
6 MIN and MAX if NUMBER is greater then MAX, otherwise returns NUMBER." | |
7 (if (< number min) | |
8 min | |
9 (if (> number max) | |
10 max | |
11 number))) | |
12 | |
13 (defun gaussian-random (&optional min max) | |
14 "Returns two gaussian random double floats as the primary and secondar… | |
15 optionally constrained by MIN and MAX. Gaussian random numbers form a st… | |
16 normal distribution around 0.0d0. | |
17 | |
18 Sufficiently positive MIN or negative MAX will cause the algorithm used … | |
19 take a very long time. If MIN is positive it should be close to zero, and | |
20 similarly if MAX is negative it should be close to zero." | |
21 (macrolet | |
22 ((valid (x) | |
23 `(<= (or min ,x) ,x (or max ,x)) )) | |
24 (labels | |
25 ((gauss () | |
26 (loop | |
27 for x1 = (- (random 2.0d0) 1.0d0) | |
28 for x2 = (- (random 2.0d0) 1.0d0) | |
29 for w = (+ (expt x1 2) (expt x2 2)) | |
30 when (< w 1.0d0) | |
31 do (let ((v (sqrt (/ (* -2.0d0 (log w)) w)))) | |
32 (return (values (* x1 v) (* x2 v)))))) | |
33 (guard (x) | |
34 (unless (valid x) | |
35 (tagbody | |
36 :retry | |
37 (multiple-value-bind (x1 x2) (gauss) | |
38 (when (valid x1) | |
39 (setf x x1) | |
40 (go :done)) | |
41 (when (valid x2) | |
42 (setf x x2) | |
43 (go :done)) | |
44 (go :retry)) | |
45 :done)) | |
46 x)) | |
47 (multiple-value-bind | |
48 (g1 g2) (gauss) | |
49 (values (guard g1) (guard g2)))))) | |
50 | |
51 (declaim (inline iota)) | |
52 (defun iota (n &key (start 0) (step 1)) | |
53 "Return a list of n numbers, starting from START (with numeric contagi… | |
54 from STEP applied), each consequtive number being the sum of the previou… | |
55 and STEP. START defaults to 0 and STEP to 1. | |
56 | |
57 Examples: | |
58 | |
59 (iota 4) => (0 1 2 3) | |
60 (iota 3 :start 1 :step 1.0) => (1.0 2.0 3.0) | |
61 (iota 3 :start -1 :step -1/2) => (-1 -3/2 -2) | |
62 " | |
63 (declare (type (integer 0) n) (number start step)) | |
64 (loop ;; KLUDGE: get numeric contagion right for the first element too | |
65 for i = (+ (- (+ start step) step)) then (+ i step) | |
66 repeat n | |
67 collect i)) | |
68 | |
69 (declaim (inline map-iota)) | |
70 (defun map-iota (function n &key (start 0) (step 1)) | |
71 "Calls FUNCTION with N numbers, starting from START (with numeric cont… | |
72 from STEP applied), each consequtive number being the sum of the previou… | |
73 and STEP. START defaults to 0 and STEP to 1. Returns N. | |
74 | |
75 Examples: | |
76 | |
77 (map-iota #'print 3 :start 1 :step 1.0) => 3 | |
78 ;;; 1.0 | |
79 ;;; 2.0 | |
80 ;;; 3.0 | |
81 " | |
82 (declare (type (integer 0) n) (number start step)) | |
83 (loop ;; KLUDGE: get numeric contagion right for the first element too | |
84 for i = (+ start (- step step)) then (+ i step) | |
85 repeat n | |
86 do (funcall function i)) | |
87 n) | |
88 | |
89 (declaim (inline lerp)) | |
90 (defun lerp (v a b) | |
91 "Returns the result of linear interpolation between A and B, using the | |
92 interpolation coefficient V." | |
93 ;; The correct version is numerically stable, at the expense of an | |
94 ;; extra multiply. See (lerp 0.1 4 25) with (+ a (* v (- b a))). The | |
95 ;; unstable version can often be converted to a fast instruction on | |
96 ;; a lot of machines, though this is machine/implementation | |
97 ;; specific. As alexandria is more about correct code, than | |
98 ;; efficiency, and we're only talking about a single extra multiply, | |
99 ;; many would prefer the stable version | |
100 (+ (* (- 1.0 v) a) (* v b))) | |
101 | |
102 (declaim (inline mean)) | |
103 (defun mean (sample) | |
104 "Returns the mean of SAMPLE. SAMPLE must be a sequence of numbers." | |
105 (/ (reduce #'+ sample) (length sample))) | |
106 | |
107 (defun median (sample) | |
108 "Returns median of SAMPLE. SAMPLE must be a sequence of real numbers." | |
109 ;; Implements and uses the quick-select algorithm to find the median | |
110 ;; https://en.wikipedia.org/wiki/Quickselect | |
111 | |
112 (labels ((randint-in-range (start-int end-int) | |
113 "Returns a random integer in the specified range, inclusive" | |
114 (+ start-int (random (1+ (- end-int start-int))))) | |
115 (partition (vec start-i end-i) | |
116 "Implements the partition function, which performs a partial | |
117 sort of vec around the (randomly) chosen pivot. | |
118 Returns the index where the pivot element would be located | |
119 in a correctly-sorted array" | |
120 (if (= start-i end-i) | |
121 start-i | |
122 (let ((pivot-i (randint-in-range start-i end-i))) | |
123 (rotatef (aref vec start-i) (aref vec pivot-i)) | |
124 (let ((swap-i end-i)) | |
125 (loop for i from swap-i downto (1+ start-i) do | |
126 (when (>= (aref vec i) (aref vec start-i)) | |
127 (rotatef (aref vec i) (aref vec swap-i)) | |
128 (decf swap-i))) | |
129 (rotatef (aref vec swap-i) (aref vec start-i)) | |
130 swap-i))))) | |
131 | |
132 (let* ((vector (copy-sequence 'vector sample)) | |
133 (len (length vector)) | |
134 (mid-i (ash len -1)) | |
135 (i 0) | |
136 (j (1- len))) | |
137 | |
138 (loop for correct-pos = (partition vector i j) | |
139 while (/= correct-pos mid-i) do | |
140 (if (< correct-pos mid-i) | |
141 (setf i (1+ correct-pos)) | |
142 (setf j (1- correct-pos)))) | |
143 | |
144 (if (oddp len) | |
145 (aref vector mid-i) | |
146 (* 1/2 | |
147 (+ (aref vector mid-i) | |
148 (reduce #'max (make-array | |
149 mid-i | |
150 :displaced-to vector)))))))) | |
151 | |
152 (declaim (inline variance)) | |
153 (defun variance (sample &key (biased t)) | |
154 "Variance of SAMPLE. Returns the biased variance if BIASED is true (th… | |
155 and the unbiased estimator of variance if BIASED is false. SAMPLE must b… | |
156 sequence of numbers." | |
157 (let ((mean (mean sample))) | |
158 (/ (reduce (lambda (a b) | |
159 (+ a (expt (- b mean) 2))) | |
160 sample | |
161 :initial-value 0) | |
162 (- (length sample) (if biased 0 1))))) | |
163 | |
164 (declaim (inline standard-deviation)) | |
165 (defun standard-deviation (sample &key (biased t)) | |
166 "Standard deviation of SAMPLE. Returns the biased standard deviation if | |
167 BIASED is true (the default), and the square root of the unbiased estima… | |
168 for variance if BIASED is false (which is not the same as the unbiased | |
169 estimator for standard deviation). SAMPLE must be a sequence of numbers." | |
170 (sqrt (variance sample :biased biased))) | |
171 | |
172 (define-modify-macro maxf (&rest numbers) max | |
173 "Modify-macro for MAX. Sets place designated by the first argument to … | |
174 maximum of its original value and NUMBERS.") | |
175 | |
176 (define-modify-macro minf (&rest numbers) min | |
177 "Modify-macro for MIN. Sets place designated by the first argument to … | |
178 minimum of its original value and NUMBERS.") | |
179 | |
180 ;;;; Factorial | |
181 | |
182 ;;; KLUDGE: This is really dependant on the numbers in question: for | |
183 ;;; small numbers this is larger, and vice versa. Ideally instead of a | |
184 ;;; constant we would have RANGE-FAST-TO-MULTIPLY-DIRECTLY-P. | |
185 (defconstant +factorial-bisection-range-limit+ 8) | |
186 | |
187 ;;; KLUDGE: This is really platform dependant: ideally we would use | |
188 ;;; (load-time-value (find-good-direct-multiplication-limit)) instead. | |
189 (defconstant +factorial-direct-multiplication-limit+ 13) | |
190 | |
191 (defun %multiply-range (i j) | |
192 ;; We use a a bit of cleverness here: | |
193 ;; | |
194 ;; 1. For large factorials we bisect in order to avoid expensive bignum | |
195 ;; multiplications: 1 x 2 x 3 x ... runs into bignums pretty soon, | |
196 ;; and once it does that all further multiplications will be with b… | |
197 ;; | |
198 ;; By instead doing the multiplication in a tree like | |
199 ;; ((1 x 2) x (3 x 4)) x ((5 x 6) x (7 x 8)) | |
200 ;; we manage to get less bignums. | |
201 ;; | |
202 ;; 2. Division isn't exactly free either, however, so we don't bisect | |
203 ;; all the way down, but multiply ranges of integers close to each | |
204 ;; other directly. | |
205 ;; | |
206 ;; For even better results it should be possible to use prime | |
207 ;; factorization magic, but Nikodemus ran out of steam. | |
208 ;; | |
209 ;; KLUDGE: We support factorials of bignums, but it seems quite | |
210 ;; unlikely anyone would ever be able to use them on a modern lisp, | |
211 ;; since the resulting numbers are unlikely to fit in memory... but | |
212 ;; it would be extremely unelegant to define FACTORIAL only on | |
213 ;; fixnums, _and_ on lisps with 16 bit fixnums this can actually be | |
214 ;; needed. | |
215 (labels ((bisect (j k) | |
216 (declare (type (integer 1 #.most-positive-fixnum) j k)) | |
217 (if (< (- k j) +factorial-bisection-range-limit+) | |
218 (multiply-range j k) | |
219 (let ((middle (+ j (truncate (- k j) 2)))) | |
220 (* (bisect j middle) | |
221 (bisect (+ middle 1) k))))) | |
222 (bisect-big (j k) | |
223 (declare (type (integer 1) j k)) | |
224 (if (= j k) | |
225 j | |
226 (let ((middle (+ j (truncate (- k j) 2)))) | |
227 (* (if (<= middle most-positive-fixnum) | |
228 (bisect j middle) | |
229 (bisect-big j middle)) | |
230 (bisect-big (+ middle 1) k))))) | |
231 (multiply-range (j k) | |
232 (declare (type (integer 1 #.most-positive-fixnum) j k)) | |
233 (do ((f k (* f m)) | |
234 (m (1- k) (1- m))) | |
235 ((< m j) f) | |
236 (declare (type (integer 0 (#.most-positive-fixnum)) m) | |
237 (type unsigned-byte f))))) | |
238 (if (and (typep i 'fixnum) (typep j 'fixnum)) | |
239 (bisect i j) | |
240 (bisect-big i j)))) | |
241 | |
242 (declaim (inline factorial)) | |
243 (defun %factorial (n) | |
244 (if (< n 2) | |
245 1 | |
246 (%multiply-range 1 n))) | |
247 | |
248 (defun factorial (n) | |
249 "Factorial of non-negative integer N." | |
250 (check-type n (integer 0)) | |
251 (%factorial n)) | |
252 | |
253 ;;;; Combinatorics | |
254 | |
255 (defun binomial-coefficient (n k) | |
256 "Binomial coefficient of N and K, also expressed as N choose K. This i… | |
257 number of K element combinations given N choises. N must be equal to or | |
258 greater then K." | |
259 (check-type n (integer 0)) | |
260 (check-type k (integer 0)) | |
261 (assert (>= n k)) | |
262 (if (or (zerop k) (= n k)) | |
263 1 | |
264 (let ((n-k (- n k))) | |
265 ;; Swaps K and N-K if K < N-K because the algorithm | |
266 ;; below is faster for bigger K and smaller N-K | |
267 (when (< k n-k) | |
268 (rotatef k n-k)) | |
269 (if (= 1 n-k) | |
270 n | |
271 ;; General case, avoid computing the 1x...xK twice: | |
272 ;; | |
273 ;; N! 1x...xN (K+1)x...xN | |
274 ;; -------- = ---------------- = ------------, N>1 | |
275 ;; K!(N-K)! 1x...xK x (N-K)! (N-K)! | |
276 (/ (%multiply-range (+ k 1) n) | |
277 (%factorial n-k)))))) | |
278 | |
279 (defun subfactorial (n) | |
280 "Subfactorial of the non-negative integer N." | |
281 (check-type n (integer 0)) | |
282 (if (zerop n) | |
283 1 | |
284 (do ((x 1 (1+ x)) | |
285 (a 0 (* x (+ a b))) | |
286 (b 1 a)) | |
287 ((= n x) a)))) | |
288 | |
289 (defun count-permutations (n &optional (k n)) | |
290 "Number of K element permutations for a sequence of N objects. | |
291 K defaults to N" | |
292 (check-type n (integer 0)) | |
293 (check-type k (integer 0)) | |
294 (assert (>= n k)) | |
295 (%multiply-range (1+ (- n k)) n)) |