sbcl stuff
This commit is contained in:
Ian Keane
parent
1d1dbc34df
commit
5d91dbb667
335 changed files with 119806 additions and 1 deletions
|
|
@ -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))
|
||||
Loading…
Add table
Add a link
Reference in a new issue