Факторизация с использованием первого (неоптимизированного) метода Монте-Карло (1-й p-алгоритм Полларда)

Разложить число на простые множители.

Неформальное описание алгоритма факторизации
1. Вызвать рекурсивную факторизацию для числа.
2. Выход из рекурсии - переданное в факторизацию число оказалось простым. Это число заносится в массив простых чисел.
3. В факторизации вызвать метод Монте-Карло для получения случайного множителя.
4. В факторизации, вызванной в пункте 1, вызвать факторизацию для полученного множителя и для исходного числа, разделенного на этот множитель.

Теоретическое обоснование алгоритма Монте-Карло (получение случайного нетривиального множителя)
Алгоритм позволяет получить нетривиальный множитель. Время работы в худшем случае - корень квадратный из исходного числа m.

Рассмотрим отображение f кольца Z(m) в кольцо Z(m):
f: Z(m) --> Z(M)
f(x) = x ^ 2 + 1(mod m)

Случайным образом выбирается элемент b0 из кольца Z(m).
Рассмотрим бесконечную последовательность элементов кольца Z(m):
b0, b1 = f(b0), b2 = f(b1), b3 = f(b2), ...

Эта последовательность представляет собой орбиту элемента b0 при отображении f. Так как все элементы последовательности принадлежат конечному множеству, то эта последовательность циклическая. Другими словами, она содержит начальный апериодичный отрезок с бесконечным периодом повторения.

Пусть p = делитель m. Рассмотрим элементы последовательности b0, b1 = f(b0), b2 = f(b1), b3 = f(b2), ... по модулю p (образы элементов b_i при каноническом полиморфизме Z(m)-->Z(p)):
с0 = b0(mod p), c1 = b1(mod p), c2 = b2(mod p), ...

Так как в Z(p) элементов меньше, чем в Z(m) (так как делитель p числа m меньше этого числа, то и замкнутое относительно теоретико-множественных операций кольцо будет содержать меньше элементов), то с большой вероятностью период последовательности с0 = b0(mod p), c1 = b1(mod p), c2 = b2(mod p), ... меньше периода последовательности b0, b1 = f(b0), b2 = f(b1), b3 = f(b2), ... .

Следовательно, существует такая пара индексов i и j, что c_i = b_j при условии b_i != b_j.

Это означает, что b_i == b_j(mod p), b_i != b_j(mod m).

Таким образом, делаем вывод: (b_i - b_j) делится на p, но не делится на m. Следовательно, НОД((b_i - b_j), m) нетривиален, и нам удалось разложить исходное число m на множители.

Указанный выше процесс формально описывается следующим образом - 1-й алгоритм Монте-Карло сводится к поиску цикла в бесконечной рекурсивной последовательности, состоящей из элементов конечного множества. При этом вместо сравнения между собой двух элементов вычисляется НОД разности этих элементов и исходного числа m. Алгоритм завершается, когда НОД нетривиален.

Пошаговое описание 1-го метода Монте-Карло
Пусть надо разложить на множители число P
Шаг 1. Выбирается многочлен f(x) с целочисленными коэффициентами, степени не выше 2. Обычно берется многочлен вида y = x ^ 2 + c(mod P).
Шаг 2. Случайно выбирается x0 = y0 меньше P.
Шаг 2. Вычисляются значения x_i = f(x_i − 1)(modP),y_i = f(f(y_i − 1))(modP).
Шаг 3. Находится d = ( | x_i − y_i |, P).
Шаг 4. Если d = 1, происходит переход на шаг 2, если d = P, происходит остановка - разложение на множители провести не удалось. Если 1 < d < P, то найдено разложение числа P.

P. S. Важно понимать что указанный выше метод генерирует случайный множитель числа, а не простой (как зачастую ошибочно считают). Хотя при достаточном количестве проверок алгоритма на контрольном числе простые множители имеют место быть.

Реализация алгоритма:

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Разложение числа на простые множители                                  ++++ 
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
#include <ctime>
#include <iostream>
 
 
//----------------------------------------------------------------------1----
using namespace std;
 
#define SIZE 1000
//---------------------------------------------------------------------------
 
 
//----------------------------------------------------------------------2----
//Простые множители.
int prime_factors[SIZE];                                   
 
//Их количество.
int count = 0;                                             
//---------------------------------------------------------------------------
 
 
//----------------------------------------------------------------------3----
//НОД.
int GreatestCommonDivisor(int arg1, int arg2);             
 
//Число по модулю.
int Mod(int operand, int module);                         
 
//Генератор случайных чисел.
int RndGen(int num);                                       
 
//Получение случайного множителя num.
int GenerateFactor(int num);                           
 
//Разложение на простые множители.
void Factorization(int num, int prime_factors[]);         
 
//Проверка числа на простоту.
bool IsPrime(int num);                                                                                               
 
//Вывод результатов.
void PrintResult(int num, int prime_factors[], int count);  
 
//---------------------------------------------------------------------------
 
 
//----------------------------------------------------------------------4----
int main()
{
	int num;
 
	cout << "Enter number: ";
	cin >> num;
	cout << endl;
 
	//Факторизация.
	Factorization(num, prime_factors);					   
 
	//Вывод результатов.
	PrintResult(num, prime_factors, count);	                
 
	system("PAUSE");
 
	return 0;
}
//---------------------------------------------------------------------------
 
 
//----------------------------------------------------------------------5----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Получение НОД arg1 и arg2 бинарным алгоритмом Евклида.                 ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
int GreatestCommonDivisor(int arg1, int arg2)
{
	int shift;
	int difference;
 
	//НОД(0, n) = n; НОД(m, 0) = m.
	if (arg1 == 0 || arg2 == 0)
		return arg1 | arg2;
 
	/*Если m, n чётные, тогда НОД(m, n) = 2 * НОД(m / 2, n / 2).
	Пусть shift := lg K, где K наибольшая степень 2 такая, что оба числа 
	целиком делятся на 2 в степени K. */
	for (shift = 0; ((arg1 | arg2) & 1) == 0; ++shift) 
	{
		arg1 >>= 1;
		arg2 >>= 1;
	}
 
	// Если m чётное, тогда НОД(m, n) = НОД(m / 2, n).
	while ((arg1 & 1) == 0)
		arg1 >>= 1;
 
	//Далее считается, что arg1 нечетное. 
	do 
	{
		//Если n чётное, тогда НОД(m, n) = НОД(m, n / 2).		
		while ((arg2 & 1) == 0)  
			arg2 >>= 1;
 
		/*Далее считается, что arg1 и arg2 нечетные. Следовательно, 
		arg1 - arg2 четно.
		Пусть arg1 = min(arg1, arg2), arg2 = (arg1 - arg2) / 2.
		Если m, n нечётные и m < n, тогда НОД(m, n) = НОД(n - m, m).*/
		if (arg1 < arg2) 
		{
			arg2 -= arg1;
		} 
		//Если m, n нечётные и m > n, тогда НОД(m, n) = НОД(m - n, m).        
		else 
		{
			difference = arg1 - arg2;
			arg1 = arg2;
			arg2 = difference;
		}
 
		arg2 >>= 1;
	} 
	while (arg2 != 0);
 
	return arg1 << shift;
}
 
//----------------------------------------------------------------------6----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Отображение множества целых числе в кольцо [0; module]                 ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
int Mod(int operand, int module)
{
	while (operand < 0) 
	{
		operand += module;
	}
 
	operand %= module;
 
	return operand;
}
//---------------------------------------------------------------------------
 
//----------------------------------------------------------------------7----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Генерация случайного числа в интервале [0; num]. Отображение в       ++++
//интервал осуществляется с помощью функции Мod(int num).                ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
int RndGen(int num)
{
	srand((unsigned)time(NULL));
 
	int result;
 
	result = rand();
	result = Mod(result, num);
 
	return result;
}
//---------------------------------------------------------------------------
 
//----------------------------------------------------------------------8----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Факторизация числа. Все операции                                       ++++
//являются операциями в кольце по модулю num.                            ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void Factorization(int num, int prime_factors[])
{															
	//Выход из рекурсии.
	if (IsPrime(num))     
	{   
		//Проверка числа на простоту. Если простое - то записать его.
		prime_factors[count++] = num;					   
 
		return;
	}
 
	//Текущий множитель.
	int current_multiplier;								   
 
	//В бесконечном цикле генерируем случайные множители числа.
	while (true)                                          
	{                                                   
		current_multiplier = GenerateFactor(num);          
 
		/*Получив отличный от числа множитель, выходим из цикла, чтобы 
		рекурсивно вызвать факторизацию для полученного множителя и 
		для "оставшейся" части числа.*/
		if (current_multiplier != num)					 
			break;                                      
	}													 
 
	Factorization(num / current_multiplier, prime_factors);
 
	Factorization(current_multiplier, prime_factors);
}
//---------------------------------------------------------------------------
 
//----------------------------------------------------------------------9----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Получения случайного множителя первым методом Монте-Карло.             ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
int GenerateFactor(int num)
{
	//Глубина поиска случайного множителя.
	int max_steps			= 100000;					   
 
	//Текущая глубина.
	int current_step		= 0;                      
	//Постоянная c функции отображения f(n)= n ^ 2 + c.
	int polynomial_const	= 2;						   
 
	//Случайноe число в интервале [0; num].
 
	int x = RndGen(num);								  
 
	//Применим отображение f(n)= n ^ 2 + c.
	int y = Mod(x * x + polynomial_const, num);			  
 
	/*Пока НОД тривиален и глубина поиска не исчерпана применить отображение
	один раз к x и два раза к y.*/
	int divisor = GreatestCommonDivisor(abs(x - y), num);  
 
	while (current_step < max_steps && divisor == 1)	   
	{													  
		x = Mod(x * x + polynomial_const, num);			  
		y = Mod(y * y + polynomial_const, num);			   
		y = Mod(y * y + polynomial_const, num);			   
 
		divisor = GreatestCommonDivisor(abs(x - y), num);  
 
		current_step++;
	}
	//Успех == НОД нетривиален.
	return divisor;										   					
}
//---------------------------------------------------------------------------
 
//----------------------------------------------------------------------10----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Проверка числа на простоту "классическим" переборным алгоритмом.       ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
bool IsPrime(int num)									   
{	
	/*Брать числа, квадраты которых меньше исследуемого num, и делить
	на них num. Если num не делится нацело ни на одно из них, то оно
	простое.*/
	for (int i = 2; i * i <= num; i++)                    
		if (num % i == 0)								   
			return false;								 
 
	return true;
}
//---------------------------------------------------------------------------
 
//----------------------------------------------------------------------11----
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Вывод полученных простых множителей.                                   ++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void PrintResult(int num, int prime_factors[], int count)
{
	if (count > 1)
		cout << num << " = " << prime_factors[0] << " * " << prime_factors[1];
	else 
		cout << num << " = " << prime_factors[0];
 
	for (int j = 2; j < count; j++)
		cout << " * " << prime_factors[j];
 
	cout << endl << endl;
}
//---------------------------------------------------------------------------

Оценка сложности
Сложность реализованного выше алгоритма факторизации с использованием получения случайного множителя первым методом Монте-Карло: O(n ^ (1/2)) + k * O(m_i ^ (1/2)), где n - исходное число, k - количество генерируемых методом Монте-Карло множителей, m_i - эти множители. Первое слагаемое отвечает за сложность первого метода Монте-Карло в худшем случае, второе - за суммарную сложность проверок k получаемых случайных множителей на простоту "классическим" переборным алгоритмом.

Критическое замечание
Узким местом указанной выше реализации является проверка числа на простоту IsPrime(). В приложениях реального времени (взлом шифров с большими ключами шифрования, к примеру) необходимо реализовать более эффективные методы (к примеру, вероятностный тест простоты Миллера-Рабина).

Ключевые слова: 
факторизация, первый метод Монте-Карло, 1-й p-алгоритм Полларда, разложение на простые множители
ВложениеРазмер
www.opita.net_monte-karlo_1.rar6.86 кб