Test pierwszości liczby metodą probabilistyczną (Solovaya - Strassena)

Poniższy kod służy do testowania pierwszości liczby metodą probabilistyczną Solovaya - Strassena. Krótki opis: algorytm z prawdopodobieństwem 1/2 określa prawidłowo, czy liczba jest pierwsza, czy złożona. Należy przeprowadzić k testów, żeby z prawdopodobieństwem błędu 2^-k określić, że liczba jest pierwsza. W praktyce oznacza to wywołanie statycznej metody IsPrime z drugim parametrem o wartości równej ilości powtórzeń testu. 10 oznacza prawdopodobieństwo pomyłki równe 2^-10 = 0,90234375%, bezpieczną wartością parametru jest 64. Teoretyczna złożoność obliczeniowa algorytmu to O(k × log³n).
Kod jest tylko demonstracją sposobu działania algorytmu i działa jedynie na typie long (64b). Jeśli potrzebujesz wyższego zakresu liczb (np. do zaimplementowania RSA) użyj jednej z gotowych bibliotek, np. BigNumber.

public class SolovayStrassenPrimeTest
{
	private static Random r = new Random((int)DateTime.Now.Ticks);

	private static int L(long a, long p)
	{
		if (p < 2)
			throw new ArgumentOutOfRangeException("p", "p must not be < 2");
		if (a == 0)
			return 0;
		if (a == 1)
			return 1;

		int result;
		if ((a & 1) == 0)
		{
			result = L(a >> 1, p);
			if (((p * p - 1) & 8) != 0)
				result = -result;
		}
		else
		{
			result = L(p % a, a);
			if (((a - 1) * (p - 1) & 4) != 0)
				result = -result;
		}
		return result;
	}


	private static long RandomLong()
	{
		return r.Next() >> 32 | r.Next();
	}

	
	private static long GCD(long a, long b)
	{
		long Remainder;

		while (b != 0)
		{
			Remainder = a % b;
			a = b;
			b = Remainder;
		}

		return a;
	}

	
	private static long ModExp(long x, long e, long m)
	{
		long result = 1;
		x = x % m;
		while (e > 0)
		{
			if ((e & 1) == 1)
				result = (result * x) % m;
			e = e >> 1;
			x = (x * x) % m;
		}

		return result;
	}


	public static bool IsPrimeBruteForce(long n)
	{
		for (int i = 2; i <= (int)Math.Sqrt(n); i++)
			if (n % i == 0)
				return false;

		return true;
	}


	public static bool IsPrime(long n, int k)
	{
		if (n <= 2)
			throw new ArgumentOutOfRangeException();
		if ((n & 1) == 0)
			return false;

		for (int i = 0; i < k; i++)
		{
			long a = RandomLong() % (n - 2) + 2;
			long x = L(a, n);
			if (x == 0 || ModExp(a, (n - 1) / 2, n) != (x + n) % n)
				return false;
		}
		return true;
	}
}