Merge commit '95aeb2ebae32a01ff79644daa523bda5d8552863' as 'third_party/lisp/alexandria' r/394

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 000000000000..1c06f71d508f
--- /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))