-- file FermatBest.Ada

-- primes by Fermat's method. Best version

-- 08/05/24


with Ada.Numerics.Big_Numbers.Big_Integers;
package Big renames Ada.Numerics.Big_Numbers.Big_Integers;


with Big; use Big;
package Numbers is
	subtype Number is Big_Integer;
end Numbers;


with Ada.Text_IO; use Ada.Text_IO;
with Big; use Big;
procedure Put_Num(N: Big_Integer) is
begin
	Put(To_String(N, 0));
end Put_Num;


with Big; use Big;
with Numbers; use Numbers;
function Mersenne(P: Integer) return Number is
begin
	return To_Big_Integer(2)**P - 1;
end Mersenne;


with Ada.Text_IO; use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Big; use Big;
with Numbers; use Numbers;
with Mersenne;
with Put_Num;
procedure Get_Num(N: out Number) is
	C: Character;
	P: Integer;
	EOL: Boolean;
begin
	Look_Ahead(C, EOL);
	loop
		exit when not EOL; 
		Skip_Line;
		Look_Ahead(C, EOL); 
	end loop;
	if C = 'M' or C = 'm' then
		Get(C); Get(P);
		N := Mersenne(P);
		Put("  ="); Put_Num(N); New_Line;
	else
		Get(P);  N := To_Big_Integer(P);   -- sort of Get(N);
	end if;
end Get_Num;


package Primes 
	with Elaborate_Body is
	Max: constant := 20000;
	Prime: array (1..Max) of Integer;
end;


package body Primes is
	N: Integer := 2;
	Index: Integer := Prime'First;
	Found: Boolean := False;
begin
	-- initialization part to build prime table
	Prime(Prime'First) := 2;
	loop
		N := N+1; Found := True;
		for I in Prime'First .. Index loop
			if N rem Prime(I) = 0 then -- divides by existing prime
				Found := False; exit;
			end if;
		end loop;
		if Found then

			-- found a new prime
			Index := Index+1; Prime(Index) := N;
			exit when Index = Max;
		end if;
	end loop;
end Primes;


package Squares
	with Elaborate_Body is
	Square_Digits: Integer := 4;
	Square_Ten_Power: Integer := 10**Square_Digits;
	Poss_Last_Digits: array(0..Square_Ten_Power-1) of Boolean := 
       (others => False);
end;


package body Squares is
begin
	-- make Poss_Last_Digits array
	for I in 1 .. Square_Ten_Power loop
		Poss_Last_Digits(I*I mod Square_Ten_Power) := True;
	end loop;
end Squares;		


with Big; use Big;
with Numbers; use Numbers;
procedure Square_Root(XX: in Number; Try: Number;
        X: out Number; R: out Number) is
	-- X is the largest X such that X*X + R equals XX with R >= zero
	-- Try is initial guess 
	K, KK: Number;
	DK: Number;
	Three: constant Number := 3;
begin
	K := Try;
	loop
		KK := K*K;
		DK := (XX - KK)/(K+K);

		-- iterate until nearly there
		if abs DK < Three then -- nearly there
			if XX < KK then
				K := XX/K;		-- ensure K*K is less than XX
			end if; 
			exit;
		end if;
		-- do another iterate
		K := K+DK;
	end loop;

	-- now loop from below

	loop
		KK := K*K;
		if KK >= XX then 
			if KK = XX then
				X := K; R := 0; return;
			end if;
			X := K - 1; R := XX - X*X; return;
		end if;
		K := K + 1;
	end loop;
end Square_Root;


with Square_Root;
with Big; use Big;
with Numbers; use Numbers;
with Squares; use Squares;
procedure Fermat_Best(N: in Number; Min_Prime: in Number;
        P, Q: out Number) is
	-- we know that factors up to Min_Prime have been removed,
	-- so max square to try is
	-- roughly (N/Min_Prime + Min_Prime)/2; we add 2 for luck
	Two: constant Number := 2;
	Number_Square_Ten_Power: constant Number := 
						To_Big_Integer(Square_Ten_Power);
	X: Number;
	Y: Number;
	R: Number;
	K: Number;
	DK: Number;
	N2, X2, K2: Integer;
	Last_Digits: Integer;
	Try: Number := 1;
	Max_Square : constant Number := 
						(N/Min_Prime + Min_Prime)/ Two + Two;
begin
	Square_Root(N, 1, X, R);
	if R = 0 then
		-- N was a perfect square
		P := X; Q := X;
		return;
	end if;

		K := X*X-N;
		DK := X+X+1;
		N2 := To_Integer(N rem Number_Square_Ten_Power);
		X2 := To_Integer(X rem Number_Square_Ten_Power);

	loop
		X := X + 1;
		if X > Max_Square then
			-- must be prime
			P := N; Q := 1;
			return;
		end if;

		X2 := (X2 + 1) rem Square_Ten_Power;
		K := K + DK;
		DK := DK + Two;

		K2 := (X2*X2-N2) mod Square_Ten_Power;
		Last_Digits := K2;
		if Poss_Last_Digits(Last_Digits) then

			Square_Root(K, Try, Y, R);
			if R = 0 then
				-- X*X-N was a perfect square
				P := X+Y; Q := X-Y;
				return;
			end if;
			Try := Y;
		end if;
	end loop;
exception
	when others => P := 0; Q := 0;
end Fermat_Best;


with Numbers; use Numbers;
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Get_Num;
with Put_Num;
with Big; use Big;
with Primes; use Primes;
with Fermat_Best;
with Ada.Calendar; use Ada.Calendar;
procedure Do_Fermat_Best is
	NN, PP, QQ: Number;
	T_Start, T_End: Time;
	package Duration_IO is new Fixed_IO(Duration);
	use Duration_IO;
begin	
	Put_Line("Welcome to Fermat's best method");
	loop
<<Again>>
		begin	
			New_Line;
			Put("Insert number N = "); Get_Num(NN);
		exception
			when others =>
				Put_Line("Not a number or too big");  Skip_Line;
				goto Again;
		end;

		exit when NN = 0;

			-- check to see if N divides by a known prime

		for I in Prime'range loop
			loop
				if NN rem To_Big_Integer(Prime(I)) = 0 then
					Put("N divides by "); Put(Prime(I), 0);
					Put_Line(" so removing factor");
					NN := NN / To_Big_Integer(Prime(I));
				else
					exit;
				end if;
			end loop;
		end loop;

		if NN = 1 then
			Put_Line("all factors removed");
		else

			T_Start := Clock;
			Fermat_Best(NN, To_Big_Integer(Prime(Primes.Max)),
								 PP, QQ);
			T_End := Clock;
			if PP = 0 and QQ = 0 then
				Put_Line("Overflowed internally");
				goto Again;
			end if;
			Put("Two factors are   ");
			Put_Num(PP); Put("   "); Put_Num(QQ); New_Line;
			Put("Time = ");  Put(T_End-T_Start, 2, 1);  New_Line;
		end if;
	end loop;
	Put_line("Goodbye");  Skip_Line(2);
end Do_Fermat_Best;