117 lines
4.2 KiB
Racket
117 lines
4.2 KiB
Racket
#lang racket
|
|
|
|
(require math/flonum)
|
|
|
|
(provide random-exponential random-laplace
|
|
random-flonum random-single-flonum
|
|
random-integer->random-exact-rational
|
|
random-integer->random-flonum
|
|
random-integer->random-single-flonum
|
|
random-integer->random-real)
|
|
|
|
;; Returns a rational flonum from an exponential distribution
|
|
;; About half the numbers from this distribution are < scale, and half are > scale
|
|
(define (random-exponential [scale 1.0])
|
|
(let loop ()
|
|
(define r (random))
|
|
;; appx. probability 1.11e-16 this loops:
|
|
(cond [(= r 0) (loop)]
|
|
[else (exact->inexact (* (- (log (random))) scale))])))
|
|
|
|
;; Randomly signs a flonum from random-exponential
|
|
(define (random-laplace [scale 1.0])
|
|
(define r (random-exponential scale))
|
|
(if ((random) . < . 0.5) r (- r)))
|
|
|
|
;; Returns an integer from Uniform(0,2^64-1)
|
|
(define (random-flonum-bits)
|
|
(for/fold ([i 0]) ([_ (in-range 4)])
|
|
(+ (* i #x10000) (random #x10000))))
|
|
|
|
;; Returns a rational flonum with bits uniformly distributed
|
|
(define (random-flonum)
|
|
(let loop ()
|
|
(define x (bit-field->flonum (random-flonum-bits)))
|
|
;; appx. 0.0005 probability this loops
|
|
(if (rational? x) x (loop))))
|
|
|
|
;; Returns a single flonum with bits uniformly distributed; sometimes -inf.f or +inf.f
|
|
(define (random-single-flonum*)
|
|
(define s (random #b10)) ; sign bit
|
|
(define e (random #b100000000)) ; exponent
|
|
(define f (random #b100000000000000000000000)) ; fractional part
|
|
;; Formula from Wikipedia page on IEEE 754 binary32 format
|
|
(real->single-flonum
|
|
(* (if (= s 0) 1 -1)
|
|
(+ 1 (/ f #b100000000000000000000000))
|
|
(expt 2 (- e 127)))))
|
|
|
|
;; Returns a rational single flonum with bits uniformly distributed
|
|
(define (random-single-flonum)
|
|
(let loop ()
|
|
(define x (random-single-flonum*))
|
|
(if (rational? x) x (loop))))
|
|
|
|
;; Converts random integers to random exact rationals
|
|
(define (random-integer->random-exact-rational E)
|
|
(cond
|
|
[(= E 0) 0] ; code below would pick 0/0
|
|
[else
|
|
;; random fraction
|
|
(define n (exact-ceiling (random-exponential E)))
|
|
(define d (exact-ceiling (random-exponential E)))
|
|
(cond
|
|
[(= d 0) 0] ; appx. probability 1.11e-16
|
|
[else
|
|
(define x (/ n d))
|
|
(if ((random) . < . 0.5) x (- x))])]))
|
|
|
|
;; Converts random integers to random double flonums
|
|
;; Often returns very large or small numbers (near the limits of floating-point range)
|
|
;; Sometimes returns +nan.0, +inf.0 or -inf.0
|
|
(define (random-integer->random-flonum E)
|
|
(define r (random))
|
|
(cond
|
|
;; probability 0.25: random flonum, laplace-distributed with scale E
|
|
[(r . < . 0.25) (random-laplace (abs E))]
|
|
;; probability 0.25: random flonum with uniform bits
|
|
[(r . < . 0.50) (random-flonum)]
|
|
;; probability 0.35: very small or very large flonum
|
|
[(r . < . 0.85)
|
|
(define r (random))
|
|
(cond [(r . < . 0.5)
|
|
(define x (ordinal->flonum E))
|
|
(cond [(= x 0.0) (if ((random) . < . 0.5) 0.0 -0.0)]
|
|
[else x])]
|
|
[(r . < . 0.75)
|
|
(flstep -inf.0 (abs E))]
|
|
[else
|
|
(flstep +inf.0 (- (abs E)))])]
|
|
;; probability 0.05 each: +nan.0, +inf.0, -inf.0
|
|
[(r . < . 0.90) +nan.0]
|
|
[(r . < . 0.95) +inf.0]
|
|
[else -inf.0]))
|
|
|
|
;; Converts random integers to random single flonums
|
|
;; Sometimes returns +nan.f, +inf.f or -inf.f
|
|
(define (random-integer->random-single-flonum E)
|
|
(define r (random))
|
|
(cond
|
|
;; probability 0.50: random single flonum, laplace-distributed with scale E
|
|
[(r . < . 0.50) (real->single-flonum (random-laplace (abs E)))]
|
|
;; probability 0.35: random single flonum with uniform bits
|
|
[(r . < . 0.85) (random-single-flonum)]
|
|
;; probability 0.05 each: +nan.f, +inf.f, -inf.f
|
|
[(r . < . 0.90) +nan.f]
|
|
[(r . < . 0.95) +inf.f]
|
|
[else -inf.f]))
|
|
|
|
;; Converts random integers to random reals
|
|
(define (random-integer->random-real E)
|
|
(define r (random))
|
|
;; probability 0.25 each
|
|
(cond [(r . < . 0.25) E]
|
|
[(r . < . 0.50) (random-integer->random-exact-rational E)]
|
|
[(r . < . 0.75) (random-integer->random-flonum E)]
|
|
[else (random-integer->random-single-flonum E)]))
|