<div dir="ltr">FWIW, I ended up with this function:<div><br></div><div><div>(define (pick-a-maze maze-count)</div><div>  (+ (if (zero? (random 2))</div><div>         (/ maze-count 2)</div><div>         0)</div><div>     (random-natural (/ maze-count 4))</div>
<div>     (random-natural (/ maze-count 4))))</div></div><div><br></div><div style>which seems to be working out great.</div><div style><br></div><div style>Robby</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">
On Thu, Apr 11, 2013 at 12:53 PM, Joe Gilray <span dir="ltr">&lt;<a href="mailto:jgilray@gmail.com" target="_blank">jgilray@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">Really fun stuff.  Thanks Neil.<div><br></div><div>-Joe</div></div><div class="gmail_extra"><br><br><div class="gmail_quote"><div><div class="h5">On Mon, Apr 8, 2013 at 2:10 PM, Robby Findler <span dir="ltr">&lt;<a href="mailto:robby@eecs.northwestern.edu" target="_blank">robby@eecs.northwestern.edu</a>&gt;</span> wrote:<br>

</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><div dir="ltr">Thanks! What a great message!<span><font color="#888888"><div><br></div>
<div>Robby</div></font></span></div>
<div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Apr 8, 2013 at 2:47 PM, Neil Toronto <span dir="ltr">&lt;<a href="mailto:neil.toronto@gmail.com" target="_blank">neil.toronto@gmail.com</a>&gt;</span> wrote:<br>


<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><br>
<br>
On 04/02/2013 11:08 AM, Robby Findler wrote:<br>
</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<br>
than 604707155844263340879799987806<u></u>365649019991479003765974057968<u></u>835437533310870017962569786982<u></u>4,<br>
or other large integers?<br>
<br>
I wouldn&#39;t mind uniformly random (but if I get to make wishes, I&#39;d like<br>
to avoid small numbers, large numbers, and numbers near n/2).<br>
<br>
I tried looking in the finite distribution and integer distribution<br>
sections of the math library but didn&#39;t see anything.<br>
</blockquote>
<br></div>
It&#39;s not in `math/distributions&#39; because that module doesn&#39;t have integer distributions yet. (Well, it does, in the style of R, where the &quot;integers&quot; are actually flonums. They&#39;re fast, but probably not exact enough, or may not return random numbers big enough, for what you want them for.)<br>



<br>
The random big integer functions are sort of hidden in `math/base&#39;, as `random-natural&#39;, `random-integer&#39; and `random-bits&#39;. The last is the fastest; use it if you want a natural number less than some power of 2.<br>



<br>
The distributions for those are uniform, so to get your wishes, you&#39;ll have to be a bit tricky. Here&#39;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:<br>



<br>
#lang racket<br>
<br>
(require plot)<br>
<br>
(define is (build-list 20000 (λ (_) (random 100))))<br>
(define js (build-list 20000 (λ (_) (random 100))))<br>
(plot (density (map + is js)))<br>
<br>
<br>
Here&#39;s a more general function that chooses large natural numbers &lt; n by adding some number of uniformly random naturals:<br>
<br>
(require math)<br>
<br>
(define random-natural-parts 3)<br>
<br>
(define (my-random-natural n)<br>
  (define p (min n random-natural-parts))<br>
  (define m (quotient n p))<br>
  (define m-last (- n (* m (- p 1))))<br>
  (+ (apply + (build-list (- p 1) (λ (_) (random-natural (+ m 1)))))<br>
     (random-natural m-last)))<br>
<br>
(define xs<br>
  (build-list 20000 (λ (_) (my-random-natural<br>
                            151223462345987123498759238861<u></u>3))))<br>
(plot (density xs))<br>
(printf &quot;max xs = ~v~n&quot; (apply max xs))<br>
(printf &quot;min xs = ~v~n&quot; (apply min xs))<br>
<br>
<br>
It should work for all n &gt;= 1. Use `random-natural-parts&#39; 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.<br>



<br>
--<br>
<br>
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:<br>



<br>
(define (my-random-natural* n)<br>
  (define x (my-random-natural n))<br>
  (cond [(or (not (power-of-two? x)) ((random) . &lt; . 1/4))  x]<br>
        [else  (my-random-natural* n)]))<br>
<br>
;; Plot a histogram<br>
(define-values (ys cs)<br>
  (let*-values<br>
      ([(ys)     (build-list 20000 (λ (_) (my-random-natural* 40)))]<br>
       [(ys cs)  (count-samples ys)])<br>
    (sort-samples &lt; ys cs)))<br>
<br>
(plot (discrete-histogram (map list ys cs)))<br>
<br>
<br>
Here, &quot;rejection&quot; means trying again, since you always want an answer.<br>
<br>
To avoid numbers *near* a power of two, make the acceptance probability a function of the distance to the nearest power of two. We&#39;ll need a &quot;nearest power of two&quot; function first:<br>
<br>
(define (floor-log2 x)<br>
  (- (integer-length x) 1))<br>
<br>
(define (ceiling-log2 x)<br>
  (integer-length (- x 1)))<br>
<br>
(define (nearest-pow2 x)<br>
  (cond [(zero? x)  1]<br>
        [else<br>
         (define x0 (expt 2 (floor-log2 x)))<br>
         (define x1 (expt 2 (ceiling-log2 x)))<br>
         (if ((abs (- x x0)) . &lt; . (abs (- x x1))) x0 x1)]))<br>
<br>
<br>
Then make and use a function that returns acceptance probabilities:<br>
<br>
(define pow2-prob 1/4)<br>
(define pow2-spread 4)<br>
<br>
(define (accept-prob x)<br>
  (define pow2 (nearest-pow2 x))<br>
  (define d (min 1 (/ (abs (- x pow2)) pow2-spread)))<br>
  ;; A convex combination of 1 and pow2-prob, with fraction d<br>
  (+ d (* (- 1 d) pow2-prob)))<br>
<br>
(define (my-random-natural* n)<br>
  (define x (my-random-natural n))<br>
  (cond [((random) . &lt; . (accept-prob x))  x]<br>
        [else  (my-random-natural* n)]))<br>
<br>
<br>
Here, `pow2-prob&#39; is the probability a power of two is accepted, and `pow2-spread&#39; is the minimum distance an x can be from a power of two to always be accepted.<br>
<br>
&gt; (map accept-prob &#39;(16 17 18 19 20 21))<br>
&#39;(1/4 7/16 5/8 13/16 1 1)<br>
<br>
The `accept-prob&#39; function can be anything, as long as it returns probabilities, and doesn&#39;t return small probabilities so often that `my-random-natural*&#39; spends a long time looping.<span><font color="#888888"><br>



<br>
Neil ⊥<br>
<br>
</font></span></blockquote></div><br></div>
</div></div><br></div></div>____________________<br>
  Racket Users list:<br>
  <a href="http://lists.racket-lang.org/users" target="_blank">http://lists.racket-lang.org/users</a><br>
<br></blockquote></div><br></div>
</blockquote></div><br></div>