--
-- Test pierwszosci - test Millera-Rabina
--
-- www.algorytm.org
-- (c)2007 Tomasz Lubinski
--

with Text_IO;
use Text_IO;

with Interfaces;
use Interfaces;

with Ada.Numerics.Discrete_Random;

procedure Miller_Rabin is
	
   powerOf2 : array (0..30) of Integer :=
   ( 1,  2,  4,  8,  16,  32,  64,
     128,  256,  512,  1024, 2048, 4096, 8192,
     16384, 32768, 65536, 131072, 262144, 524288, 1048576,
     2097152, 4194304, 8388608, 16777216, 33554432, 67108864, 134217728,
     268435456, 536870912, 1073741824 );   
       
   -- calculates a^b mod n
   function power_modulo_fast(a: Integer; b: Unsigned_32; m: Integer) return Integer is
      i: Unsigned_32;
      result: Integer := 1;
	 x: Integer := a mod m;
   begin
      i := 1;
      while i<=b loop
         x := x mod m;
         if (b and i) /= 0 then
            result := result * x;
            result := result mod m;
         end if;
         x := x*x;
         i := i * 2;
      end loop;
      return (result);
   end;
   
   -- Miller-Rabin test
   -- retrun: true - probably prime
   --         false - not prime
   function MillerRabinTest(n: Positive; k: Positive) return Boolean is
      subtype Rand_Range is Positive range 2 .. n-1;
      package Random_Gen is new Ada.Numerics.Discrete_Random (Rand_Range);
      use Random_Gen;
      a, d: Positive;
      g : Generator;
      s: Natural := 0;
      s2: Natural := 1;      
      prime: Boolean;
   begin
   
      if n < 4 then
         return TRUE;
      end if;

      if n mod 2 = 0 then
         return FALSE;
      end if;

      Reset(g);

      -- calculate s and d
      while (Unsigned_32(s2) and Unsigned_32(n-1)) = 0 loop
         s  := s + 1;
         s2 := s2 * 2;
      end loop;
      d := n/s2;

      -- try k times
      for i in 1 .. k loop
         a := Random(g);
         if power_modulo_fast(a, Unsigned_32(d), n) /= 1 then
            prime := FALSE;
            for r in 0 .. s-1 loop
               if power_modulo_fast(a, Unsigned_32(powerOf2(r)*d), n) = n - 1 then 
                  prime := TRUE;
                  exit;
               end if;
            end loop;
            if prime = FALSE then
               return FALSE;
            end if;
         end if;
      end loop;

      return TRUE;
   end;

begin

   -- check 13
   if MillerRabinTest(13, 3) then
      Put_Line("13 jest pradopodobnie liczba pierwsza");
   else
      Put_Line("13 nie jest liczba pierwsza");
   end if;

   -- check 27
   if MillerRabinTest(27, 3) then
      Put_Line("27 jest pradopodobnie liczba pierwsza");
   else
      Put_Line("27 nie jest liczba pierwsza");
   end if;

end;