GopherProxy

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