/*
Fast prime number generation routines by M Phillips - 2008.

This code is provided as is with no warranties or guarantees of
any kind.

Please send an email to M Phillips (mbp2@i4free.co.nz)
  - if you use this file in a released product, or
  - if you find any bugs, or
  - if you have any suggestions
*/
#ifndef FAST_PRIMES_H
#define FAST_PRIMES_H

#include <cmath>
#include <vector>

template<typename T> void PrimesLessThanN(const T n, std::vector<T> &primes) {
	primes.clear();
	primes.push_back(2);
	std::vector<T> limits;
	limits.push_back(4);
	size_t rootIdx = 0;
	for (T numToTest = 3U; numToTest < n; numToTest += 2U) {
		bool prime = true;
		while (primes[rootIdx] * primes[rootIdx] < numToTest)
			++rootIdx;
		for (size_t i = 1; i <= rootIdx; ++i) {
			T &limit = limits[i];
			while (limit < numToTest)
				limit += primes[i];
			if (limit == numToTest) {
				prime = false;
				break;
			}
		}
		if (prime) {
			primes.push_back(numToTest);
			limits.push_back(numToTest*numToTest);
		}
	}
}

template<typename T> void PrimesLessThanN_Modulus(const T n, std::vector<T> &primes) {
	primes.clear();
	primes.push_back(2);
	size_t rootIdx = 0;
	for (T numToTest = 3U; numToTest < n; numToTest += 2U) {
		bool prime = true;
		while (primes[rootIdx] * primes[rootIdx] < numToTest)
			++rootIdx;
		for (size_t i = 1; i <= rootIdx; ++i) {
			if (!(numToTest % primes[i])) {
				prime = false;
				break;
			}
		}
		if (prime)
			primes.push_back(numToTest);
	}
}

template<typename T> void PrimesLessThanN_Eratosthenes(const T n, std::vector<T> &primes) {
	primes.clear();
	std::vector<char> notPrimes(n);
	std::fill(notPrimes.begin(), notPrimes.end(), false);
	T sqrtn = (T)sqrt((double)n);
	for (T numToTest = 2U; numToTest < n; ++numToTest) {
		if (!notPrimes[numToTest]) {
			primes.push_back(numToTest);
			if (numToTest <= sqrtn) {
				T pos = numToTest * numToTest;
				while (pos < n) {
					notPrimes[pos] = true;
					pos += numToTest;
				}
			}
		}
	}
}

#endif