numbers.lisp 9.02 KB
 Nikodemus Siivola committed Oct 15, 2006 1 2 3 4 ``````(in-package :alexandria) (declaim (inline clamp)) (defun clamp (number min max) `````` Nikodemus Siivola committed Mar 09, 2010 5 `````` "Clamps the NUMBER into [min, max] range. Returns MIN if NUMBER is lesser then `````` Nikodemus Siivola committed Oct 15, 2006 6 7 8 9 10 11 12 13 14 15 ``````MIN and MAX if NUMBER is greater then MAX, otherwise returns NUMBER." (if (< number min) min (if (> number max) max number))) (defun gaussian-random (&optional min max) "Returns two gaussian random double floats as the primary and secondary value, optionally constrained by MIN and MAX. Gaussian random numbers form a standard `````` Nikodemus Siivola committed Jan 26, 2013 16 17 18 19 20 ``````normal distribution around 0.0d0. Sufficiently positive MIN or negative MAX will cause the algorithm used to take a very long time. If MIN is positive it should be close to zero, and similarly if MAX is negative it should be close to zero." `````` aarvid committed Feb 16, 2014 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 `````` (macrolet ((valid (x) `(<= (or min ,x) ,x (or max ,x)) )) (labels ((gauss () (loop for x1 = (- (random 2.0d0) 1.0d0) for x2 = (- (random 2.0d0) 1.0d0) for w = (+ (expt x1 2) (expt x2 2)) when (< w 1.0d0) do (let ((v (sqrt (/ (* -2.0d0 (log w)) w)))) (return (values (* x1 v) (* x2 v)))))) (guard (x) (unless (valid x) (tagbody :retry (multiple-value-bind (x1 x2) (gauss) (when (valid x1) (setf x x1) (go :done)) (when (valid x2) (setf x x2) (go :done)) (go :retry)) :done)) x)) (multiple-value-bind (g1 g2) (gauss) (values (guard g1) (guard g2)))))) `````` Nikodemus Siivola committed Oct 15, 2006 50 51 52 53 54 `````` (declaim (inline iota)) (defun iota (n &key (start 0) (step 1)) "Return a list of n numbers, starting from START (with numeric contagion from STEP applied), each consequtive number being the sum of the previous one `````` Attila Lendvai committed Mar 10, 2008 55 ``````and STEP. START defaults to 0 and STEP to 1. `````` Nikodemus Siivola committed Oct 15, 2006 56 57 `````` Examples: `````` Nikodemus Siivola committed Oct 16, 2006 58 `````` `````` Nikodemus Siivola committed Jan 26, 2013 59 `````` (iota 4) => (0 1 2 3) `````` Nikodemus Siivola committed Oct 15, 2006 60 61 62 63 64 `````` (iota 3 :start 1 :step 1.0) => (1.0 2.0 3.0) (iota 3 :start -1 :step -1/2) => (-1 -3/2 -2) " (declare (type (integer 0) n) (number start step)) (loop repeat n `````` Nikodemus Siivola committed Jan 26, 2013 65 66 67 `````` ;; KLUDGE: get numeric contagion right for the first element too for i = (+ (- (+ start step) step)) then (+ i step) collect i)) `````` Nikodemus Siivola committed Oct 15, 2006 68 `````` `````` Nikodemus Siivola committed Jun 01, 2007 69 70 71 72 ``````(declaim (inline map-iota)) (defun map-iota (function n &key (start 0) (step 1)) "Calls FUNCTION with N numbers, starting from START (with numeric contagion from STEP applied), each consequtive number being the sum of the previous one `````` Attila Lendvai committed Mar 10, 2008 73 ``````and STEP. START defaults to 0 and STEP to 1. Returns N. `````` Nikodemus Siivola committed Jun 01, 2007 74 75 76 `````` Examples: `````` Attila Lendvai committed Mar 10, 2008 77 78 79 80 `````` (map-iota #'print 3 :start 1 :step 1.0) => 3 ;;; 1.0 ;;; 2.0 ;;; 3.0 `````` Nikodemus Siivola committed Jun 01, 2007 81 82 83 84 85 86 87 88 ``````" (declare (type (integer 0) n) (number start step)) (loop repeat n ;; KLUDGE: get numeric contagion right for the first element too for i = (+ start (- step step)) then (+ i step) do (funcall function i)) n) `````` Nikodemus Siivola committed Oct 15, 2006 89 90 91 92 93 94 ``````(declaim (inline lerp)) (defun lerp (v a b) "Returns the result of linear interpolation between A and B, using the interpolation coefficient V." (+ a (* v (- b a)))) `````` Nikodemus Siivola committed Mar 30, 2012 95 96 97 98 ``````(declaim (inline mean)) (defun mean (sample) "Returns the mean of SAMPLE. SAMPLE must be a sequence of numbers." (/ (reduce #'+ sample) (length sample))) `````` Nikodemus Siivola committed Oct 15, 2006 99 `````` `````` Nikodemus Siivola committed Mar 30, 2012 100 101 102 103 ``````(declaim (inline median)) (defun median (sample) "Returns median of SAMPLE. SAMPLE must be a sequence of real numbers." (let* ((vector (sort (copy-sequence 'vector sample) #'<)) `````` Nikodemus Siivola committed Oct 15, 2006 104 105 106 107 `````` (length (length vector)) (middle (truncate length 2))) (if (oddp length) (aref vector middle) `````` Johan Ur Riise committed Sep 16, 2009 108 `````` (/ (+ (aref vector middle) (aref vector (1- middle))) 2)))) `````` Nikodemus Siivola committed Oct 15, 2006 109 110 `````` (declaim (inline variance)) `````` 111 112 113 114 ``````(defun variance (sample &key (biased t)) "Variance of SAMPLE. Returns the biased variance if BIASED is true (the default), and the unbiased estimator of variance if BIASED is false. SAMPLE must be a sequence of numbers." `````` Nikodemus Siivola committed Oct 15, 2006 115 116 117 118 119 120 121 122 `````` (let ((mean (mean sample))) (/ (reduce (lambda (a b) (+ a (expt (- b mean) 2))) sample :initial-value 0) (- (length sample) (if biased 0 1))))) (declaim (inline standard-deviation)) `````` 123 124 125 126 127 ``````(defun standard-deviation (sample &key (biased t)) "Standard deviation of SAMPLE. Returns the biased standard deviation if BIASED is true (the default), and the square root of the unbiased estimator for variance if BIASED is false (which is not the same as the unbiased estimator for standard deviation). SAMPLE must be a sequence of numbers." `````` Nikodemus Siivola committed Oct 15, 2006 128 129 130 131 132 133 134 135 136 `````` (sqrt (variance sample :biased biased))) (define-modify-macro maxf (&rest numbers) max "Modify-macro for MAX. Sets place designated by the first argument to the maximum of its original value and NUMBERS.") (define-modify-macro minf (&rest numbers) min "Modify-macro for MIN. Sets place designated by the first argument to the minimum of its original value and NUMBERS.") `````` Nikodemus Siivola committed Jul 01, 2007 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 `````` ;;;; Factorial ;;; KLUDGE: This is really dependant on the numbers in question: for ;;; small numbers this is larger, and vice versa. Ideally instead of a ;;; constant we would have RANGE-FAST-TO-MULTIPLY-DIRECTLY-P. (defconstant +factorial-bisection-range-limit+ 8) ;;; KLUDGE: This is really platform dependant: ideally we would use ;;; (load-time-value (find-good-direct-multiplication-limit)) instead. (defconstant +factorial-direct-multiplication-limit+ 13) (defun %multiply-range (i j) ;; We use a a bit of cleverness here: ;; ;; 1. For large factorials we bisect in order to avoid expensive bignum ;; multiplications: 1 x 2 x 3 x ... runs into bignums pretty soon, ;; and once it does that all further multiplications will be with bignums. ;; ;; By instead doing the multiplication in a tree like ;; ((1 x 2) x (3 x 4)) x ((5 x 6) x (7 x 8)) ;; we manage to get less bignums. ;; ;; 2. Division isn't exactly free either, however, so we don't bisect ;; all the way down, but multiply ranges of integers close to each ;; other directly. ;; ;; For even better results it should be possible to use prime ;; factorization magic, but Nikodemus ran out of steam. ;; ;; KLUDGE: We support factorials of bignums, but it seems quite ;; unlikely anyone would ever be able to use them on a modern lisp, ;; since the resulting numbers are unlikely to fit in memory... but ;; it would be extremely unelegant to define FACTORIAL only on ;; fixnums, _and_ on lisps with 16 bit fixnums this can actually be ;; needed. (labels ((bisect (j k) (declare (type (integer 1 #.most-positive-fixnum) j k)) (if (< (- k j) +factorial-bisection-range-limit+) (multiply-range j k) (let ((middle (+ j (truncate (- k j) 2)))) (* (bisect j middle) (bisect (+ middle 1) k))))) (bisect-big (j k) (declare (type (integer 1) j k)) (if (= j k) j `````` Luís Oliveira committed Jul 11, 2007 184 `````` (let ((middle (+ j (truncate (- k j) 2)))) `````` Nikodemus Siivola committed Jul 01, 2007 185 186 187 188 189 190 191 192 193 194 `````` (* (if (<= middle most-positive-fixnum) (bisect j middle) (bisect-big j middle)) (bisect-big (+ middle 1) k))))) (multiply-range (j k) (declare (type (integer 1 #.most-positive-fixnum) j k)) (do ((f k (* f m)) (m (1- k) (1- m))) ((< m j) f) (declare (type (integer 0 (#.most-positive-fixnum)) m) `````` Stelian Ionescu committed Aug 06, 2007 195 `````` (type unsigned-byte f))))) `````` Jan Moringen committed May 12, 2014 196 197 198 `````` (if (and (typep i 'fixnum) (typep j 'fixnum)) (bisect i j) (bisect-big i j)))) `````` Nikodemus Siivola committed Jul 01, 2007 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 `````` (declaim (inline factorial)) (defun %factorial (n) (if (< n 2) 1 (%multiply-range 1 n))) (defun factorial (n) "Factorial of non-negative integer N." (check-type n (integer 0)) (%factorial n)) ;;;; Combinatorics (defun binomial-coefficient (n k) "Binomial coefficient of N and K, also expressed as N choose K. This is the number of K element combinations given N choises. N must be equal to or greater then K." (check-type n (integer 0)) (check-type k (integer 0)) (assert (>= n k)) (if (or (zerop k) (= n k)) 1 (let ((n-k (- n k))) `````` Nikodemus Siivola committed Sep 23, 2010 223 224 225 226 `````` ;; Swaps K and N-K if K < N-K because the algorithm ;; below is faster for bigger K and smaller N-K (when (< k n-k) (rotatef k n-k)) `````` Nikodemus Siivola committed Jul 01, 2007 227 228 229 230 231 232 233 234 235 236 237 238 239 `````` (if (= 1 n-k) n ;; General case, avoid computing the 1x...xK twice: ;; ;; N! 1x...xN (K+1)x...xN ;; -------- = ---------------- = ------------, N>1 ;; K!(N-K)! 1x...xK x (N-K)! (N-K)! (/ (%multiply-range (+ k 1) n) (%factorial n-k)))))) (defun subfactorial (n) "Subfactorial of the non-negative integer N." (check-type n (integer 0)) `````` Mason Smith committed Oct 30, 2011 240 241 242 243 244 245 `````` (if (zerop n) 1 (do ((x 1 (1+ x)) (a 0 (* x (+ a b))) (b 1 a)) ((= n x) a)))) `````` Nikodemus Siivola committed Jul 01, 2007 246 247 248 `````` (defun count-permutations (n &optional (k n)) "Number of K element permutations for a sequence of N objects. `````` Nikodemus Siivola committed Sep 23, 2010 249 250 251 252 253 ``````K defaults to N" (check-type n (integer 0)) (check-type k (integer 0)) (assert (>= n k)) (%multiply-range (1+ (- n k)) n))``````