/*
cc timerbit.c -o pbit -lm
*/
/* 
timer1.txt:   cc timerbit.c -o pbit -lm         on molokai
timer2.txt:   cc timerbit.c -o pbit -lm -fast   on molokai
      ....50/50 random bits really are that fast! 300 reps of 20K takes 0.21 secs!
timer3.txt:   cc timerbit.c -o pbit -lm         on manoa (compaq laptop)
timer4.txt:   cc timerbit.c -o pbit -lm -O4   on manoa
*/

#include	<stdlib.h>
#include	<stdio.h>
#include	<time.h>

#define	NRUSED	20000
#define	NREPS	50
#define	NSKIP	10
#define	NP 	127
#define	NPMNQ 	(NP-30)
#define	NRANDS	(NRUSED+NP)
#define	NPROBZ	(NRUSED)

long unsigned	rands1[NRANDS];
long unsigned	filt[32];
typedef struct probarray
{
	long unsigned	prob[NPROBZ];
	long unsigned	probfilt;
	int		probbit;
	int		probword;
	double		prob_val;
} Probz;
Probz	probbits;

long unsigned multishift();

main()
{
	int iter, wordnum, bitnum, oldbit, newbit;
	long icnt, yescnt1, yescnt2, yescnt3;
	int seed;
	long	probcnter1;
	long unsigned	*prob1, probfilt1;
	double	drandtime, rbittime, pbittime;
	double bee, cee, dee, eee;

	for(iter=NSKIP;iter<=NREPS;iter+=NSKIP)
	{

		bee = clock();

		seed=1465739853;
		srand48(seed);
		yescnt1 = 0;
		for(icnt=0;icnt<32*NRUSED*iter;icnt++)
			if(drand48()<0.5) yescnt1++;

		cee = clock();	





		seed=1465739853;
		srand48(seed);

		InitFilt();
		InitRands(rands1);
		yescnt2 = 0;
		for(icnt=0;icnt<iter;icnt++)
		{
			for(bitnum=0;bitnum<32;bitnum++)
				for(wordnum=0;wordnum<NRUSED;wordnum++)
					if(rands1[wordnum] & filt[bitnum]) yescnt2++;
			RedoRands(rands1);
		}

		dee = clock();	






		seed=1465739853;
		srand48(seed);
		InitFilt();
		InitRands(rands1);
		probbits.prob_val = 0.5000000;
		InitProb(&probbits);
		yescnt3 = 0;
		for(icnt=0;icnt<iter;icnt++)
		{
			CheckProbz(&probbits,&probcnter1,&prob1,&probfilt1);
			for(bitnum=0;bitnum<32;bitnum++)
			{
				for(wordnum=0;wordnum<NRUSED;wordnum++)
					if(prob1[wordnum] & probfilt1) yescnt3++;
				probcnter1 = 0;
				AdvanceProbz(&probbits,&probcnter1,&probfilt1);
			}
		}

		eee = clock();	

		drandtime = (double)(cee-bee)/(double)CLOCKS_PER_SEC;
		rbittime = (double)(dee-cee)/(double)CLOCKS_PER_SEC;
		pbittime = (double)(eee-dee)/(double)CLOCKS_PER_SEC;


		printf("%d \t %e \t %e \t %e      \t %ld \t %ld \t %ld\n",
  			iter, drandtime, rbittime, pbittime, 
			yescnt1, yescnt2, yescnt3);
		fflush(stdout);

	}


}

InitFilt()
{
	int		i;
	long unsigned	j=1;
	
	for(i=0;i<32;i++,j*=2)
		filt[i] = j;
}

InitRands(array)
long unsigned *array;
{
	int	i, j;
	
	for(i=0;i<NP;i++)
	{
		array[i] = 0x0;
		for(j=0;j<32;j++)
			if(drand48()<0.5) 
				array[i] = array[i] | filt[j];
	}
	for(i=0;i<NRUSED;i++)
		array[NP+i] = array[i] ^ array[NPMNQ+i];
}

RedoRands(array)
long unsigned *array;
{
	int	i;
	long unsigned *pnter, *pnter2;
	
	pnter = &array[NRUSED];
	for(i=0;i<NP;i++)
		array[i] = pnter[i];

	pnter = &array[NP];
	pnter2 = &array[NPMNQ];
	for(i=0;i<NRUSED;i++)
		pnter[i] = array[i] ^ pnter2[i];		
}

InitProb(probbie)
Probz *probbie;
{
	int i, npp, ibit;
	long	iword;
	long unsigned *array;

	array = probbie->prob;

	if(probbie->prob_val==0.0)
	{
		for(i=0;i<NPROBZ;i++)
			array[i] = 0x0;
	}
	else if(probbie->prob_val==1.0)
	{
		for(i=0;i<NPROBZ;i++)
			array[i] = 0xffffffff;
	}
	else
	{
		npp = probbie->prob_val * 32 * NPROBZ;
		
		if(probbie->prob_val >= 0.5)
			npp = 32 * NPROBZ - npp;
		
		for(i=0;i<NPROBZ;i++)
			array[i] = 0x0;
	
		for(i=0;i<npp;i++)
		{
			do
			{
				iword = NPROBZ * drand48();
				ibit = 32 * drand48();
			}
			while((array[iword]&filt[ibit])==filt[ibit]);
			array[iword] |= filt[ibit];
		}
	
		if(probbie->prob_val >= 0.5)
			for(i=0;i<NPROBZ;i++)
				array[i] = ~array[i];
	}
	probbie->probbit = 0;
	probbie->probword = NPROBZ;
	probbie->probfilt = filt[0];
}

CheckProbz(probbie,numleft,probs,filter)
Probz	*probbie;
long	*numleft;
long unsigned  **probs, *filter;
{
	*numleft = probbie->probword;
	*filter = probbie->probfilt;
	*probs = probbie->prob;
}

AdvanceProbz(probbie,numleft,filter)
Probz	*probbie;
long	*numleft;
long unsigned  *filter;
{
	probbie->probword = *numleft;
	if(probbie->probword==0)
	{
		probbie->probbit++;
		if(probbie->probbit==32)
		{
			ScrambleProbz(probbie->prob);
			probbie->probbit = 0;
		}
		probbie->probword = NPROBZ;
		probbie->probfilt = filt[probbie->probbit];
		*numleft = (long)NPROBZ;
		*filter = probbie->probfilt;
	}
}

ScrambleProbz(array)
long unsigned	array[NPROBZ];
{
	int	i, j, 
		ii, jj, 
		isii;
	long unsigned itemp1, itemp2, iswitch;

	for(i=0;i<NPROBZ;i++)
	{
		isii = 32 * drand48() - 16;
		array[i] = multishift(array[i],isii);
			
		jj = NPROBZ * drand48();

		itemp1 = array[i];
		itemp2 = array[jj];
		iswitch = rands1[i] & (itemp1 ^ itemp2);
		array[i]  = itemp2 ^ iswitch;
		array[jj] = itemp1 ^ iswitch;
	}
	RedoRands(rands1);
}

long unsigned multishift(word,bits)	/*  multi-bit shifter  */
long unsigned word;
int bits;
{
	long unsigned temp1, temp2;
	int i;
	
	temp2 = word;
	if(bits>=0)
	{
		for(i=0;i<bits;i++)
		{
			temp1 = temp2 & 0x80000000;
			temp2 = (temp2 << 1) & 0xfffffffe;
			if(temp1==0x80000000)
				temp2 = temp2 | 0x1;
		}
	}
	else
	{
		for(i=0;i<-bits;i++)
		{
			temp1 = temp2 & 0x1;
			temp2 = (temp2 >> 1) & 0x7fffffff;
			if(temp1==0x1)
				temp2 = temp2 | 0x80000000;
		}
	}
	return(temp2);
}