main.js (6420B)
1 /** 2 * @license Apache-2.0 3 * 4 * Copyright (c) 2018 The Stdlib Authors. 5 * 6 * Licensed under the Apache License, Version 2.0 (the "License"); 7 * you may not use this file except in compliance with the License. 8 * You may obtain a copy of the License at 9 * 10 * http://www.apache.org/licenses/LICENSE-2.0 11 * 12 * Unless required by applicable law or agreed to in writing, software 13 * distributed under the License is distributed on an "AS IS" BASIS, 14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 * See the License for the specific language governing permissions and 16 * limitations under the License. 17 */ 18 19 'use strict'; 20 21 // MODULES // 22 23 var factory = require( './factory.js' ); 24 25 26 // MAIN // 27 28 /** 29 * Returns a pseudorandom number drawn from a discrete uniform distribution with minimum support `a` and maximum support `b`. 30 * 31 * ## Method 32 * 33 * - Let \\( R \\) be a pseudorandom generator (PRNG) which yields integers on the interval \\( \[ A, B ] \\). 34 * 35 * - If \\( a = b \\), then \\( rv = a \\). 36 * 37 * - Let \\( r1 = b - a \\) and \\( r2 = B - A \\). If \\( r2 = r1 \\) (likely a rare occurrence), then 38 * 39 * ```tex 40 * rv = ( R() - B ) + a 41 * ``` 42 * 43 * where, for real integer types, operation order is important in order to avoid overflow. 44 * 45 * - If \\( r2 < r1 \\), use rejection sampling to map random variates from \\( R \\) to a larger domain (e.g., \\( {0,1,2,3} \rightarrow {0,1,2,3,4} \\)). For real integer types (and floating-point integer values), we must take extra care to avoid overflow. During sampling, the following conditions will hold: 46 * 47 * - First, consider the post-condition: \\( \textrm{result} \leq r2 \\), thus avoiding overflow. 48 * 49 * - Begin with definition of \\( \textrm{limit} \\) 50 * 51 * ```tex 52 * \textrm{limit} = \lfloor{\frac{r2+1}{r1+1}\rfloor 53 * ``` 54 * 55 * thus, 56 * 57 * ```tex 58 * \textrm{limit}\ \cdot (r1+1) \leq r2+1 59 * ``` 60 * 61 * - Let \\( m \\) be a random factor where the loop condition is defined as 62 * 63 * ```tex 64 * m \leq \textrm{limit} 65 * ``` 66 * 67 * - Let \\( \textrm{result} \\) be the generator output, which is expressed base \\( r2+1 \\) and obeys the loop invariant \\( \textrm{result} < m \\). 68 * 69 * - Let \\( rv \\) be a realization of the PRNG. Then, 70 * 71 * ```tex 72 * rv-A \leq r1 73 * ``` 74 * 75 * and, by the loop condition, \\( m \leq \textrm{limit} \\). 76 * 77 * - Therefore, 78 * 79 * ```tex 80 * m \cdot (rv - A + 1) \leq r2+1 81 * ``` 82 * 83 * - Rearranging terms, 84 * 85 * ```tex 86 * m + m \cdot (rv - A) \leq r2+1 87 * ``` 88 * 89 * - Since \\( \textrm{result} < m \\), 90 * 91 * ```tex 92 * \textrm{result} + m \cdot (rv - A) < r2+1 93 * ``` 94 * 95 * - Next, consider the post-condition: \\( \textrm{result} < m \cdot (r2+1) \\). 96 * 97 * - Since \\( \textrm{result} < m \\) and \\( rv - A \leq r1 \\), 98 * 99 * ```tex 100 * \textrm{result} + m \cdot (rv - A) < m + m \cdot (rv - A) 101 * ``` 102 * 103 * - Therefore, 104 * 105 * ```tex 106 * \textrm{result} + m \cdot (rv - A) < m + m \cdot r1 107 * ``` 108 * 109 * - Therefore, 110 * 111 * ```tex 112 * \textrm{result} + m \cdot (rv - A) < m \cdot (r1+1) 113 * ``` 114 * 115 * - Next, consider the post-condition: \\( m \leq r2 \\). 116 * 117 * - According to the definition of \\( \textrm{limit} \\) and the loop condition \\( m \leq \textrm{limit} \\), 118 * 119 * ```tex 120 * m \cdot (r1+1) \leq r2+1 121 * ``` 122 * 123 * - If \\( r2 \\) is **not** an integer power of the generator range \\( r1 \\), i.e., 124 * 125 * ```tex 126 * m \cdot (r1+1) \neq r2+1 127 * ``` 128 * 129 * then 130 * 131 * ```tex 132 * m \cdot (r1+1) < r2+1 133 * ``` 134 * 135 * - Thus, \\( \textrm{result} < m \\). 136 * 137 * - Next, consider the post-condition: \\( r2/m < r1+1 \\). 138 * 139 * - To show this is true, let us try to prove its opposite. Given the loop condition \\( m > \textrm{limit} \\), assume 140 * 141 * ```tex 142 * r2/m > r1+1 143 * ``` 144 * 145 * - Accordingly, 146 * 147 * ```tex 148 * r2 \geq m \cdot (r1+1) 149 * ``` 150 * 151 * - Hence, 152 * 153 * ```tex 154 * r2+1 > m \cdot (r1+1) 155 * ``` 156 * 157 * - Using the loop condition, 158 * 159 * ```tex 160 * r2+1 > (\textrm{limit}+1) \cdot (r1+1) 161 * ``` 162 * 163 * - Rearranging terms, 164 * 165 * ```tex 166 * \frac{r2+1}{r1+1} > \textrm{limit} + 1 167 * ``` 168 * 169 * - Hence, 170 * 171 * ```tex 172 * \textrm{limit} < \lfloor{\frac{r2+1}{r1+1}} \rfloor 173 * ``` 174 * 175 * - But the definition of \\( \textrm{limit} \\) is 176 * 177 * ```tex 178 * \textrm{limit} = \lfloor{\frac{r2+1}{r1+1}} 179 * ``` 180 * 181 * - Thus, our assumption cannot be true, providing the post-condition by reductio ad absurdum. 182 * 183 * - Next, consider the post-condition 184 * 185 * ```tex 186 * r2 \leq \frac{r2}{m} \cdot m + (m - 1) 187 * ``` 188 * 189 * - Recall the identity 190 * 191 * ```tex 192 * r2 = \frac{r2}{m} \cdot m + r2 \mod m 193 * ``` 194 * 195 * - By the definition of the modulus 196 * 197 * ```tex 198 * r2 \mod m < m 199 * ``` 200 * 201 * - Therefore, 202 * 203 * ```tex 204 * r2 < \frac{r2}{m} \cdot m + m 205 * ``` 206 * 207 * - Hence, 208 * 209 * ```tex 210 * r2 \leq \frac{r2}{m} \cdot m + (m - 1) 211 * ``` 212 * 213 * At this point, the maximum value \\( \textrm{result} \\) is \\( m-1 \\). Hence, we can generate numbers that can be at least as large as \\( r2 \\), but we must be careful to avoid overflow during addition and in the sampling rejection. Anything which overflows is larger than \\( r2 \\) and can thus be rejected. 214 * 215 * - If \\( r1 > r2 \\), use rejection sampling to map random variates from \\( R \\) to a smaller domain (e.g., \\( {0,1,2,3,4} \rightarrow {0,1,2,3} \\)) by defining "buckets" in which multiple random variates in \\( R \\) map to a single random variate in the smaller domain. We are safe in adding 1 to \\( r2 \\); however, we need to be careful to not cause overflow when adding 1 to \\( r1 \\). 216 * 217 * @name discreteUniform 218 * @type {PRNG} 219 * @param {integer} a - minimum support 220 * @param {integer} b - maximum support 221 * @returns {integer} pseudorandom number 222 * 223 * @example 224 * var v = discreteUniform( 1, 10 ); 225 * // returns <number> 226 */ 227 var discreteUniform = factory(); 228 229 230 // EXPORTS // 231 232 module.exports = discreteUniform;