cryrand.cal

早期freebsd实现

 *
 *    This is a reserved seed for sshufrand(0) and srand(0).  One should
 *    not directly call srand(-1).
 *
 * else:
 *
 *    Reserved for future use.
 *
 * Anyone comfortable with seed == 0 should also be comfortable with
 * non-zero seeds.  A non-zero seeded sequence will produce a values
 * that have the exact same pseudo-random properties as the algorithm
 * described by Knuth.  I.e., the sequence, while different, is as good
 * (or bad) as the sequence produced by a seed of 0.
 *
 * This function updates both the cry_seed and a55_seed64 global values.
 *
 * This function is called from the shuffle generator seed function sshufrand().
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you need a fast generator and do not need a crypto strong one,
 *       you should consider using the shuffle generator (see sshufrand()
 *	 and srand()).  Direct use of this function is not recommended.
 */
define
sa55rand(seed)
{
    local oldseed;		/* previous seed */
    local junk;			/* discards the first few numbers */
    local j;

    /* firewall */
    if (!isint(seed)) {
	quit "bad arg: arg must be an integer";
    }
    if (seed < -1) {
	quit "bad arg: seed < 0 is reserved for future use";
    }

    /* misc table setup */
    oldseed = cry_seed;				/* remember the previous seed */
    cry_seed = seed;				/* save the new seed */
    if (cry_seed == -1) {
    	/* since -1 was a special case, pretend it really was zero */
    	cry_seed = 0;
    }
    add55_tbl = add55_init_tbl;	/* reload the table */
    add55_j = 23;               /* reset first walking table pointer */
    add55_k = 54;               /* reset second walking table pointer */

    /* obtain our 64 bit xor seed */
    add55_seed64 = randreseed64(cry_seed);

    /* spin the pseudo-random number generator a while */
    if (seed == 0) {
    	/* we will act as if sshufrand(0) or srand(0) had been called */
	for (j=0; j < 513; ++j) {
	    junk = a55rand();
	}
    } else {
	for (j=0; j < 128; ++j) {
	    junk = a55rand();
	}
    }

    /* return the old seed */
    return(oldseed);
}


/*
 * shufrand - implement the shuffle pseudo-random number generator
 *
 * returns:
 *	A number in the half open interval [0,2^64)
 *
 * This function implements the shuffle number generator.
 *
 * This is a generator based on the shuffle generator as described in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
 *
 * The function rand() calls this function.
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator,
 *	 this function may very well serve your needs.
 */
define
shufrand()
{
    local j;		/* table index to replace */

    /*
     * obtain a new random value
     * determine the table entry to shuffle
     * shuffle out the value we will return
     */
    shuf_y = shuf_tbl[j = (shuf_y >> shuf_shift)];

    /* shuffle in the new random value */
    shuf_tbl[j] = a55rand();

    /* return the shuffled out value */
    return (shuf_y);
}


/*
 * sshufrand - seed the shuffle pseudo-random generator
 *
 * given:
 *	a seed
 *
 * returns:
 *	the previous seed
 *
 * This function implements the shuffle number generator.
 *
 * This is a generator based on the shuffle generator as described in
 * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
 * vol 2, 2nd edition (1981), Section 3.2.2, page 32, Algorithm B.
 *
 * The low 64 bits of seed are used by the additive 55 generator.
 * See the sa55rand() function for details.  The remaining bits of seed
 * are used to perform an initial shuffle on the shuffle state table.
 * The size of the seed also determines the size of the shuffle table.
 *
 * The shuffle table size is always a power of 2, and is at least 128
 * entries long.  Let the table size be:
 *
 *	shuf_size = 2^shuf_pow
 *
 * The number of ways one could shuffle that table is:
 *
 *	(2^shuf_pow)!
 *
 * We select the smallest 'shuf_pow' (and thus the size of the shuffle table)
 * such that the following are true:
 *
 *	(2^shuf_pow)! >= (seed / 2^64)    and    2^shuf_pow >= 128
 *
 * Given that we now have the table size of 'shuf_size', we must go about
 * loading the table and shuffling it.
 *
 * Loading is easy, we will generate random values via the additive 55
 * generator (a55rand()) and load them into successive entries.
 *
 * We enumerate the (2^shuf_pow)! values via:
 *
 *	shuf_seed = 2*(U[2] + 3*(U[3] + 4*(U[4] + ...
 *			  + (U[shuf_pow-1]*(shuf_pow-1)) ... )))
 *	0 <= U[j] < j
 *
 * We swap the swap table entries shuf_tbl[U[j]] & shuf_tbl[j-1] for all
 * 1 < j < shuf_pow.
 *
 * We will convert 'seed / 2^64' into shuf_seed, by applying the 64 bit
 * scramble function on 64 bit chunks of 'seed / 2^64'.
 *
 * The function srand() calls this function.
 *
 * There is no limit on the size of a seed.  On the other hand,
 * extremely large seeds require large tables and long seed times.
 * Using a seed in the range of [2^64, 2^64 * 128!) should be
 * sufficient for most purposes.  An easy way to stay within this
 * range to to use seeds that are between 21 and 215 digits long, or
 * 64 to 780 bits long.
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator,
 *	 this function may very well serve your needs.
 */
define
sshufrand(seed)
{
    local shuf_pow;		/* power of two - log2(shuf_size) */
    local shuf_seed;		/* seed bits beyond low 64 bits */
    local oldseed;		/* previous seed */
    local shift;		/* shift factor to access 64 bit chunks */
    local swap_indx;		/* what to swap shuf_tbl[0] with */
    local rval;			/* random value form additive 55 generator */
    local j;

    /* firewall */
    if (!isint(seed)) {
	quit "bad arg: seed must be an integer";
    }
    if (seed < 0) {
	quit "bad arg: seed < 0 is reserved for future use";
    }

    /*
     * seed the additive 55 generator
     */
    if (seed == 0) {
	/* allow sshufrand(0) and srand(0) to arrive at the same state */
	oldseed = sa55rand(-1);
    } else {
	oldseed = sa55rand(seed);
    }

    /*
     * form the shuffle table size and constants
     */
    shuf_seed = seed >> 64;
    for (shuf_pow = 7; shuf_seed > (j=fact(1<<(shuf_pow))) && shuf_pow < 64; \
	 ++shuf_pow) {
    }
    shuf_size = (1 << shuf_pow);
    shuf_shift = 64 - shuf_pow;
    /* reallocate the shuffle table */
    mat shuf_tbl[shuf_size];

    /*
     * scramble the seed above the low 64 bits
     */
    if (shuf_seed > 0) {
	j = 0;
	for (shift=0; shift < highbit(shuf_seed)+1; shift += 64) {
	    j |= (randreseed64(shuf_seed >> shift) << shift);
	}
	shuf_seed = j;
    }

    /*
     * load the shuffle table
     */
    for (j=0; j < shuf_size; ++j) {
	shuf_tbl[j] = a55rand();
    }
    shuf_y = a55rand();		/* get the next Y value */

    /*
     * We will shuffle based the process outlined in:
     *
     * Knuth's "The Art of Computer Programming - Seminumerical Algorithms",
     * vol 2, 2nd edition (1981), Section 3.4.2, page 139, Algorithm P.
     *
     * Here, we will let j run over the range [0,shuf_size) instead of
     * [shuf_size,0) as outlined in algorithm P.  We will also generate
     * U values from shuf_seed.
     */
    j = 0;
    while (shuf_seed > 0 && ++j < shuf_size) {

	/* determine what index we will swap with the '0' index */
	quomod(shuf_seed, (j+1), shuf_seed, swap_indx);

	/* swap table entries, if needed */
	if (swap_indx != j) {
	    swap(shuf_tbl[j], shuf_tbl[swap_indx]);
	}
    }

    /*
     * run the generator for twice the table size
     */
    for (j=0; j < shuf_size*2; ++j) {
	rval = shufrand();
    }

    /* return the old seed */
    return (oldseed);
}


/*
 * rand - generate a pseudo-random value over a given range via additive 55
 *
 * usage:
 *	rand() 		- generate a pseudo-random integer >=0 and < 2^64
 *	rand(a) 	- generate a pseudo-random integer >=0 and < a
 *	rand(a,b)	- generate a pseudo-random integer >=a and <= b
 *
 * returns:
 *	a large pseudo-random integer over a give range (see usage)
 *
 * When no arguments are given, a pseudo-random number in the half open
 * interval [0,2^64) is produced.  This form is identical to calling
 * shufrand().
 *
 * When 1 argument is given, a pseudo-random number in the half open interval
 * [0,a) is produced.
 *
 * When 2 arguments are given, a pseudo-random number in the closed interval
 * [a,b] is produced.
 *
 * The input values a and b, if given, must be integers.
 *
 * This generator is simply a different interface to the shuffle generator.
 * calling shufrand(), or seeding via sshufrand() will affect the output
 * of this function.
 *
 * NOTE: Unlike cryrand(), this function does not preserve unused bits for
 *	 use by the next function call.
 *
 * NOTE: The Un*x rand() function returns only 16 bit or 31 bits, while we
 *	 return a number of any given size (default is 64 bits).
 *
 * NOTE: This is NOT Blum's method.  This is used by Blum's method to
 *       generate some internal values.
 *
 * NOTE: If you do not need a crypto strong pseudo-random generator
 *	 this function may very well serve your needs.
 */
define
rand(a,b)
{
    local range;		/* we must generate [0,range) first */
    local offset;		/* what to add to get a adjusted range */
    local ret;			/* pseudo-random value */
    local fullwords;		/* number of full 64 bit chunks in ret */
    local finalmask;		/* mask of bits in final chunk of range */
    local j;

    /*
     * setup and special cases
     */
    /* deal with the rand() case */
    if (isnull(a) && isnull(b)) {
	/* no args means same range as additive 55 generator */
	return(a55rand());
    }
    /* firewall - args, if given must be in integers */
    if (!isint(a) || (!isnull(b) && !isint(b))) {
	quit "bad args: args, if given, must be integers";
    }
    /* convert rand(x) into rand(0,x-1) */
    if (isnull(b)) {
	/* convert call into a closed interval */
	b = a-1;
	a = 0;
	/* firewall - rand(0) should act like rand(0,0) */
	if (b == -1) {
	    return(0);
	}
    }
    /* determine the range and offset */
    if (a >= b) {
	/* deal with the case of rand(a,a) */
	if (a == b) {
	    /* not very random, but it is true! */
	    return(a);
	}
	range = a-b+1;
	offset = b;
    } else {
	/* convert rand(a,b), where a < b into rand(b,a) */
	range = b-a+1;
	offset = a;
    }
    /*
     * At this point, we seek a pseudo-random number [0,range) to which
     * we will add offset to produce a number [offset,range+offset).
     */
    fullwords = highbit(range-1)//64;
    finalmask = (1 << (1+(highbit(range-1)%64)))-1;

    /*
     * loop until we get a value that is in range
     *
     * A note in modulus biasing:
     *
     * We will not fall into the trap of thinking that we can simply take
     * a value mod 'range'.  Consider the case where 'range' is '80'
     * and we are given pseudo-random numbers [0,100).  If we took them
     * mod 80, then the numbers [0,20) would be produced more frequently
     * because the numbers [81,100) mod 80 wrap back into [0,20).
     */
    do {
	/* load up all lower full 64 bit chunks with pseudo-random bits */
	ret = 0;
	for (j=0; j < fullwords; ++j) {