#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <assert.h>

#define BITNUM						(ull)8
#define BITNUM_MASK					7
#define BITNUM_SHIFT				3
#define DEFAULT_WORKARR_NUM			(ull)700000
#define DEFAULT_DISPARR_NUM			(unsigned long)100000

typedef unsigned long long ull;

ull findprimes_upthrough(ull limit, char *workarr, ull walen, ull **parr_out)
{
	ull i, j, sq, left, *parr, wablen;
	
	//initiate
	sq = (ull)ceill(sqrtl((long double)limit));
	wablen = walen * BITNUM;
	
	//clear
	for(i = 0; i < walen; i++) workarr[i] = 0;
	
	//we only need to cancel out through the square root
	for(i = 2, left = limit - 1; i <= sq; i++)
	{
		if(workarr[i >> BITNUM_SHIFT] & 1 << (i & BITNUM_MASK)) continue;
		
		for(j = i * 2; j <= limit; j += i)
		{
			if(workarr[j >> BITNUM_SHIFT] & 1 << (j & BITNUM_MASK)) continue;
			
			left--;
			workarr[j >> BITNUM_SHIFT] |= 1 << (j & BITNUM_MASK);
		}
	}
	
	parr = malloc(left * sizeof(ull));
	assert(parr);
	
	//copy prime data
	for(i = 2, j = 0; i <= limit; i++)
	{
		if(workarr[i >> BITNUM_SHIFT] & 1 << (i & BITNUM_MASK)) continue;
		
		parr[j++] = i;
	}
	
	*parr_out = parr;
	
	return j;
}

//j = offset + num % fact[i]
ull findprimes(FILE *fd, ull start, ull end)
{
	ull sq, cf, j, i, lim, total, wablen, walen, offset, pflen, untildlen, *pfarr;
	unsigned long curdisp, idisp, curdlen;
	char *workarr, *disparr;
	
	//set up all the info
	sq = (ull)ceill(sqrtl((long double)end));
	wablen = sq > DEFAULT_WORKARR_NUM ? sq : DEFAULT_WORKARR_NUM;
	walen = wablen / BITNUM + 1;
	
	//allocate
	workarr = malloc(walen);
	assert(workarr);
	
	//set up our display array (to minimize i/o)
	//calc max number of digits in each number, + 1 for newline
	curdlen = (unsigned long)ceill(logl((long double)end + 1) / logl((long double)10)) + 1;
	disparr = malloc(curdlen * DEFAULT_DISPARR_NUM + 1);
	assert(disparr);
	
	//find primes <= square root (in other words, the prime factors)
	pflen = findprimes_upthrough(sq, workarr, walen, &pfarr);
	printf("found prime factors, starting on %lu segment(s) of %lu\n", 
		(unsigned long)ceill((long double)(end - start + 1) / (long double)wablen), (unsigned long)wablen);
	
	if(pfarr[pflen - 1] >= start)
	{
		//find position at which this occurs
		for(i = 0; i < pflen && pfarr[i] < start; i++) ;
		
		//print
		//these are few enough in number we don't need to bother with a buffer array
		for(; i < pflen; i++)
			fprintf(fd, "%s\n", ulltoa(pfarr[i], disparr, 10));
		
		//reset range
		start = pfarr[pflen - 1] + 1;
	}
	
	//initialize display array info
	curdisp = idisp = 0;
	
	//we need to know how far away 'start' is from the subsequent power of ten
	for(untildlen = 1, curdlen = 0; untildlen <= start; untildlen *= 10, curdlen++) ;
	untildlen -= start;
	
	//now step through segments
	for(offset = start, total = 0; offset <= end; offset += wablen)
	{
		//clear workarr
		//for(i = 0; i < walen; i++) workarr[i] = 0;
		memset(workarr, 0, walen);
		
		//cancel out all multiples of each factor
		for(i = 0; i < pflen; i++)
		{
			cf = pfarr[i];
			
			if(cf * cf > offset + wablen) break;
			//j = (cf - offset % cf) % cf;
			j = offset % cf;
			if(j) j = cf - j;
			
			for(; j < wablen; j += cf)
				workarr[j >> BITNUM_SHIFT] |= 1 << (j & BITNUM_MASK);
			
		}
		
		//in logarithmic steps
		for(i = 0; i < wablen; )
		{
			//move up to next power of ten
			if(!untildlen)
			{
				untildlen = 9 * (ull)pow(10, curdlen);
				curdlen++;
			}
			
			//set limit
			if(untildlen + i < wablen)
			{
				lim = i + untildlen;
				untildlen = 0;
			}
			else
			{
				lim = wablen;
				untildlen -= wablen;
			}
			
			//print primes
			for(; i < lim; i++)
			{
				//not a prime
				if(workarr[i >> BITNUM_SHIFT] & 1 << (i & BITNUM_MASK)) continue;
				
				//bigger than end
				if((j = i + offset) > end) goto continue_mainloop;
				
				//print if we've reached our limit
				if(curdisp == DEFAULT_DISPARR_NUM)
				{
					disparr[idisp] = 0;
					fprintf(fd, "%s", disparr);
					
					idisp = curdisp = 0;
				}
				
				//convert to string
				ulltoa(j, &disparr[idisp], 10);
				
				//move our pointer along
				idisp += curdlen;
				disparr[idisp++] = '\n';
				
				//fputs(ulltoa(j, disparr, 10), fd);
				
				curdisp++;
				total++;
			}
			
		}
		
		continue_mainloop:
		continue;
	}
	
	//print what's left over
	disparr[idisp] = 0;
	fprintf(fd, "%s", disparr);
	
	free(disparr);
	free(pfarr);
	free(workarr);
	
	return total;
}

int main(int argc, char *argv[])
{
	unsigned long total;
	clock_t t1, t2;
	FILE *fd;
	
	if(argc != 4) return 0;
	t1 = clock();
	
	fd = fopen(argv[3], "wb");
	assert(fd);
	printf("file opened\n");
	
	total = (unsigned long)findprimes(fd, atoll(argv[1]), atoll(argv[2]));
	printf("total: %lu\n", total);
	
	fclose(fd);
	t2 = clock();
	
	printf("time elapsed: %f\n", (double)(t2 - t1) / (double)CLOCKS_PER_SEC);
	
	return 0;
}