#include <iostream>
#include <math.h>
#include <NTL/ZZ.h>
#include <NTL/mat_ZZ.h>
using namespace NTL;
using namespace std;

void fermat(ZZ n, ZZ p)
{
	ZZ m, d1, d2, A, B, b, x;
	cout << "p/n ";
	if (n%p == 0)
		cout << "division succeeded. Divider is " << p;
	else
	{
		cout << "division failed. Factorization: " << endl;
		m = SqrRoot(n^2);
		for (x = 0; x < m; x++)
		{
			A = m + x;
			b = (A*A) - n;
			B = SqrRoot(b^2);
			d1 = A + B;
			d2 = A - B;
			/*cout << "m = " << m << endl;
			cout << "A = " << A << endl;
			cout << "b = " << b << endl;
			cout << "B = " << B << endl;
			cout << "d1 = " << d1 << endl;
			cout << "d2 = " << d2 << endl;
			cout << "d1*d2 = " << d1*d2 << endl;*/
			if (d1*d2 == n)
				cout << d1 << " and " << d2 << endl;
		}
	}
}

int main()
{
	ZZ n, p;
	cout << "n = ";
	cin >> n;
	cout << "(Must be prime!) p = ";
	cin >> p;
	fermat(n, p);
	cout << endl;
	system("pause");
	return 0;
}