(frob %random-single-float single-float)
(frob %random-double-float double-float))
-;;; Mersenne Twister RNG
-;;;
-;;; FIXME: It's unpleasant to have RANDOM functionality scattered
-;;; through the code this way. It would be nice to move this into the
-;;; same file as the other RANDOM definitions.
+;;; Return an expression to generate an integer of N-BITS many random
+;;; bits, using the minimal number of random chunks possible.
+(defun generate-random-expr-for-power-of-2 (n-bits state)
+ (declare (type (integer 1 #.sb!vm:n-word-bits) n-bits))
+ (multiple-value-bind (n-chunk-bits chunk-expr)
+ (cond ((<= n-bits n-random-chunk-bits)
+ (values n-random-chunk-bits `(random-chunk ,state)))
+ ((<= n-bits (* 2 n-random-chunk-bits))
+ (values (* 2 n-random-chunk-bits) `(big-random-chunk ,state)))
+ (t
+ (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))
+ (if (< n-bits n-chunk-bits)
+ `(logand ,(1- (ash 1 n-bits)) ,chunk-expr)
+ chunk-expr)))
+
+;;; This transform for compile-time constant word-sized integers
+;;; generates an accept-reject loop to achieve equidistribution of the
+;;; returned values. Several optimizations are done: If NUM is a power
+;;; of two no loop is needed. If the random chunk size is half the word
+;;; size only one chunk is used where sufficient. For values of NUM
+;;; where it is possible and results in faster code, the rejection
+;;; probability is reduced by accepting all values below the largest
+;;; multiple of the limit that fits into one or two chunks and and doing
+;;; a division to get the random value into the desired range.
(deftransform random ((num &optional state)
- ((integer 1 #.(expt 2 sb!vm::n-word-bits)) &optional *))
- ;; FIXME: I almost conditionalized this as #!+sb-doc. Find some way
- ;; of automatically finding #!+sb-doc in proximity to DEFTRANSFORM
- ;; to let me scan for places that I made this mistake and didn't
- ;; catch myself.
- "use inline (UNSIGNED-BYTE 32) operations"
- (let ((type (lvar-type num))
- (limit (expt 2 sb!vm::n-word-bits))
- (random-chunk (ecase sb!vm::n-word-bits
- (32 'random-chunk)
- (64 'sb!kernel::big-random-chunk))))
- (if (numeric-type-p type)
- (let ((num-high (numeric-type-high (lvar-type num))))
- (aver num-high)
- (cond ((constant-lvar-p num)
- ;; Check the worst case sum absolute error for the
- ;; random number expectations.
- (let ((rem (rem limit num-high)))
- (unless (< (/ (* 2 rem (- num-high rem))
- num-high limit)
- (expt 2 (- sb!kernel::random-integer-extra-bits)))
- (give-up-ir1-transform
- "The random number expectations are inaccurate."))
- (if (= num-high limit)
- `(,random-chunk (or state *random-state*))
- #!-(or x86 x86-64)
- `(rem (,random-chunk (or state *random-state*)) num)
- #!+(or x86 x86-64)
- ;; Use multiplication, which is faster.
- `(values (sb!bignum::%multiply
- (,random-chunk (or state *random-state*))
- num)))))
- ((> num-high random-fixnum-max)
- (give-up-ir1-transform
- "The range is too large to ensure an accurate result."))
- #!+(or x86 x86-64)
- ((< num-high limit)
- `(values (sb!bignum::%multiply
- (,random-chunk (or state *random-state*))
- num)))
- (t
- `(rem (,random-chunk (or state *random-state*)) num))))
- ;; KLUDGE: a relatively conservative treatment, but better
- ;; than a bug (reported by PFD sbcl-devel towards the end of
- ;; 2004-11.
- (give-up-ir1-transform
- "Argument type is too complex to optimize for."))))
+ ((constant-arg (integer 1 #.(expt 2 sb!vm:n-word-bits)))
+ &optional *)
+ *
+ :policy (and (> speed compilation-speed)
+ (> speed space)))
+ "optimize to inlined RANDOM-CHUNK operations"
+ (let ((num (lvar-value num)))
+ (if (= num 1)
+ 0
+ (flet ((chunk-n-bits-and-expr (n-bits)
+ (cond ((<= n-bits n-random-chunk-bits)
+ (values n-random-chunk-bits
+ '(random-chunk (or state *random-state*))))
+ ((<= n-bits (* 2 n-random-chunk-bits))
+ (values (* 2 n-random-chunk-bits)
+ '(big-random-chunk (or state *random-state*))))
+ (t
+ (error "Unexpectedly small N-RANDOM-CHUNK-BITS")))))
+ (if (zerop (logand num (1- num)))
+ ;; NUM is a power of 2.
+ (let ((n-bits (integer-length (1- num))))
+ (multiple-value-bind (n-chunk-bits chunk-expr)
+ (chunk-n-bits-and-expr n-bits)
+ (if (< n-bits n-chunk-bits)
+ `(logand ,(1- (ash 1 n-bits)) ,chunk-expr)
+ chunk-expr)))
+ ;; Generate an accept-reject loop.
+ (let ((n-bits (integer-length num)))
+ (multiple-value-bind (n-chunk-bits chunk-expr)
+ (chunk-n-bits-and-expr n-bits)
+ (if (or (> (* num 3) (expt 2 n-chunk-bits))
+ (logbitp (- n-bits 2) num))
+ ;; Division can't help as the quotient is below 3,
+ ;; or is too costly as the rejection probability
+ ;; without it is already small (namely at most 1/4
+ ;; with the given test, which is experimentally a
+ ;; reasonable threshold and cheap to test for).
+ `(loop
+ (let ((bits ,(generate-random-expr-for-power-of-2
+ n-bits '(or state *random-state*))))
+ (when (< bits num)
+ (return bits))))
+ (let ((d (truncate (expt 2 n-chunk-bits) num)))
+ `(loop
+ (let ((bits ,chunk-expr))
+ (when (< bits ,(* num d))
+ (return (values (truncate bits ,d)))))))))))))))
+
\f
;;;; float accessors
(if (< x ,most-negative)
,most-negative
(coerce x ',type)))
- (numeric-type-low num)))
+ (numeric-type-low num)
+ nil))
(hi (bound-func (lambda (x)
(if (< ,most-positive x )
,most-positive
(coerce x ',type)))
- (numeric-type-high num))))
+ (numeric-type-high num)
+ nil)))
(specifier-type `(,',type ,(or lo '*) ,(or hi '*)))))
(defoptimizer (,fun derive-type) ((num))
;; Process the intersection.
(let* ((low (interval-low intersection))
(high (interval-high intersection))
- (res-lo (or (bound-func fun (if increasingp low high))
+ (res-lo (or (bound-func fun (if increasingp low high) nil)
default-low))
- (res-hi (or (bound-func fun (if increasingp high low))
+ (res-hi (or (bound-func fun (if increasingp high low) nil)
default-high))
(format (case (numeric-type-class arg)
((integer rational) 'single-float)
(int-hi (if hi
(ceiling (type-bound-number hi))
'*))
- (f-lo (or (bound-func #'float lo)
+ (f-lo (or (bound-func #'float lo nil)
'*))
- (f-hi (or (bound-func #'float hi)
+ (f-hi (or (bound-func #'float hi nil)
'*)))
(specifier-type `(or (rational ,int-lo ,int-hi)
(single-float ,f-lo, f-hi)))))
(float
;; A positive integer to a float power is a float.
- (modified-numeric-type y-type
- :low (interval-low bnd)
- :high (interval-high bnd)))
+ (let ((format (numeric-type-format y-type)))
+ (aver format)
+ (modified-numeric-type
+ y-type
+ :low (coerce-numeric-bound (interval-low bnd) format)
+ :high (coerce-numeric-bound (interval-high bnd) format))))
(t
;; A positive integer to a number is a number (for now).
(specifier-type 'number))))
(int-hi (if hi
(ceiling (type-bound-number hi))
'*))
- (f-lo (or (bound-func #'float lo)
+ (f-lo (or (bound-func #'float lo nil)
'*))
- (f-hi (or (bound-func #'float hi)
+ (f-hi (or (bound-func #'float hi nil)
'*)))
(specifier-type `(or (rational ,int-lo ,int-hi)
(single-float ,f-lo, f-hi)))))
(float
;; A positive rational to a float power is a float.
- (modified-numeric-type y-type
- :low (interval-low bnd)
- :high (interval-high bnd)))
+ (let ((format (numeric-type-format y-type)))
+ (aver format)
+ (modified-numeric-type
+ y-type
+ :low (coerce-numeric-bound (interval-low bnd) format)
+ :high (coerce-numeric-bound (interval-high bnd) format))))
(t
;; A positive rational to a number is a number (for now).
(specifier-type 'number))))
((or integer rational)
;; A positive float to an integer or rational power is
;; always a float.
- (make-numeric-type
- :class 'float
- :format (numeric-type-format x-type)
- :low (interval-low bnd)
- :high (interval-high bnd)))
+ (let ((format (numeric-type-format x-type)))
+ (aver format)
+ (make-numeric-type
+ :class 'float
+ :format format
+ :low (coerce-numeric-bound (interval-low bnd) format)
+ :high (coerce-numeric-bound (interval-high bnd) format))))
(float
;; A positive float to a float power is a float of the
;; higher type.
- (make-numeric-type
- :class 'float
- :format (float-format-max (numeric-type-format x-type)
- (numeric-type-format y-type))
- :low (interval-low bnd)
- :high (interval-high bnd)))
+ (let ((format (float-format-max (numeric-type-format x-type)
+ (numeric-type-format y-type))))
+ (aver format)
+ (make-numeric-type
+ :class 'float
+ :format format
+ :low (coerce-numeric-bound (interval-low bnd) format)
+ :high (coerce-numeric-bound (interval-high bnd) format))))
(t
;; A positive float to a number is a number (for now)
(specifier-type 'number))))
:complexp :real
:low (numeric-type-low type)
:high (numeric-type-high type))))))
-#-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
+
(defoptimizer (realpart derive-type) ((num))
(one-arg-derive-type num #'realpart-derive-type-aux #'realpart))
+
(defun imagpart-derive-type-aux (type)
(let ((class (numeric-type-class type))
(format (numeric-type-format type)))
:complexp :real
:low (numeric-type-low type)
:high (numeric-type-high type))))))
-#-sb-xc-host ; (See CROSS-FLOAT-INFINITY-KLUDGE.)
+
(defoptimizer (imagpart derive-type) ((num))
(one-arg-derive-type num #'imagpart-derive-type-aux #'imagpart))
;; exactly the same way as the functions themselves do
;; it.
(if (csubtypep arg domain)
- (let ((res-lo (bound-func fun (numeric-type-low arg)))
- (res-hi (bound-func fun (numeric-type-high arg))))
+ (let ((res-lo (bound-func fun (numeric-type-low arg) nil))
+ (res-hi (bound-func fun (numeric-type-high arg) nil)))
(unless increasingp
(rotatef res-lo res-hi))
(make-numeric-type
projects / sbcl.git / blobdiff