numbers.lisp 9.46 KB
Newer Older
Nikodemus Siivola's avatar
Nikodemus Siivola committed
1 2 3 4
(in-package :alexandria)

(declaim (inline clamp))
(defun clamp (number min max)
5
  "Clamps the NUMBER into [min, max] range. Returns MIN if NUMBER is lesser then
Nikodemus Siivola's avatar
Nikodemus Siivola committed
6 7 8 9 10 11 12 13 14 15
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
16 17 18 19 20
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."
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
  (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))))))
Nikodemus Siivola's avatar
Nikodemus Siivola committed
50 51 52 53 54

(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
55
and STEP. START defaults to 0 and STEP to 1.
Nikodemus Siivola's avatar
Nikodemus Siivola committed
56 57

Examples:
Nikodemus Siivola's avatar
Nikodemus Siivola committed
58

59
  (iota 4)                      => (0 1 2 3)
Nikodemus Siivola's avatar
Nikodemus Siivola committed
60 61 62 63 64
  (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 repeat n
65 66 67
        ;; KLUDGE: get numeric contagion right for the first element too
        for i = (+ (- (+ start step) step)) then (+ i step)
        collect i))
Nikodemus Siivola's avatar
Nikodemus Siivola committed
68

69 70 71 72
(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
73
and STEP. START defaults to 0 and STEP to 1. Returns N.
74 75 76

Examples:

77 78 79 80
  (map-iota #'print 3 :start 1 :step 1.0) => 3
    ;;; 1.0
    ;;; 2.0
    ;;; 3.0
81 82 83 84 85 86 87 88
"
  (declare (type (integer 0) n) (number start step))
  (loop repeat n
        ;; KLUDGE: get numeric contagion right for the first element too
        for i = (+ start (- step step)) then (+ i step)
        do (funcall function i))
  n)

Nikodemus Siivola's avatar
Nikodemus Siivola committed
89 90 91 92
(declaim (inline lerp))
(defun lerp (v a b)
  "Returns the result of linear interpolation between A and B, using the
interpolation coefficient V."
93 94 95 96 97 98 99 100
  ;; 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)))
Nikodemus Siivola's avatar
Nikodemus Siivola committed
101

102 103 104 105
(declaim (inline mean))
(defun mean (sample)
  "Returns the mean of SAMPLE. SAMPLE must be a sequence of numbers."
  (/ (reduce #'+ sample) (length sample)))
Nikodemus Siivola's avatar
Nikodemus Siivola committed
106

107 108 109 110
(declaim (inline median))
(defun median (sample)
  "Returns median of SAMPLE. SAMPLE must be a sequence of real numbers."
  (let* ((vector (sort (copy-sequence 'vector sample) #'<))
Nikodemus Siivola's avatar
Nikodemus Siivola committed
111 112 113 114
         (length (length vector))
         (middle (truncate length 2)))
    (if (oddp length)
        (aref vector middle)
Johan Ur Riise's avatar
Johan Ur Riise committed
115
        (/ (+ (aref vector middle) (aref vector (1- middle))) 2))))
Nikodemus Siivola's avatar
Nikodemus Siivola committed
116 117

(declaim (inline variance))
118 119 120 121
(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."
Nikodemus Siivola's avatar
Nikodemus Siivola committed
122 123 124 125 126 127 128 129
  (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))
130 131 132 133 134
(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."
Nikodemus Siivola's avatar
Nikodemus Siivola committed
135 136 137 138 139 140 141 142 143
  (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.")
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190

;;;; 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
191
                 (let ((middle (+ j (truncate (- k j) 2))))
192 193 194 195 196 197 198 199 200 201
                   (* (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)
202
                        (type unsigned-byte f)))))
203 204 205
    (if (and (typep i 'fixnum) (typep j 'fixnum))
        (bisect i j)
        (bisect-big i j))))
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229

(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)))
230 231 232 233
        ;; 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))
234 235 236 237 238 239 240 241 242 243 244 245 246
        (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))
247 248 249 250 251 252
  (if (zerop n)
      1
      (do ((x 1 (1+ x))
           (a 0 (* x (+ a b)))
           (b 1 a))
          ((= n x) a))))
253 254 255

(defun count-permutations (n &optional (k n))
  "Number of K element permutations for a sequence of N objects.
256 257 258 259 260
K defaults to N"
  (check-type n (integer 0))
  (check-type k (integer 0))
  (assert (>= n k))
  (%multiply-range (1+ (- n k)) n))