about summary refs log tree commit diff
path: root/third_party/lisp/alexandria/numbers.lisp
blob: 1c06f71d508f32c7e17b90621d2f686b72d37c30 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
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))