#include <stdlib.h>
#include <string.h>

// returns list of primes up to limit; count stored in first element
extern unsigned long* primes(unsigned long limit)
{
	unsigned long i = 2; // start sieving with 2
	unsigned long pos, x;
	char* numbers = (char*) malloc(limit*sizeof(char));
	unsigned long* primes;
	unsigned long primecount = 0;
	memset(numbers, 1, limit);
	numbers[0] = 0;

	while(i*i <= limit)
	{
		// sieve current prime
		x = i*i;
		while(x <= limit)
		{
			numbers[x-1] = 0;
			x += i;
		}
		
		// search next prime
		pos = i;
		while(++pos <= limit)
		{
			if(numbers[pos-1] == 1)
				break;
		}
		i = pos;
	}

	// count primes
	pos = 0;
	while(++pos <= limit)
	{
		if(numbers[pos-1] == 1)
			primecount++;
	}

	// create list of primes
	primes = (unsigned long*) malloc((primecount+1)*sizeof(unsigned long)); // store count at first element
	primes[0] = primecount;
	pos = 0;
	x = 1;
	while(++pos <= limit)
	{
		if(numbers[pos-1] == 1)
			primes[x++] = pos;
	}

	free(numbers);
	return primes;
}

extern int isprime(unsigned long number, const unsigned long* primes)
{
	unsigned long count = primes[0];
	unsigned long pos = count/2 + 1;
	unsigned long dist = pos/2 + 1;
	unsigned long minpos = 0, maxpos = count;

	if(number > primes[count] || number < 2)
		return 0;

	// some kind of binary search...
	while(1)
	{
		if(number > primes[pos])
		{
			minpos = pos;
			pos += dist;
		}
		else if(number < primes[pos])
		{
			maxpos = pos;
			pos -= dist;
		}
		else
			return 1;

		dist /= 2;

		if(dist == 1)
		{
			while(++minpos <= maxpos)
				if(primes[minpos] == number)
					return 1;
			break;
		}
	}

	return 0;
}

extern unsigned long phi(unsigned long n, unsigned long* primelist)
{
	double product = (double) n;
	unsigned long i;

	//if(isprime(n, primelist))
	//	return n-1;

	for(i=1; i<=primelist[0]; i++)
	{
		unsigned long p = primelist[i];
		if(p > n)
			break;
		if(n % p != 0)
			continue;
		product *= (1 - 1.0/p);
	}
	return (unsigned long) product;
}

extern int ispermutation(unsigned long a, unsigned long b)
{
	char digits1[10];
	char digits2[10];
	memset(digits1, 0, 10);
	memset(digits2, 0, 10);

	while(a > 0 && b > 0)
	{
		char d1 = a % 10;
		char d2 = b % 10;
		digits1[d1]++;
		digits2[d2]++;
		a /= 10;
		b /= 10;
	}

	if(a > 0 || b > 0)
		return 0;

	if(memcmp(digits1, digits2, 10) == 0)
		return 1;

	return 0;
}

/*int main()
{
	long* p = primes(10000);
	//printf("%i\n", isprime(9973, p));
	int i;
	for(i=0; i<p[0]; i++)
		printf("%d: %d\n", p[i+1], isprime(p[i+1], p));
	free(p);
	return 0;
}*/