about summary refs log tree commit diff
path: root/third_party/lisp/alexandria/numbers.lisp
diff options
context:
space:
mode:
Diffstat (limited to 'third_party/lisp/alexandria/numbers.lisp')
-rw-r--r--third_party/lisp/alexandria/numbers.lisp295
1 files changed, 0 insertions, 295 deletions
diff --git a/third_party/lisp/alexandria/numbers.lisp b/third_party/lisp/alexandria/numbers.lisp
deleted file mode 100644
index 1c06f71d508f..000000000000
--- a/third_party/lisp/alexandria/numbers.lisp
+++ /dev/null
@@ -1,295 +0,0 @@
-(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))