zrand.c

Calc Software Package for Number Calc

	 */
	for (i=0; i < SHUFCNT; ++i) {

		/*
		 * skip values if required
		 *
		 * See:
		 *	Knuth's "The Art of Computer Programming -
		 *	Seminumerical Algorithms", Vol 2, 3rd edition (1998),
		 *	Section 3.6, page 188".
		 */
		if (s100.need_to_skip <= 0) {
			int sk;

			/* skip the require number of values */
			for (sk=0; sk < RAND_SKIP; ++sk) {

				/* bump j and k */
				if (++s100.j >= S100) {
					s100.j = 0;
				}
				if (++s100.k >= S100) {
					s100.k = 0;
				}

				/* slot[k] -= slot[j] */
				SSUB(s100, s100.k, s100.j);

				/* NOTE: don't shuffle, no shuffle table yet */
			}

			/* reset the skip count */
			s100.need_to_skip = RAND_CONSEQ_USE;
			if (conf->calc_debug & CALCDBG_RAND) {
			    printf("rand: skipped %d states\n", RAND_SKIP);
			}

		/* no skip, just descrment use counter */
		} else {
			--s100.need_to_skip;
		}

		/* bump j and k */
		if (++s100.j >= S100) {
			s100.j = 0;
		}
		if (++s100.k >= S100) {
			s100.k = 0;
		}

		/* slot[k] -= slot[j] */
		SSUB(s100, s100.k, s100.j);

		/* shuf[i] = slot[k] */
		SSHUF(s100, i, s100.k);
	}

	/*
	 * note that we are seeded
	 */
	s100.seeded = 1;

	/*
	 * return the previous state
	 */
	return ret;
}


/*
 * zsetrand - set the s100 generator state
 *
 * given:
 *	state - the state to copy
 *
 * returns:
 *	previous s100 state
 */
RAND *
zsetrand(CONST RAND *state)
{
	RAND *ret;		/* previous s100 state */

	/*
	 * initialize state if first call
	 */
	if (!s100.seeded) {
		s100 = init_s100;
	}

	/*
	 * save the current state to return later
	 */
	ret = randcopy(&s100);

	/*
	 * load the new state
	 */
	s100 = *state;

	/*
	 * return the previous state
	 */
	return ret;
}


/*
 * slotcp - copy up to 64 bits from a 64 bit array of FULLs to some HALFs
 *
 * We will copy data from an array of FULLs into an array of HALFs.
 * The destination within the HALFs is some bit location found in bitstr.
 * We will attempt to copy 64 bits, but if there is not enough room
 * (bits not yet loaded) in the destination bit string we will transfer
 * what we can.
 *
 * The src slot is 64 bits long and is stored as an array of FULLs.
 * When FULL_BITS is 64 the element is 1 FULL, otherwise FULL_BITS
 * is 32 bits and the element is 2 FULLs.  The most significant bit
 * in the array (highest bit in the last FULL of the array) is to
 * be transfered to the most significant bit in the destination.
 *
 * given:
 *	bitstr	- most significant destination bit in a bit string
 *	src	- low order FULL in a 64 bit slot
 *	count	- number of bits to transfer (must be 0 < count <= 64)
 *
 * returns:
 *	number of bits transfered
 */
static int
slotcp(BITSTR *bitstr, FULL *src, int count)
{
	HALF *dh;		/* most significant HALF in dest */
	int dnxtbit;		/* next bit beyond most signif in dh */
	int need;		/* number of bits we need to transfer */
	int ret;		/* bits transfered */

	/*
	 * determine how many bits we actually need to transfer
	 */
	dh = bitstr->loc;
	dnxtbit = bitstr->bit+1;
	count &= (SBITS-1);
	need = (bitstr->len < count) ? bitstr->len : count;

	/*
	 * prepare for the return
	 *
	 * Note the final bitstr location after we have moved the
	 * position down 'need' bits.
	 */
	bitstr->len -= need;
	bitstr->loc -= need / BASEB;
	bitstr->bit -= need % BASEB;
	if (bitstr->bit < 0) {
		--bitstr->loc;
		bitstr->bit += BASEB;
	}
	ret = need;

	/*
	 * deal with aligned copies quickly
	 */
	if (dnxtbit == BASEB) {
		if (need == SBITS) {
#if 2*FULL_BITS == SBITS
			*dh-- = (HALF)(src[SLEN-1] >> BASEB);
			*dh-- = (HALF)(src[SLEN-1]);
#endif
			*dh-- = (HALF)(src[0] >> BASEB);
			*dh = (HALF)(src[0]);
#if 2*FULL_BITS == SBITS
		} else if (need > FULL_BITS+BASEB) {
			*dh-- = (HALF)(src[SLEN-1] >> BASEB);
			*dh-- = (HALF)(src[SLEN-1]);
			*dh-- = (HALF)(src[0] >> BASEB);
			*dh = ((HALF)src[0] &
				     highhalf[need-FULL_BITS-BASEB]);
		} else if (need > FULL_BITS) {
			*dh-- = (HALF)(src[SLEN-1] >> BASEB);
			*dh-- = (HALF)(src[SLEN-1]);
			*dh = ((HALF)(src[0] >> BASEB) &
				     highhalf[need-FULL_BITS]);
#endif
		} else if (need > BASEB) {
			*dh-- = (HALF)(src[SLEN-1] >> BASEB);
			*dh = ((HALF)(src[SLEN-1]) & highhalf[need-BASEB]);
		} else {
			*dh = ((HALF)(src[SLEN-1] >> BASEB) & highhalf[need]);
		}
		return ret;
	}

	/*
	 * load the most significant HALF
	 */
	if (need >= dnxtbit) {
		/* fill up the most significant HALF */
		*dh-- |= (HALF)(src[SLEN-1] >> (FULL_BITS-dnxtbit));
		need -= dnxtbit;
	} else if (need > 0) {
		/* we exhaust our need before 1st half is filled */
		*dh |= (HALF)((src[SLEN-1] >> (FULL_BITS-need)) <<
				(dnxtbit-need));
		return ret;	/* our need has been filled */
	} else {
		return ret;	/* our need has been filled */
	}

	/*
	 * load the 2nd most significant HALF
	 */
	if (need > BASEB) {
		/* fill	 up the 2nd most significant HALF */
		*dh-- = (HALF)(src[SLEN-1] >> (BASEB-dnxtbit));
		need -= BASEB;
	} else if (need > 0) {
		/* we exhaust our need before 2nd half is filled */
		*dh |= ((HALF)(src[SLEN-1] >> (BASEB-dnxtbit)) &
			 highhalf[need]);
		return ret;	/* our need has been filled */
	} else {
		return ret;	/* our need has been filled */
	}

	/*
	 * load the 3rd most significant HALF
	 *
	 * At this point we know that our 3rd HALF will force us
	 * to cross into a second FULL for systems with 32 bit FULLs.
	 * We know this because the aligned case has been previously
	 * taken care of above.
	 *
	 * For systems that have 64 bit FULLs (and 32 bit HALFs) this
	 * is will be our least significant HALF.  We also know that
	 * the need must be < BASEB.
	 */
#if FULL_BITS == SBITS
	*dh |= (((HALF)src[0] & highhalf[dnxtbit+need]) << dnxtbit);
#else
	if (need > BASEB) {
		/* load the remaining bits from the most signif FULL */
		*dh-- = ((((HALF)src[SLEN-1] & lowhalf[BASEB-dnxtbit])
			   << dnxtbit) | (HALF)(src[0] >> (FULL_BITS-dnxtbit)));
		need -= BASEB;
	} else if (need > 0) {
		/* load the remaining bits from the most signif FULL */
		*dh-- |= (((((HALF)src[SLEN-1] & lowhalf[BASEB-dnxtbit])
			 << dnxtbit) | (HALF)(src[0] >> (FULL_BITS-dnxtbit))) &
			  highhalf[need]);
		return ret;	/* our need has been filled */
	} else {
		return ret;	/* our need has been filled */
	}

	/*
	 * load the 4th most significant HALF
	 *
	 * At this point, only 32 bit FULLs are operating.
	 */
	if (need > BASEB) {
		/* fill	 up the 2nd most significant HALF */
		*dh-- = (HALF)(src[0] >> (BASEB-dnxtbit));
		/* no need todo: need -= BASEB, because we are nearly done */
	} else if (need > 0) {
		/* we exhaust our need before 2nd half is filled */
		*dh |= ((HALF)(src[0] >> (BASEB-dnxtbit)) &
			 highhalf[need]);
		return ret;	/* our need has been filled */
	} else {
		return ret;	/* our need has been filled */
	}

	/*
	 * load the 5th and least significant HALF
	 *
	 * At this point we know that the need will be satisfied.
	 */
	*dh |= (((HALF)src[0] & lowhalf[BASEB-dnxtbit]) << dnxtbit);
#endif
	return ret;	/* our need has been filled */
}


/*
 * slotcp64 - copy 64 bits from a 64 bit array of FULLs to some HALFs
 *
 * We will copy data from an array of FULLs into an array of HALFs.
 * The destination within the HALFs is some bit location found in bitstr.
 * Unlike slotcp(), this function will always copy 64 bits.
 *
 * The src slot is 64 bits long and is stored as an array of FULLs.
 * When FULL_BITS is 64 this array is 1 FULL, otherwise FULL_BITS
 * is 32 bits and the array is 2 FULLs.	 The most significant bit
 * in the array (highest bit in the last FULL of the array) is to
 * be transfered to the most significant bit in the destination.
 *
 * given:
 *	bitstr	- most significant destination bit in a bit string
 *	src	- low order FULL in a 64 bit slot
 *
 * returns:
 *	number of bits transfered
 */
static void
slotcp64(BITSTR *bitstr, FULL *src)
{
	HALF *dh = bitstr->loc;		/* most significant HALF in dest */
	int dnxtbit = bitstr->bit+1;	/* next bit beyond most signif in dh */

	/*
	 * prepare for the return
	 *
	 * Since we are moving the point 64 bits down, we know that
	 * the bit location (bitstr->bit) will remain the same.
	 */
	bitstr->len -= SBITS;
	bitstr->loc -= SBITS/BASEB;

	/*
	 * deal with aligned copies quickly
	 */
	if (dnxtbit == BASEB) {
#if 2*FULL_BITS == SBITS
		*dh-- = (HALF)(src[SLEN-1] >> BASEB);
		*dh-- = (HALF)(src[SLEN