Filosofy3 / Filosofy3 / Filosofy3

// Философы.cpp: определяет точку входа для консольного приложения.
//

#include "stdafx.h"
#include <windows.h>
#include <stdio.h>
#include <ctime>
#include <iostream>

using namespace std;

//Класс кривой
class curve
{
public:
	int a;
	int b;
	curve (int, int);
};

//Класс точки
class point
{
public:
	int x;
	int y;
	point (int, int);
};

//Конструктор кривой
curve :: curve (int a1, int b1)
{
	a = a1;
	b = b1;
}

//Конструктор точки
point :: point (int x1, int y1)
{
	x = x1;
	y = y1;
}

//Н.О.Д.
int NOD (int a1, int a2)
{
      int NODn = (a1 < a2 ? a1 : a2);
      for (int mNOD = 2; mNOD <= NODn; mNOD++)
      {
           if (a1 % mNOD == 0 && a2 % mNOD == 0)
           {
				return mNOD;
           }
      }
      return 1;
}

//Возведение в степень
long unsigned VozvedenieVStepModP (int multiplier, int step, int mod)
        {
            multiplier = multiplier % mod;
            long ans = 1;
           
            while (step > 0)
            {
                if (step % 2 == 0)
                {
                    step /= 2;
                    multiplier = multiplier * multiplier % mod;
                }
                else
                {
                    step--;
                    ans = ans * multiplier % mod;
                }
            }

            return ans;
        }

//Тест простоты Миллера-Рабина
bool Miller_Rabin (int n)
        {
			if (n == 3 || n == 5 || n == 7 || n == 11)
				return true;
            if (n > 2)
            {
                if (n % 2 == 1)
                {
                    int n1 = n - 1;
                    int s = 0, d = 1;
                    while (n1 % 2 == 0)
                    {
                        n1 = n1 / 2;
                        s++;
                    }
                    d = n1;
                    bool f = true;


                    for (int a = 2; a <= 11; a++)
                    {
                        if (f)
                        {
                            f = false;
                            int *x = new int[s];

                            x[0] = VozvedenieVStepModP (a, d, n);

                            if (x[0] == 1 || x[0] == n - 1)
                            {
                                f = true;
                                continue;
                            }
                            else
                            {
                                for (int i = 1; i <= s - 1; i++)
                                {
                                    x[i] = VozvedenieVStepModP(x[i - 1], 2, n);
                                    if (x[i] == (n - 1))
                                    {
                                        f = true;
                                        break;
                                    }
                                }
                            }
                        }
                        else
                            return f;
                    }
                    return f;
                }
                return false;
            }
            return true;
        }

//Функция суммы двух точек на кривой
point sum_of_points (point point1, point point2, curve c)
{
	int lambda = 0;
	if (point1.x != point2.x || point1.y != point2.y)
		lambda = (point2.y - point1.y) / (point2.x - point1.x);
	else
		lambda = (3 * point1.x * point1.x + c.a ) / 2 * point1.y;
	int x3 = lambda * lambda - point1.x - point2.x;
	int y3 = lambda * (point1.x - x3) - point1.y;
	return point (x3, y3);
}

//Вычисление кратной точки
point point_multiplication (point point1, int k, curve c)
{
	point result_point = point1;
	for (int i = 0; i < k - 1; i++)
		result_point = sum_of_points (result_point, point1, c);
	return result_point;
}

//Отыскание наибольшей степени r такой, что p^r < B1
int find_degree (int p, int B1)
{
	int r = 1;
	int p1 = p;
	while (p1 < B1)
	{
		p1 = p1 * p;
		if (p1 < B1)
			r++;
	}
	return r;
}

//Алгоритм Ленстры
long Lenstra (long n)
{
	int B1 = 10000;
	srand(time(NULL));
	STEP2:
	int a = rand() % n;
	cout<<"a = "<<a<<endl;
	int x = rand() % n;
	cout<<"x = "<<x<<endl;
	int y = rand() % n;
	cout<<"y = "<<y<<endl;
	int b = (y * y - x * x * x - a * x) % n;
	if (b < 0)
		b = n + b;
	cout<<"b = "<<b<<endl;
	int g = NOD (n, 4 * a * a * a + 27 * b * b);
	cout<<"g = "<<g<<endl;
	curve c = curve (a,b);
	cout<<"Elliptic curve: y^2 = x^3 + "<<c.a<<"x + "<<c.b<<endl;
	cout<<"Point P("<<x<<", "<<y<<")"<<endl;
	/*if (g == n) 
		goto STEP2;*/
	if (1 < g && g < n)
		return g;
	point P0 = point (x,y);
	point P = P0;
	int r = 0;
	for (int i = 2; i < B1; i++)
	{
		if (Miller_Rabin (i))
		{
			r = find_degree (i, B1);				
			for (int j = 0; j < r; j++)
			{
				P = point_multiplication (P, j, c);
			}
			if (NOD(n, P.x) > 1 || NOD(n, P.y) > 1)
				return i;
		}	
	}

	//Вторая стадия алгоритма
	int B2 = B1 * 2;
	for (int i = B1; i < B2; i++)
	{
		if (Miller_Rabin (i))
		{
			P = point_multiplication (P, i, c);
			if (NOD(n, P.x) > 1 || NOD(n, P.y) > 1)
				return i;
		}
	}

	return 1;
}

int _tmain(int argc, _TCHAR* argv[])
{
	int n, b, k;
	/*scanf("%ld", &g);*/
	cout<<"Hello!"<<endl;
	cout<<"Enter n: ";
	cin>>n;
	/*cout<<"Enter b:"<<endl;
	cin>>b;
	curve c = curve(a,b);
	cout<<"Elliptic curve: y^2 = x^3 + "<<c.a<<"x + "<<c.b<<endl;*/
	long p = Lenstra (n);
	if (p)
		cout<<"The answer is "<<p<<", "<<n<<" = "<<p<<" * "<<n/p<<endl;
	return 0;

}