// Copyright (2014) Staci White

__global__ void BF(double C, double epsilonbf, double gammabf, double
m, const int H, const double * Brv, const double *
unif1, const double * unif2, double * Bpvec){ 

	   double k = m/(gammabf*epsilonbf);
	   int i = 1;
	   int R = 1;
	   int n1 = 0;
	   int n2 = 0;
	   int bidx = blockIdx.x;

        while( i>0 && R > 0){

	      while( i>0 && i<k){

	       	      n1+=1; 

		      int G=lrint(floor(log(unif1[bidx+H*(n1-1)])/log(1/C)))+1;

       		     i+=(1-Brv[bidx+H*(n1-1)])*G-1;

	 	}

			
	     if(i>=k){

		n2+=1;

		if(unif2[bidx+H*(n2-1)]>exp(-i*log(1+gammabf*epsilonbf))){

			R=0;

		}

		C=C*(1+gammabf*epsilonbf);
		epsilonbf = (1-gammabf)*epsilonbf;
		k = k/(1-gammabf);

	     }

	}


	if(i==0){
			
		Bpvec[bidx]=1;

	}



}