std::vector<unsigned int> prime;

void atkin(unsigned int limit)
{
	unsigned int* bitmap = (unsigned int*) calloc(limit / 32 + 1 + ((limit & 31) ? 1 : 0), sizeof(unsigned int));
	bitmap[0] = 12;
	unsigned int sqr_lim = (unsigned int)sqrt((long double)limit);
 
	// Предположительно простые - это целые с нечетным числом
	// представлений в данных квадратных формах.
	// x2 и y2 - это квадраты i и j (оптимизация).
	unsigned int x2 = 0;
	for (unsigned int i = 1; i <= sqr_lim; i++) {
		x2 += 2 * i - 1;
		unsigned int y2 = 0;
		for (unsigned int j = 1; j <= sqr_lim; j++) {
			y2 += 2 * j - 1;
	 
			unsigned int n = 4 * x2 + y2;
			if ((n <= limit) && (n % 12 == 1 || n % 12 == 5))
				bitmap[n >> 5] ^= (1 << (n & 31));
	 
			n -= x2; // Оптимизация n = 3 * x2 + y2; 
			if ((n <= limit) && (n % 12 == 7)) {
				bitmap[n >> 5] ^= (1 << (n & 31));
			}
	 
			n -= 2 * y2; // Оптимизация n = 3 * x2 - y2;
			if ((i > j) && (n <= limit) && (n % 12 == 11)) {
				bitmap[n >> 5] ^= (1 << (n & 31));
			}
		}
	}
	 
	// Отсеиваем кратные квадратам простых чисел в интервале [5, sqrt(limit)].
	// (основной этап не может их отсеять)
	for (unsigned int i = 5; i <= sqr_lim; i++) {
		if (bitmap[i >> 5] & (1 << (i & 31))) {
			unsigned int n = i * i;
			for (unsigned int j = n; j <= limit; j += n) {
				bitmap[j >> 5] &= ~(1 << (j & 31));
			}
		}
	}
	 
	// Вывод списка простых чисел в консоль.
	//printf("2, 3, 5"); 
	prime.push_back(2);
	prime.push_back(3);
	prime.push_back(5);
	for(unsigned int i = 7; i <= limit; i += 2) {  
		if ((bitmap[i >> 5] & (1 << (i & 31)))){
		  // printf(", %d", i);
		   prime.push_back(i);
		}
	}
	free(bitmap);
}

unsigned int next_prime(unsigned int x)
{
	static std::vector<unsigned int> cache;
	static unsigned int max_prime = 0;
	if(max_prime == 0) {
		cache.push_back(2);
		cache.push_back(3);
		max_prime = 3;
	}
	if(x >= max_prime) { // нет в кэше
		unsigned int limit = (x < 1000) ? 1024 : x /4 * 5; // x+20%

		unsigned int start = max_prime + 2; // должно быть четное
		unsigned int* bitmap = (unsigned int*) calloc((limit - start) / 64 + 1, sizeof(unsigned int));
		limit -= start;

		if(max_prime > 3) { // инициализация предыдущим расчетом
			unsigned int* end = &cache[cache.size() - 1];
			for(unsigned int* cur = &cache[1];cur <= end; cur++) {
				unsigned int i = *cur * *cur; 
				if(i > limit + start) break;
				unsigned int step = *cur << 1;
				if(i < start) { // выравнивание под cur*cur+2N*cur
					i = start - (start - i) % step;
					if(i < start) i += step;
				}
				i -= start;
				//printf("start=%d i=%d max=%d\n", start, i, max_prime);
				for(; i < limit; i += step) { // Вычеркиваем кратные max_prime
					bitmap[i >> 6] |= (1 << ((i >> 1) & 31));
				}
			}
		}

		bool need_fill = true;
		while(need_fill) {
			if(max_prime >= 65536 || max_prime * max_prime >= limit + start) {
				need_fill = false; // дальше заполнять не надо
			} else {
				unsigned int step = max_prime << 1;
				for(unsigned int i = max_prime * max_prime - start; i < limit; i += step) { // Вычеркиваем кратные max_prime
					bitmap[i >> 6] |= (1 << ((i >> 1) & 31));
				}
			}
			// вычитываем следущую порцию
			for(unsigned int i = max_prime + 2 - start; i < limit; i += 2) {
				if((bitmap[i >> 6] & (1 << ((i >> 1) & 31))) == 0) {
					cache.push_back(i + start);
					if(need_fill) {
						break;
					}
				}
			}
			max_prime = cache.back();
		}
		free(bitmap);
	}

	static unsigned int prev_id = 0;
	static unsigned int prev_prime = 0;
	if(x < 2) {
		if(x == 0) { // освобождение памяти
			max_prime = 0;
			cache.clear();
		}
		prev_id = 0;
		prev_prime = 2;
	} else if(x == prev_prime) {
		prev_id++;
		prev_prime = cache[prev_id];
	} else { // поиск в кэше
		unsigned int* min = &cache[0];
		unsigned int* max = &cache[cache.size() - 1];
		while(max - min > 1) {
			unsigned int* mid = min + (max - min) / 2;
			if (*mid > x) {
				max = mid;
			} else {
				min = mid;
			}
		}
		prev_id = max - &cache[0];
		prev_prime = *max;
	}
	return prev_prime;
}

// Пример использования
int k = 1;
while(k < 1000) {
   k = next_prime(k);
   printf("%d ", k);
}

unsigned int next_prime(unsigned int x)
{
	static unsigned int* bitmap = NULL; // биткарта
	static unsigned int bm_size = 0; // размер биткарты в байтах
	static unsigned int limit = 0; // до скольки заполнять решето

	unsigned int ret = 0;
	while(true) {
		if(x < 3) {
			if(x == 0) { // освобождение памяти
				limit = 0;
				bm_size = 0;
				if(bitmap) {
					free(bitmap);
					bitmap = NULL;
				}
			}
			ret = (x < 2) ? 2 : 3;
			break;
		}
		// поиск следующего после x
		for(unsigned int i = (x | 1) + 2; i < limit; i += 2) {
			if((bitmap[i >> 6] & (1 << ((i >> 1) & 31))) == 0) {
				ret = i;
				break;
			}
		}
		if(ret != 0) break;
		// нехватило биткарты, дорасчет
		unsigned int old_limit = limit;
		limit += 1048576; // 64Kb
		if(limit/16 >= bm_size) { // довыделение памяти
			unsigned int bm_size_new = bm_size / 4 * 5;
			if(limit/16 > bm_size_new) bm_size_new = limit/16;
			if((x + 1048576)/16 > bm_size_new) bm_size_new = ((x + 1048576)/16) + 4; 
			if(bm_size_new > 0x10000000) bm_size_new = 0x10000000; // 2^(sizeof(unsigned int)/16)
			unsigned int* p = (unsigned int*) realloc(bitmap, bm_size_new); 
			if(!p) {
				printf("No memory. New size %d byte\n", bm_size_new);
				ret = 0xFFFFFFFF;
				break;
			}
			bitmap = p;
			memset(bitmap + bm_size/4, 0, bm_size_new - bm_size);
			bm_size = bm_size_new;
		}
		unsigned int max_prime = 3; // максимальное простое
		while(true) {
			unsigned int i = max_prime * max_prime; 
			if(i >= limit) break; // дальше заполнять не надо
			unsigned int step = max_prime << 1;
			if(i < old_limit) { // выравнивание под max_prime*max_prime+2N*max_prime
				i = old_limit - (old_limit - i) % step;
				if(i < old_limit) i += step;
			}
			for(; i < limit; i += step) { // Вычеркиваем кратные max_prime
				bitmap[i >> 6] |= (1 << ((i >> 1) & 31));
			}
			// ищем следущее простое
			for(unsigned int i = max_prime + 2; i < limit; i += 2) {
				if((bitmap[i >> 6] & (1 << ((i >> 1) & 31))) == 0) {
					max_prime = i;
					break;
				}
			}
		}
	}
	return ret;
}

		unsigned int max = next_prime(4000000000);
		unsigned int xe = 2;
		while(xe < 4000000000) {
			if(max % xe == 0) {printf("\nerror max = %u xe=%u\n", max, xe); break;}
			xe = next_prime(xe);
		}

typedef unsigned long long uint64_t;

uint64_t next_prime_64(uint64_t x)
{
	static unsigned int* bitmap = NULL; // биткарта
	static unsigned int bm_size = 0; // размер биткарты в байтах
	static uint64_t limit = 0; // до скольки заполнять решето

	uint64_t ret = 0;
	while(true) {
		if(x < 3) {
			if(x == 0) { // освобождение памяти
				limit = 0;
				bm_size = 0;
				if(bitmap) {
					free(bitmap);
					bitmap = NULL;
				}
			}
			ret = (x < 2) ? 2 : 3;
			break;
		}
		// поиск следующего после x
		for(uint64_t i = (x | 1) + 2; i < limit; i += 2) {
			if((bitmap[i >> 6] & (1 << ((i >> 1) & 31))) == 0) {
				ret = i;
				break;
			}
		}
		if(ret != 0) break;
		// нехватило биткарты, дорасчет
		
		uint64_t old_limit = limit;
		limit += 1048576; // 64Kb
		if(limit/16 >= bm_size) { // довыделение памяти
			unsigned int bm_size_new = bm_size / 4 * 5;
			if((unsigned int)limit/16 > bm_size_new) bm_size_new = (unsigned int)limit/16;
			if((unsigned int)(x + 1048576)/16 > bm_size_new) bm_size_new = (unsigned int)((x + 1048576)/16) + 4; 
			if(bm_size_new < bm_size) bm_size_new = 0xFFFFFFFF; // превышение разрядности
			unsigned int* p = (unsigned int*) realloc(bitmap, bm_size_new); //calloc(limit / 64 + ((limit & 63) ? 1 : 0), sizeof(unsigned int));
			if(!p) {
				printf("No memory. New size %u byte\n", bm_size_new);
				ret = 0xFFFFFFFF;
				break;
			}
			//if(bitmap != p)	printf("Move %d bytes\n", old_limit/16);
			bitmap = p;
			memset(bitmap + bm_size/4, 0, bm_size_new - bm_size);
			bm_size = bm_size_new;
		}
		uint64_t max_prime = 3; // максимальное простое
		while(true) {
			uint64_t i = max_prime * max_prime; 
			if(i >= limit) break; // дальше заполнять не надо
			uint64_t step = max_prime << 1;
			if(i < old_limit) { // выравнивание под max_prime*max_prime+2N*max_prime
				i = old_limit - (old_limit - i) % step;
				if(i < old_limit) i += step;
			}
			for(; i < limit; i += step) { // Вычеркиваем кратные max_prime
				bitmap[i >> 6] |= (1 << ((i >> 1) & 31));
			}
			// ищем следущее простое
			for(uint64_t i = max_prime + 2; i < limit; i += 2) {
				if((bitmap[i >> 6] & (1 << ((i >> 1) & 31))) == 0) {
					max_prime = i;
					break;
				}
			}
		}
	}
	return ret;
}

		limit += 1048576; // 64Kb

		if(limit < old_limit) { // превышение разрядности
			printf("No more primes :(\n");
			break;
		}

	next_prime(0); // очистка памяти
	unsigned int x = next_prime(1); // выделение памяти
	clock_t t = clock();
	x = next_prime(10000000);
	printf("10^7=%u time %d msec\n", x, clock() - t);

	next_prime(0);
	x = next_prime(1);
	t = clock();
	x = next_prime(100000000);
	printf("10^8=%u time %d msec\n", x, clock() - t);

	next_prime(0);
	x = next_prime(1);
	t = clock();
	x = next_prime(1000000000);
	printf("10^9=%u time %d msec\n", x, clock() - t);

	next_prime(0);
	x = next_prime(1);
	t = clock();
	while(x < 1000000000) x = next_prime(x);
	printf("all to 10^9=%u time %d msec\n", x, clock() - t);

		unsigned int x = next_prime(4000000000);
		while(x) {
			unsigned int y = 2;
			while(y < 65536) {
				if(x % y == 0) printf("error %u %u\n", x, y);
				y = next_prime(y);
			}
			x = next_prime(x);
		}