diff options
Diffstat (limited to 'third_party/lisp/alexandria/numbers.lisp')
| -rw-r--r-- | third_party/lisp/alexandria/numbers.lisp | 295 |
1 files changed, 295 insertions, 0 deletions
diff --git a/third_party/lisp/alexandria/numbers.lisp b/third_party/lisp/alexandria/numbers.lisp new file mode 100644 index 0000000000..1c06f71d50 --- /dev/null +++ b/third_party/lisp/alexandria/numbers.lisp @@ -0,0 +1,295 @@ +(in-package :alexandria) + +(declaim (inline clamp)) +(defun clamp (number min max) + "Clamps the NUMBER into [min, max] range. Returns MIN if NUMBER is lesser then +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 +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." + (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)))))) + +(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 +and STEP. START defaults to 0 and STEP to 1. + +Examples: + + (iota 4) => (0 1 2 3) + (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 ;; KLUDGE: get numeric contagion right for the first element too + for i = (+ (- (+ start step) step)) then (+ i step) + repeat n + collect i)) + +(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 +and STEP. START defaults to 0 and STEP to 1. Returns N. + +Examples: + + (map-iota #'print 3 :start 1 :step 1.0) => 3 + ;;; 1.0 + ;;; 2.0 + ;;; 3.0 +" + (declare (type (integer 0) n) (number start step)) + (loop ;; KLUDGE: get numeric contagion right for the first element too + for i = (+ start (- step step)) then (+ i step) + repeat n + do (funcall function i)) + n) + +(declaim (inline lerp)) +(defun lerp (v a b) + "Returns the result of linear interpolation between A and B, using the +interpolation coefficient V." + ;; The correct version is numerically stable, at the expense of an + ;; extra multiply. See (lerp 0.1 4 25) with (+ a (* v (- b a))). The + ;; unstable version can often be converted to a fast instruction on + ;; a lot of machines, though this is machine/implementation + ;; specific. As alexandria is more about correct code, than + ;; efficiency, and we're only talking about a single extra multiply, + ;; many would prefer the stable version + (+ (* (- 1.0 v) a) (* v b))) + +(declaim (inline mean)) +(defun mean (sample) + "Returns the mean of SAMPLE. SAMPLE must be a sequence of numbers." + (/ (reduce #'+ sample) (length sample))) + +(defun median (sample) + "Returns median of SAMPLE. SAMPLE must be a sequence of real numbers." + ;; Implements and uses the quick-select algorithm to find the median + ;; https://en.wikipedia.org/wiki/Quickselect + + (labels ((randint-in-range (start-int end-int) + "Returns a random integer in the specified range, inclusive" + (+ start-int (random (1+ (- end-int start-int))))) + (partition (vec start-i end-i) + "Implements the partition function, which performs a partial + sort of vec around the (randomly) chosen pivot. + Returns the index where the pivot element would be located + in a correctly-sorted array" + (if (= start-i end-i) + start-i + (let ((pivot-i (randint-in-range start-i end-i))) + (rotatef (aref vec start-i) (aref vec pivot-i)) + (let ((swap-i end-i)) + (loop for i from swap-i downto (1+ start-i) do + (when (>= (aref vec i) (aref vec start-i)) + (rotatef (aref vec i) (aref vec swap-i)) + (decf swap-i))) + (rotatef (aref vec swap-i) (aref vec start-i)) + swap-i))))) + + (let* ((vector (copy-sequence 'vector sample)) + (len (length vector)) + (mid-i (ash len -1)) + (i 0) + (j (1- len))) + + (loop for correct-pos = (partition vector i j) + while (/= correct-pos mid-i) do + (if (< correct-pos mid-i) + (setf i (1+ correct-pos)) + (setf j (1- correct-pos)))) + + (if (oddp len) + (aref vector mid-i) + (* 1/2 + (+ (aref vector mid-i) + (reduce #'max (make-array + mid-i + :displaced-to vector)))))))) + +(declaim (inline variance)) +(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." + (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)) +(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." + (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.") + +;;;; 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 + (let ((middle (+ j (truncate (- k j) 2)))) + (* (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) + (type unsigned-byte f))))) + (if (and (typep i 'fixnum) (typep j 'fixnum)) + (bisect i j) + (bisect-big i j)))) + +(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))) + ;; 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)) + (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)) + (if (zerop n) + 1 + (do ((x 1 (1+ x)) + (a 0 (* x (+ a b))) + (b 1 a)) + ((= n x) a)))) + +(defun count-permutations (n &optional (k n)) + "Number of K element permutations for a sequence of N objects. +K defaults to N" + (check-type n (integer 0)) + (check-type k (integer 0)) + (assert (>= n k)) + (%multiply-range (1+ (- n k)) n)) |