numbers.lisp - clic - Clic is an command line interactive client for gopher written in Common LISP (HTM) git clone git://bitreich.org/clic/ git://enlrupgkhuxnvlhsf6lc3fziv5h2hhfrinws65d7roiv6bfj7d652fid.onion/clic/ (DIR) Log (DIR) Files (DIR) Refs (DIR) Tags (DIR) README (DIR) 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 lesser then 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 secondary value, 15 optionally constrained by MIN and MAX. Gaussian random numbers form a standard 16 normal distribution around 0.0d0. 17 18 Sufficiently positive MIN or negative MAX will cause the algorithm used to 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 contagion 54 from STEP applied), each consequtive number being the sum of the previous one 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 contagion 72 from STEP applied), each consequtive number being the sum of the previous one 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 (the default), 155 and the unbiased estimator of variance if BIASED is false. SAMPLE must be a 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 estimator 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 the 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 the 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 bignums. 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 is the 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))