numbers.lisp

   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 repeat n
  65         ;; KLUDGE: get numeric contagion right for the first element too
  66         for i = (+ (- (+ start step) step)) then (+ i step)
  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 repeat n
  84         ;; KLUDGE: get numeric contagion right for the first element too
  85         for i = (+ start (- step step)) then (+ i step)
  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    (+ a (* v (- b a))))
  94
  95 (declaim (inline mean))
  96 (defun mean (sample)
  97   "Returns the mean of SAMPLE. SAMPLE must be a sequence of numbers."
  98   (/ (reduce #'+ sample) (length sample)))
  99
 100 (declaim (inline median))
 101 (defun median (sample)
 102   "Returns median of SAMPLE. SAMPLE must be a sequence of real numbers."
 103   (let* ((vector (sort (copy-sequence 'vector sample) #'<))
 104          (length (length vector))
 105          (middle (truncate length 2)))
 106     (if (oddp length)
 107         (aref vector middle)
 108         (/ (+ (aref vector middle) (aref vector (1- middle))) 2))))
 109
 110 (declaim (inline variance))
 111 (defun variance (sample &key (biased t))
 112   "Variance of SAMPLE. Returns the biased variance if BIASED is true (the default),
 113 and the unbiased estimator of variance if BIASED is false. SAMPLE must be a
 114 sequence of numbers."
 115   (let ((mean (mean sample)))
 116     (/ (reduce (lambda (a b)
 117                  (+ a (expt (- b mean) 2)))
 118                sample
 119                :initial-value 0)
 120        (- (length sample) (if biased 0 1)))))
 121
 122 (declaim (inline standard-deviation))
 123 (defun standard-deviation (sample &key (biased t))
 124   "Standard deviation of SAMPLE. Returns the biased standard deviation if
 125 BIASED is true (the default), and the square root of the unbiased estimator
 126 for variance if BIASED is false (which is not the same as the unbiased
 127 estimator for standard deviation). SAMPLE must be a sequence of numbers."
 128   (sqrt (variance sample :biased biased)))
 129
 130 (define-modify-macro maxf (&rest numbers) max
 131   "Modify-macro for MAX. Sets place designated by the first argument to the
 132 maximum of its original value and NUMBERS.")
 133
 134 (define-modify-macro minf (&rest numbers) min
 135   "Modify-macro for MIN. Sets place designated by the first argument to the
 136 minimum of its original value and NUMBERS.")
 137
 138 ;;;; Factorial
 139
 140 ;;; KLUDGE: This is really dependant on the numbers in question: for
 141 ;;; small numbers this is larger, and vice versa. Ideally instead of a
 142 ;;; constant we would have RANGE-FAST-TO-MULTIPLY-DIRECTLY-P.
 143 (defconstant +factorial-bisection-range-limit+ 8)
 144
 145 ;;; KLUDGE: This is really platform dependant: ideally we would use
 146 ;;; (load-time-value (find-good-direct-multiplication-limit)) instead.
 147 (defconstant +factorial-direct-multiplication-limit+ 13)
 148
 149 (defun %multiply-range (i j)
 150   ;; We use a a bit of cleverness here:
 151   ;;
 152   ;; 1. For large factorials we bisect in order to avoid expensive bignum
 153   ;;    multiplications: 1 x 2 x 3 x ... runs into bignums pretty soon,
 154   ;;    and once it does that all further multiplications will be with bignums.
 155   ;;
 156   ;;    By instead doing the multiplication in a tree like
 157   ;;       ((1 x 2) x (3 x 4)) x ((5 x 6) x (7 x 8))
 158   ;;    we manage to get less bignums.
 159   ;;
 160   ;; 2. Division isn't exactly free either, however, so we don't bisect
 161   ;;    all the way down, but multiply ranges of integers close to each
 162   ;;    other directly.
 163   ;;
 164   ;; For even better results it should be possible to use prime
 165   ;; factorization magic, but Nikodemus ran out of steam.
 166   ;;
 167   ;; KLUDGE: We support factorials of bignums, but it seems quite
 168   ;; unlikely anyone would ever be able to use them on a modern lisp,
 169   ;; since the resulting numbers are unlikely to fit in memory... but
 170   ;; it would be extremely unelegant to define FACTORIAL only on
 171   ;; fixnums, _and_ on lisps with 16 bit fixnums this can actually be
 172   ;; needed.
 173   (labels ((bisect (j k)
 174              (declare (type (integer 1 #.most-positive-fixnum) j k))
 175              (if (< (- k j) +factorial-bisection-range-limit+)
 176                  (multiply-range j k)
 177                  (let ((middle (+ j (truncate (- k j) 2))))
 178                    (* (bisect j middle)
 179                       (bisect (+ middle 1) k)))))
 180            (bisect-big (j k)
 181              (declare (type (integer 1) j k))
 182              (if (= j k)
 183                  j
 184                  (let ((middle (+ j (truncate (- k j) 2))))
 185                    (* (if (<= middle most-positive-fixnum)
 186                           (bisect j middle)
 187                           (bisect-big j middle))
 188                       (bisect-big (+ middle 1) k)))))
 189            (multiply-range (j k)
 190              (declare (type (integer 1 #.most-positive-fixnum) j k))
 191              (do ((f k (* f m))
 192                   (m (1- k) (1- m)))
 193                  ((< m j) f)
 194                (declare (type (integer 0 (#.most-positive-fixnum)) m)
 195                         (type unsigned-byte f)))))
 196     (if (and (typep i 'fixnum) (typep j 'fixnum))
 197         (bisect i j)
 198         (bisect-big i j))))
 199
 200 (declaim (inline factorial))
 201 (defun %factorial (n)
 202   (if (< n 2)
 203       1
 204       (%multiply-range 1 n)))
 205
 206 (defun factorial (n)
 207   "Factorial of non-negative integer N."
 208   (check-type n (integer 0))
 209   (%factorial n))
 210
 211 ;;;; Combinatorics
 212
 213 (defun binomial-coefficient (n k)
 214   "Binomial coefficient of N and K, also expressed as N choose K. This is the
 215 number of K element combinations given N choises. N must be equal to or
 216 greater then K."
 217   (check-type n (integer 0))
 218   (check-type k (integer 0))
 219   (assert (>= n k))
 220   (if (or (zerop k) (= n k))
 221       1
 222       (let ((n-k (- n k)))
 223         ;; Swaps K and N-K if K < N-K because the algorithm
 224         ;; below is faster for bigger K and smaller N-K
 225         (when (< k n-k)
 226           (rotatef k n-k))
 227         (if (= 1 n-k)
 228             n
 229             ;; General case, avoid computing the 1x...xK twice:
 230             ;;
 231             ;;    N!           1x...xN          (K+1)x...xN
 232             ;; --------  =  ---------------- =  ------------, N>1
 233             ;; K!(N-K)!     1x...xK x (N-K)!       (N-K)!
 234             (/ (%multiply-range (+ k 1) n)
 235                (%factorial n-k))))))
 236
 237 (defun subfactorial (n)
 238   "Subfactorial of the non-negative integer N."
 239   (check-type n (integer 0))
 240   (if (zerop n)
 241       1
 242       (do ((x 1 (1+ x))
 243            (a 0 (* x (+ a b)))
 244            (b 1 a))
 245           ((= n x) a))))
 246
 247 (defun count-permutations (n &optional (k n))
 248   "Number of K element permutations for a sequence of N objects.
 249 K defaults to N"
 250   (check-type n (integer 0))
 251   (check-type k (integer 0))
 252   (assert (>= n k))
 253   (%multiply-range (1+ (- n k)) n))