view release on metacpan or  search on metacpan

mtwist/randistrs.c  view on Meta::CPAN

#define RD_UNIFORM_THRESHOLD	((int)((double)(1u << 31) * 2.0 * RD_MAX_BIAS))
#endif /* RD_UNIFORM_THRESHOLD */

#if NVMANTBITS <= 32
#  define mts_ldrand(x) mts_drand(x)
#elif NVMANTBITS <= 64
#else
#  define mts_ldrand(x) mts_lldrand(x)
#endif

/*
 * Generate a uniform integer distribution on the open interval
 * [lower, upper).  See comments above about RD_UNIFORM_THRESHOLD.  If
 * we are above the threshold, this function is relatively expensive
 * because we may have to repeatedly draw random numbers to get a
 * one that works.
 */
int32_t rds_iuniform(
    mt_state *		state,		/* State of the MT PRNG to use */
    int32_t		lower,		/* Lower limit of distribution */
    int32_t		upper)		/* Upper limit of distribution */
    {
    uint32_t		range = upper - lower;
					/* Range of requested distribution */

    if (range <= RD_UNIFORM_THRESHOLD)
	return lower + (int32_t)(mts_ldrand(state) * range);
    else
	{
	/*
	 * Using the simple formula would produce too much bias.
	 * Instead, draw numbers until we get one within the range.
	 * To save time, we first calculate a mask so that we only
	 * look at the number of bits we actually need.  Since finding
	 * the mask is expensive, we optionally do a bit of caching
	 * here (note that the caching makes the code non-reentrant;
	 * set MT_CACHING to turn on this misfeature).
	 *
	 * Incidentally, the astute reader will note that we use the
	 * low-order bits of the PRNG output.  If the PRNG were linear
	 * congruential, using the low-order bits wouuld be a major
	 * no-no.  However, the Mersenne Twist PRNG doesn't have that
	 * drawback.
	 */
#ifdef MT_CACHING
	static uint32_t	lastrange = 0;	/* Range used last time */
	static uint32_t	rangemask = 0;	/* Mask for range */
#else /* MT_CACHING */
	uint32_t	rangemask = 0;	/* Mask for range */
#endif /* MT_CACHING */
	register uint32_t
			ranval;		/* Random value from mts_lrand */

#ifdef MT_CACHING
	if (range != lastrange)
#endif /* MT_CACHING */
	    {
	    /*
	     * Range is different from last time, recalculate mask.
	     *
	     * A few iterations could be trimmed off of the loop if we
	     * started rangemask at the next power of 2 above
	     * RD_UNIFORM_THRESHOLD.  However, I don't currently know
	     * a formula for generating that value (though there is
	     * probably one in HAKMEM).
	     */
#ifdef MT_CACHING
	    lastrange = range;
#endif /* MT_CACHING */
	    for (rangemask = 1;
	      rangemask < range  &&  rangemask != 0;
	      rangemask <<= 1)
		;

	    /*
	     * If rangemask became zero, the range is over 2^31.  In
	     * that case, subtracting 1 from rangemask will produce a
	     * full-word mask, which is what we need.
	     */
	    rangemask -= 1;
	    }

	/*
	 * Draw random numbers until we get one in the requested range.
	 */
	do
	    {
	    ranval = mts_lrand(state) & rangemask;
	    }
	    while (ranval >= range);
	return lower + ranval;
	}
    }

#ifdef INT64_MAX
/*
 * Generate a uniform integer distribution on the half-open interval
 * [lower, upper).
 */
int64_t rds_liuniform(
    mt_state *		state,		/* State of the MT PRNG to use */
    int64_t		lower,		/* Lower limit of distribution */
    int64_t		upper)		/* Upper limit of distribution */
    {
    uint64_t		range = upper - lower;
					/* Range of requested distribution */

    /*
     * Draw numbers until we get one within the range.  To save time,
     * we first calculate a mask so that we only look at the number of
     * bits we actually need.  Since finding the mask is expensive, we
     * optionally do a bit of caching here.  See rds_iuniform for more
     * information.
     */
#ifdef MT_CACHING
    static uint32_t	lastrange = 0;	/* Range used last time */
    static uint32_t	rangemask = 0;	/* Mask for range */
#else /* MT_CACHING */
    uint32_t		rangemask = 0;	/* Mask for range */
#endif /* MT_CACHING */
    register uint32_t	ranval;		/* Random value from mts_lrand */