<div dir="ltr">Really fun stuff. Thanks Neil.<div> </div><div style>-Joe</div></div><div class="gmail_extra"> <div class="gmail_quote">On Mon, Apr 8, 2013 at 2:10 PM, Robby Findler <<a href="mailto:robby@eecs.northwestern.edu" target="_blank">robby@eecs.northwestern.edu</a>> wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Thanks! What a great message!<div> </div><div>Robby</div></div> <div class="HOEnZb"><div class="h5"><div class="gmail_extra"> <div class="gmail_quote">On Mon, Apr 8, 2013 at 2:47 PM, Neil Toronto <<a href="mailto:neil.toronto@gmail.com" target="_blank">neil.toronto@gmail.com</a>> wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Moving this to the list because it seems general-interest-y.<div> On 04/02/2013 11:08 AM, Robby Findler wrote: </div><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> Hi Neil; any advice on how to pick, at random, an natural number less than 6047071558442633408797999878063656490199914790037659740579688354375333108700179625697869824, or other large integers? I wouldn't mind uniformly random (but if I get to make wishes, I'd like to avoid small numbers, large numbers, and numbers near n/2). I tried looking in the finite distribution and integer distribution sections of the math library but didn't see anything. </blockquote> </div> It's not in `math/distributions' because that module doesn't have integer distributions yet. (Well, it does, in the style of R, where the "integers" are actually flonums. They're fast, but probably not exact enough, or may not return random numbers big enough, for what you want them for.) The random big integer functions are sort of hidden in `math/base', as `random-natural', `random-integer' and `random-bits'. The last is the fastest; use it if you want a natural number less than some power of 2. The distributions for those are uniform, so to get your wishes, you'll have to be a bit tricky. Here's one way: add them. If you add some number p of uniformly distributed numbers, as p increases, the distribution of the sum approaches a bell shape, which sounds like what you want. For example, adding two results in a triangular distribution: #lang racket (require plot) (define is (build-list 20000 (λ (_) (random 100)))) (define js (build-list 20000 (λ (_) (random 100)))) (plot (density (map + is js))) Here's a more general function that chooses large natural numbers < n by adding some number of uniformly random naturals: (require math) (define random-natural-parts 3) (define (my-random-natural n) (define p (min n random-natural-parts)) (define m (quotient n p)) (define m-last (- n (* m (- p 1)))) (+ (apply + (build-list (- p 1) (λ (_) (random-natural (+ m 1))))) (random-natural m-last))) (define xs (build-list 20000 (λ (_) (my-random-natural 1512234623459871234987592388613)))) (plot (density xs)) (printf "max xs = ~v~n" (apply max xs)) (printf "min xs = ~v~n" (apply min xs)) It should work for all n >= 1. Use `random-natural-parts' to control the shape: 1 for uniform, 2 for triangular, 3 for piecewise quadratic, etc., etc. The higher this number is, the less probable the numbers near 0 and n will be. -- To probabilistically avoid some kinds of numbers, you could use rejection sampling. If you wanted all powers of two to be 1/4 the probability they normally would be, you could accept them with probability 1/4 like this: (define (my-random-natural* n) (define x (my-random-natural n)) (cond [(or (not (power-of-two? x)) ((random) . < . 1/4)) x] [else (my-random-natural* n)])) ;; Plot a histogram (define-values (ys cs) (let*-values ([(ys) (build-list 20000 (λ (_) (my-random-natural* 40)))] [(ys cs) (count-samples ys)]) (sort-samples < ys cs))) (plot (discrete-histogram (map list ys cs))) Here, "rejection" means trying again, since you always want an answer. To avoid numbers *near* a power of two, make the acceptance probability a function of the distance to the nearest power of two. We'll need a "nearest power of two" function first: (define (floor-log2 x) (- (integer-length x) 1)) (define (ceiling-log2 x) (integer-length (- x 1))) (define (nearest-pow2 x) (cond [(zero? x) 1] [else (define x0 (expt 2 (floor-log2 x))) (define x1 (expt 2 (ceiling-log2 x))) (if ((abs (- x x0)) . < . (abs (- x x1))) x0 x1)])) Then make and use a function that returns acceptance probabilities: (define pow2-prob 1/4) (define pow2-spread 4) (define (accept-prob x) (define pow2 (nearest-pow2 x)) (define d (min 1 (/ (abs (- x pow2)) pow2-spread))) ;; A convex combination of 1 and pow2-prob, with fraction d (+ d (* (- 1 d) pow2-prob))) (define (my-random-natural* n) (define x (my-random-natural n)) (cond [((random) . < . (accept-prob x)) x] [else (my-random-natural* n)])) Here, `pow2-prob' is the probability a power of two is accepted, and `pow2-spread' is the minimum distance an x can be from a power of two to always be accepted. > (map accept-prob '(16 17 18 19 20 21)) '(1/4 7/16 5/8 13/16 1 1) The `accept-prob' function can be anything, as long as it returns probabilities, and doesn't return small probabilities so often that `my-random-natural*' spends a long time looping. Neil ⊥ </blockquote></div> </div> </div></div> ____________________ Racket Users list: <a href="http://lists.racket-lang.org/users" target="_blank">http://lists.racket-lang.org/users</a> </blockquote></div> </div>