1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.rng.sampling.distribution; 19 20 import org.apache.commons.rng.UniformRandomProvider; 21 22 /** 23 * Discrete uniform distribution sampler. 24 * 25 * <p>Sampling uses {@link UniformRandomProvider#nextInt}.</p> 26 * 27 * <p>When the range is a power of two the number of calls is 1 per sample. 28 * Otherwise a rejection algorithm is used to ensure uniformity. In the worst 29 * case scenario where the range spans half the range of an {@code int} 30 * (2<sup>31</sup> + 1) the expected number of calls is 2 per sample.</p> 31 * 32 * <p>This sampler can be used as a replacement for {@link UniformRandomProvider#nextInt} 33 * with appropriate adjustment of the upper bound to be inclusive and will outperform that 34 * method when the range is not a power of two. The advantage is gained by pre-computation 35 * of the rejection threshold.</p> 36 * 37 * <p>The sampling algorithm is described in:</p> 38 * 39 * <blockquote> 40 * Lemire, D (2019). <i>Fast Random Integer Generation in an Interval.</i> 41 * <strong>ACM Transactions on Modeling and Computer Simulation</strong> 29 (1). 42 * </blockquote> 43 * 44 * <p>The number of {@code int} values required per sample follows a geometric distribution with 45 * a probability of success p of {@code 1 - ((2^32 % n) / 2^32)}. This requires on average 1/p random 46 * {@code int} values per sample.</p> 47 * 48 * @see <a href="https://arxiv.org/abs/1805.10941">Fast Random Integer Generation in an Interval</a> 49 * 50 * @since 1.0 51 */ 52 public class DiscreteUniformSampler 53 extends SamplerBase 54 implements SharedStateDiscreteSampler { 55 56 /** The appropriate uniform sampler for the parameters. */ 57 private final SharedStateDiscreteSampler delegate; 58 59 /** 60 * Base class for a sampler from a discrete uniform distribution. This contains the 61 * source of randomness. 62 */ 63 private abstract static class AbstractDiscreteUniformSampler 64 implements SharedStateDiscreteSampler { 65 /** Underlying source of randomness. */ 66 protected final UniformRandomProvider rng; 67 68 /** 69 * @param rng Generator of uniformly distributed random numbers. 70 */ 71 AbstractDiscreteUniformSampler(UniformRandomProvider rng) { 72 this.rng = rng; 73 } 74 75 @Override 76 public String toString() { 77 return "Uniform deviate [" + rng.toString() + "]"; 78 } 79 } 80 81 /** 82 * Discrete uniform distribution sampler when the sample value is fixed. 83 */ 84 private static class FixedDiscreteUniformSampler 85 extends AbstractDiscreteUniformSampler { 86 /** The value. */ 87 private final int value; 88 89 /** 90 * @param value The value. 91 */ 92 FixedDiscreteUniformSampler(int value) { 93 // No requirement for the RNG 94 super(null); 95 this.value = value; 96 } 97 98 @Override 99 public int sample() { 100 return value; 101 } 102 103 @Override 104 public String toString() { 105 // No RNG to include in the string 106 return "Uniform deviate [X=" + value + "]"; 107 } 108 109 @Override 110 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 111 // No requirement for the RNG 112 return this; 113 } 114 } 115 116 /** 117 * Discrete uniform distribution sampler when the range is a power of 2 and greater than 1. 118 * This sampler assumes the lower bound of the range is 0. 119 * 120 * <p>Note: This cannot be used when the range is 1 (2^0) as the shift would be 32-bits 121 * which is ignored by the shift operator.</p> 122 */ 123 private static class PowerOf2RangeDiscreteUniformSampler 124 extends AbstractDiscreteUniformSampler { 125 /** Bit shift to apply to the integer sample. */ 126 private final int shift; 127 128 /** 129 * @param rng Generator of uniformly distributed random numbers. 130 * @param range Maximum range of the sample (exclusive). 131 * Must be a power of 2 greater than 2^0. 132 */ 133 PowerOf2RangeDiscreteUniformSampler(UniformRandomProvider rng, 134 int range) { 135 super(rng); 136 this.shift = Integer.numberOfLeadingZeros(range) + 1; 137 } 138 139 /** 140 * @param rng Generator of uniformly distributed random numbers. 141 * @param source Source to copy. 142 */ 143 PowerOf2RangeDiscreteUniformSampler(UniformRandomProvider rng, 144 PowerOf2RangeDiscreteUniformSampler source) { 145 super(rng); 146 this.shift = source.shift; 147 } 148 149 @Override 150 public int sample() { 151 // Use a bit shift to favour the most significant bits. 152 // Note: The result would be the same as the rejection method used in the 153 // small range sampler when there is no rejection threshold. 154 return rng.nextInt() >>> shift; 155 } 156 157 @Override 158 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 159 return new PowerOf2RangeDiscreteUniformSampler(rng, this); 160 } 161 } 162 163 /** 164 * Discrete uniform distribution sampler when the range is small 165 * enough to fit in a positive integer. 166 * This sampler assumes the lower bound of the range is 0. 167 * 168 * <p>Implements the algorithm of Lemire (2019).</p> 169 * 170 * @see <a href="https://arxiv.org/abs/1805.10941">Fast Random Integer Generation in an Interval</a> 171 */ 172 private static class SmallRangeDiscreteUniformSampler 173 extends AbstractDiscreteUniformSampler { 174 /** Maximum range of the sample (exclusive). */ 175 private final long n; 176 177 /** 178 * The level below which samples are rejected based on the fraction remainder. 179 * 180 * <p>Any remainder below this denotes that there are still floor(2^32 / n) more 181 * observations of this sample from the interval [0, 2^32), where n is the range.</p> 182 */ 183 private final long threshold; 184 185 /** 186 * @param rng Generator of uniformly distributed random numbers. 187 * @param range Maximum range of the sample (exclusive). 188 */ 189 SmallRangeDiscreteUniformSampler(UniformRandomProvider rng, 190 int range) { 191 super(rng); 192 // Handle range as an unsigned 32-bit integer 193 this.n = range & 0xffffffffL; 194 // Compute 2^32 % n 195 threshold = (1L << 32) % n; 196 } 197 198 /** 199 * @param rng Generator of uniformly distributed random numbers. 200 * @param source Source to copy. 201 */ 202 SmallRangeDiscreteUniformSampler(UniformRandomProvider rng, 203 SmallRangeDiscreteUniformSampler source) { 204 super(rng); 205 this.n = source.n; 206 this.threshold = source.threshold; 207 } 208 209 @Override 210 public int sample() { 211 // Rejection method using multiply by a fraction: 212 // n * [0, 2^32 - 1) 213 // ------------- 214 // 2^32 215 // The result is mapped back to an integer and will be in the range [0, n). 216 // Note this is comparable to range * rng.nextDouble() but with compensation for 217 // non-uniformity due to round-off. 218 long result; 219 do { 220 // Compute 64-bit unsigned product of n * [0, 2^32 - 1). 221 // The upper 32-bits contains the sample value in the range [0, n), i.e. result / 2^32. 222 // The lower 32-bits contains the remainder, i.e. result % 2^32. 223 result = n * (rng.nextInt() & 0xffffffffL); 224 225 // Test the sample uniformity. 226 // Samples are observed on average (2^32 / n) times at a frequency of either 227 // floor(2^32 / n) or ceil(2^32 / n). 228 // To ensure all samples have a frequency of floor(2^32 / n) reject any results with 229 // a remainder < (2^32 % n), i.e. the level below which denotes that there are still 230 // floor(2^32 / n) more observations of this sample. 231 } while ((result & 0xffffffffL) < threshold); 232 // Divide by 2^32 to get the sample. 233 return (int)(result >>> 32); 234 } 235 236 @Override 237 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 238 return new SmallRangeDiscreteUniformSampler(rng, this); 239 } 240 } 241 242 /** 243 * Discrete uniform distribution sampler when the range between lower and upper is too large 244 * to fit in a positive integer. 245 */ 246 private static class LargeRangeDiscreteUniformSampler 247 extends AbstractDiscreteUniformSampler { 248 /** Lower bound. */ 249 private final int lower; 250 /** Upper bound. */ 251 private final int upper; 252 253 /** 254 * @param rng Generator of uniformly distributed random numbers. 255 * @param lower Lower bound (inclusive) of the distribution. 256 * @param upper Upper bound (inclusive) of the distribution. 257 */ 258 LargeRangeDiscreteUniformSampler(UniformRandomProvider rng, 259 int lower, 260 int upper) { 261 super(rng); 262 this.lower = lower; 263 this.upper = upper; 264 } 265 266 @Override 267 public int sample() { 268 // Use a simple rejection method. 269 // This is used when (upper-lower) >= Integer.MAX_VALUE. 270 // This will loop on average 2 times in the worst case scenario 271 // when (upper-lower) == Integer.MAX_VALUE. 272 while (true) { 273 final int r = rng.nextInt(); 274 if (r >= lower && 275 r <= upper) { 276 return r; 277 } 278 } 279 } 280 281 @Override 282 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 283 return new LargeRangeDiscreteUniformSampler(rng, lower, upper); 284 } 285 } 286 287 /** 288 * Adds an offset to an underlying discrete sampler. 289 */ 290 private static class OffsetDiscreteUniformSampler 291 extends AbstractDiscreteUniformSampler { 292 /** The offset. */ 293 private final int offset; 294 /** The discrete sampler. */ 295 private final SharedStateDiscreteSampler sampler; 296 297 /** 298 * @param offset The offset for the sample. 299 * @param sampler The discrete sampler. 300 */ 301 OffsetDiscreteUniformSampler(int offset, 302 SharedStateDiscreteSampler sampler) { 303 super(null); 304 this.offset = offset; 305 this.sampler = sampler; 306 } 307 308 @Override 309 public int sample() { 310 return offset + sampler.sample(); 311 } 312 313 @Override 314 public String toString() { 315 return sampler.toString(); 316 } 317 318 @Override 319 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 320 return new OffsetDiscreteUniformSampler(offset, sampler.withUniformRandomProvider(rng)); 321 } 322 } 323 324 /** 325 * This instance delegates sampling. Use the factory method 326 * {@link #of(UniformRandomProvider, int, int)} to create an optimal sampler. 327 * 328 * @param rng Generator of uniformly distributed random numbers. 329 * @param lower Lower bound (inclusive) of the distribution. 330 * @param upper Upper bound (inclusive) of the distribution. 331 * @throws IllegalArgumentException if {@code lower > upper}. 332 */ 333 public DiscreteUniformSampler(UniformRandomProvider rng, 334 int lower, 335 int upper) { 336 super(null); 337 delegate = of(rng, lower, upper); 338 } 339 340 /** {@inheritDoc} */ 341 @Override 342 public int sample() { 343 return delegate.sample(); 344 } 345 346 /** {@inheritDoc} */ 347 @Override 348 public String toString() { 349 return delegate.toString(); 350 } 351 352 /** 353 * {@inheritDoc} 354 * 355 * @since 1.3 356 */ 357 @Override 358 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 359 // Direct return of the optimised sampler 360 return delegate.withUniformRandomProvider(rng); 361 } 362 363 /** 364 * Creates a new discrete uniform distribution sampler. 365 * 366 * @param rng Generator of uniformly distributed random numbers. 367 * @param lower Lower bound (inclusive) of the distribution. 368 * @param upper Upper bound (inclusive) of the distribution. 369 * @return the sampler 370 * @throws IllegalArgumentException if {@code lower > upper}. 371 * @since 1.3 372 */ 373 public static SharedStateDiscreteSampler of(UniformRandomProvider rng, 374 int lower, 375 int upper) { 376 if (lower > upper) { 377 throw new IllegalArgumentException(lower + " > " + upper); 378 } 379 380 // Choose the algorithm depending on the range 381 382 // Edge case for no range. 383 // This must be done first as the methods to handle lower == 0 384 // do not handle upper == 0. 385 if (upper == lower) { 386 return new FixedDiscreteUniformSampler(lower); 387 } 388 389 // Algorithms to ignore the lower bound if it is zero. 390 if (lower == 0) { 391 return createZeroBoundedSampler(rng, upper); 392 } 393 394 final int range = (upper - lower) + 1; 395 // Check power of 2 first to handle range == 2^31. 396 if (isPowerOf2(range)) { 397 return new OffsetDiscreteUniformSampler(lower, 398 new PowerOf2RangeDiscreteUniformSampler(rng, range)); 399 } 400 if (range <= 0) { 401 // The range is too wide to fit in a positive int (larger 402 // than 2^31); use a simple rejection method. 403 // Note: if range == 0 then the input is [Integer.MIN_VALUE, Integer.MAX_VALUE]. 404 // No specialisation exists for this and it is handled as a large range. 405 return new LargeRangeDiscreteUniformSampler(rng, lower, upper); 406 } 407 // Use a sample from the range added to the lower bound. 408 return new OffsetDiscreteUniformSampler(lower, 409 new SmallRangeDiscreteUniformSampler(rng, range)); 410 } 411 412 /** 413 * Create a new sampler for the range {@code 0} inclusive to {@code upper} inclusive. 414 * 415 * <p>This can handle any positive {@code upper}. 416 * 417 * @param rng Generator of uniformly distributed random numbers. 418 * @param upper Upper bound (inclusive) of the distribution. Must be positive. 419 * @return the sampler 420 */ 421 private static AbstractDiscreteUniformSampler createZeroBoundedSampler(UniformRandomProvider rng, 422 int upper) { 423 // Note: Handle any range up to 2^31 (which is negative as a signed 424 // 32-bit integer but handled as a power of 2) 425 final int range = upper + 1; 426 return isPowerOf2(range) ? 427 new PowerOf2RangeDiscreteUniformSampler(rng, range) : 428 new SmallRangeDiscreteUniformSampler(rng, range); 429 } 430 431 /** 432 * Checks if the value is a power of 2. 433 * 434 * <p>This returns {@code true} for the value {@code Integer.MIN_VALUE} which can be 435 * handled as an unsigned integer of 2^31.</p> 436 * 437 * @param value Value. 438 * @return {@code true} if a power of 2 439 */ 440 private static boolean isPowerOf2(final int value) { 441 return value != 0 && (value & (value - 1)) == 0; 442 } 443 }