path: root/src/crypto-types-big_numbers-mod_utils.adb

-- This program is free software; you can redistribute it and/or
-- modify it under the terms of the GNU General Public License as
-- published by the Free Software Foundation; either version 2 of the
-- License, or (at your option) any later version.

-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-- General Public License for more details.

-- You should have received a copy of the GNU General Public License
-- along with this program; if not, write to the Free Software
-- Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
-- 02111-1307, USA.

-- As a special exception, if other files instantiate generics from
-- this unit, or you link this unit with other files to produce an
-- executable, this unit does not by itself cause the resulting
-- executable to be covered by the GNU General Public License. This
-- exception does not however invalidate any other reasons why the
-- executable file might be covered by the GNU Public License.

with Crypto.Types.Random;
with Crypto.Asymmetric.Prime_Tables;
--with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Ada.Text_IO;

separate(Crypto.Types.Big_Numbers)

package body Mod_Utils is

   pragma Optimize (Time);
   use Crypto.Asymmetric.Prime_Tables;


   ---------------------------------------------------------------------------

   function Patch(Item, N : Big_Unsigned)
                 return Big_Unsigned is
      Diff : constant Big_Unsigned:=((Big_Unsigned_Last - N) + 1) mod N;
   begin
      return Add(Item,Diff,N);
   end Patch; pragma Inline(Patch);

   ---------------------------------------------------------------------------

   function Add(Left, Right, N : Big_Unsigned) return Big_Unsigned is
      L : constant Big_Unsigned := Left mod N;
      R : constant Big_Unsigned := Right mod N;
      Result : constant Big_Unsigned := L + R;
   begin
      if Result < Max(L,R) then
         return Patch(Result,N);
      else return
        Result mod N;
      end if;
   end Add;

   ---------------------------------------------------------------------------

   function Sub(Left, Right, N : Big_Unsigned) return Big_Unsigned is
      L : constant Big_Unsigned := Left  mod N;
      R : constant Big_Unsigned := Right mod N;
   begin
      if R > L then
         return N - R + L;
      else return L-R;
      end if;
   end Sub;

   ---------------------------------------------------------------------------

   function Div(Left, Right, N : Big_Unsigned) return Big_Unsigned is
   begin
      return Mult(Left,Inverse(Right,N),N);
   end Div; pragma Inline(Div);


   ---------------------------------------------------------------------------

   --from Erik-Zenners handout "Zahlentheoretische Algorithmen"
   function Pow(Base, Exponent, N : Big_Unsigned) return Big_Unsigned is
      L : constant Big_Unsigned := Base mod N;
      R : constant Big_Unsigned := Exponent;
      Result : Big_Unsigned := Big_Unsigned_One;
   begin
      if L = Big_Unsigned_Zero or L = Big_Unsigned_One then
         return L;
      elsif R = Big_Unsigned_Zero then return Big_Unsigned_One;
      else
         -- Square_And_Muliply
         for I in reverse  0..Bit_Length(R)-1 loop
            Result :=  Mult(Result,Result,N);
            if (Shift_Right(R, I) mod 2) = Big_Unsigned_One then
               Result := Mult(Result,L,N);
            end if;
         end loop;
         return Result mod N;
      end if;
   end Pow;

   ---------------------------------------------------------------------------

   --based on  Erik-Zenners handout "Zahlentheoretische Algorithmen"
   -- (ext. Euklid)
   -- This function returns Big_unsigned_Zero if X have no inverse mod n
   function Inverse(X, N : Big_Unsigned) return Big_Unsigned is
      B : Big_Unsigned := X mod N;
      A : Big_Unsigned := N;
   begin
      -- if gcd(A,B) /= 1 then A have no inverse mod B
      if B = Big_Unsigned_Zero or A = Big_Unsigned_Zero or
        Gcd(A,B) /= Big_Unsigned_One then
         return  Big_Unsigned_Zero;
      end if;

      declare
         T : Big_Unsigned := Big_Unsigned_One;
         Tstrich,  Tempt : Big_Unsigned;
         Q, R : Big_Unsigned;
      begin
         loop
            Big_Div(A,B,Q,R);
            if(R = Big_Unsigned_Zero) then
               return T;
            end if;

            A:=B;
            B:=R;

            Tempt:=T;

            T:=Sub(Tstrich,Mult(Q,T,N),N);

            Tstrich:=Tempt;
         end loop;
      end;
   end Inverse;

   ---------------------------------------------------------------------------

   function Get_Random(N : Big_Unsigned) return Big_Unsigned is
      Result : Big_Unsigned;
   begin
      Random.Read(Result.Number);

      for I in reverse 0..N.Last_Index loop
         if Result.Number(I) /= 0 then
            Result.Last_Index := I;
            exit;
         end if;
      end loop;
      return Result mod N ;
   end Get_Random;

   ---------------------------------------------------------------------------

   -- this function returns true if X is a Mersenne prim number
   function Lucas_Lehmer_Test(X : Big_Unsigned) return Boolean is
      Is_Mp : Boolean := false;
   begin

      if X.Last_Index = 0 then
         for I in 2..Word'Size-1 loop
            if  X.Number(0) = Shift_Left(2,I)-1 then
               Is_Mp := True;
               exit;
            end if;
         end loop;
         if Is_Mp = False then return False;
         end if;
      else
         for I in 0..X.Last_Index loop
            if X.Number(I) /= Word'Last then return False;
            end if;
         end loop;
      end if;

      declare
         P : constant Word := Word(Bit_Length(X)-1);
         S : Big_Unsigned := Big_Unsigned_Two+2; --S(1) = 4;
      begin
         for I in 2..P-1 loop
            S := (Mult(S,S,X) - 2) mod X;
         end loop;

         if S = Big_Unsigned_Zero then return True;
         else return False;
         end if;
      end;
   end Lucas_Lehmer_Test;

   ---------------------------------------------------------------------------

   --from Erik-Zenners handout "Zahlentheoretische Algorithmen"
   function Is_Miller_Rabin_Witness(Wit, X : Big_Unsigned) return Boolean is

      B : constant Big_Unsigned := X-1;
      Result : Big_Unsigned := Big_Unsigned_One;
      Root : Big_Unsigned;
   begin
      for I in reverse 0..Bit_Length(B)-1 loop
         Root := Result;
         Result := Mult(Result, Result, X);
         if ((Result = Big_Unsigned_One) and
             (Root /= Big_Unsigned_One and Root /= B)) then return True;
         elsif (Shift_Right(B,I) mod 2) = Big_Unsigned_One then
            Result := Mult(Result, Wit, X);
         end if;
      end loop;
      if Result /= Big_Unsigned_One then return True;
      else return False;
      end if;
   end Is_Miller_Rabin_Witness;

   ---------------------------------------------------------------------------

   -- Test if Wit is a witness for N
   -- If Wit is a wittness then N is no prime
    function Is_Simple_Witness(Wit, N : Big_Unsigned) return Boolean is
    begin
       -- is Wit a "Miller-Rabin"-witness
       if (Wit /= (N-Big_Unsigned_One)) and  (Wit /= Big_Unsigned_One) and
         Mult(Wit,Wit,N) = Big_Unsigned_One  then return True;

       elsif Gcd(Wit,N) /= Big_Unsigned_One then  return True;

       -- is Wit a "Fermat-Witness"
       -- elsif Pow(Wit,N-1,N) /= Big_Unsigned_One then  return True;
       else return False;
       end if;
    end Is_Simple_Witness;

    ---------------------------------------------------------------------------


     -- Returns true if N passes the specified number of Miller-Rabin tests.
   function Passed_Miller_Rabin_Test(X : Big_Unsigned; S : Positive)
                                    return Boolean is
      Witness : Big_Unsigned;
   begin
      -- Do the tests
      for I in 1..S loop
         -- Generate a uniform random on (1, X)
         loop
            Witness := Get_Random(X);
            exit when Witness > Big_Unsigned_One;
         end loop;
         if Is_Miller_Rabin_Witness(Witness, X) then
            return False;
         end if;
         end loop;
      return true;
   end Passed_Miller_Rabin_Test;

   ---------------------------------------------------------------------------

   function Pass_Prime_Test(X : Big_Unsigned; Status : Hardness)
                           return Boolean is
      Rounds : Natural;
      X_Bit_Size : constant Natural := Bit_Length(X);
   begin
       if X < Big_Unsigned_Two then return False;
       elsif Is_Even(X) then
          if X = Big_Unsigned_Two then return True;
          else return False;
          end if;
       end if;

       --X is odd

       for I in One_Digit_Primes'First+1..One_Digit_Primes'Last loop
         if X = Word(One_Digit_Primes(I)) then return true;
         elsif X mod Word(One_Digit_Primes(I)) = Big_Unsigned_Zero then
            return False;
         end if;
       end loop;

      for I in Two_Digit_Primes'Range loop
         if X = Word(Two_Digit_Primes(I)) then return true;
         elsif X mod Word(Two_Digit_Primes(I)) = Big_Unsigned_Zero then
            return False;
         end if;
      end loop;

      if Lucas_Lehmer_Test(X) then
         return True;
      end if;

      for I in Three_Digit_Primes'Range loop
         if X = Word(Three_Digit_Primes(I)) then return true;
         elsif X mod Word(Three_Digit_Primes(I)) = Big_Unsigned_Zero then
            return False;
         end if;
      end loop;

       -- The relationship between the certainty and the number of rounds
      -- we perform is given in the draft standard ANSI X9.80, "PRIME
      -- NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
      -- Comment:
      -- I don't have a look on this paper. =:) I borrowed this
      -- "algorithmen" from the j2sdk1.4.1 library (java/math/BigInteger.java)
      -- If you have the permission to send me the draft standard ANSI X9.80
      -- then send it, please!
      -- I'm a student. I have no money for ANSI or IEEE drafts. :-(
      -- It's right to require money to read a draft?
      -- This really really sucks! SCNR!

      if    (X_Bit_Size <  100) then Rounds := 50;
      elsif (X_Bit_Size <  256) then Rounds := 27;
      elsif (X_Bit_Size <  512) then Rounds := 15;
      elsif (X_Bit_Size <  768) then Rounds :=  8;
      elsif (X_Bit_Size < 1024) then Rounds :=  4;
      else                           Rounds :=  2;
      end if;

      declare
         Witness : Big_Unsigned;
      begin
         if Status = Weak then
            for I in 1..Rounds loop
               loop
                  Witness := Get_Random(X);
                  exit when Witness > Big_Unsigned_Two;
               end loop;
               if Is_Simple_Witness(Witness,X) then return False;
               end if;
            end loop;
         else
            for I in 1..Rounds loop
               loop
                  Witness := Get_Random(X);
                  exit when Witness > Big_Unsigned_Two;
               end loop;
               if Is_Miller_Rabin_Witness(Witness,X) then return False;
               end if;
            end loop;
         end if;
      end;
      return True;
    end Pass_Prime_Test;

    ---------------------------------------------------------------------------


    function Is_Prime(X : Big_Unsigned) return Boolean is
    begin
       return Pass_Prime_Test(X, Strong);
    end Is_Prime; pragma Inline (Is_Prime);

    ---------------------------------------------------------------------------

    -- This function is faster then Is_prime but a lot of no strong pseudo
    -- primes pass this test
    function Looks_Like_A_Prime(X : Big_Unsigned) return Boolean is
    begin
      return Pass_Prime_Test(X, Weak);
    end Looks_Like_A_Prime; pragma Inline(Looks_Like_A_Prime);

   ---------------------------------------------------------------------------

    function Get_Prime(N : Big_Unsigned) return Big_Unsigned is
       Result : Big_Unsigned := Get_Random(N);
   begin
      if N <= Big_Unsigned_Two then
         raise Constraint_Error;
      end if;

      -- make sure that Result is odd
      Result.Number(0) := Result.Number(0) or 1;
      loop
         if Is_Prime(Result) then  return Result;
         else  Result := (Result+2) mod N ;
         end if;
      end loop;
   end Get_Prime;


   ---------------------------------------------------------------------------


   function "mod"(Left :  D_Big_Unsigned; Right : Big_Unsigned)
                 return Big_Unsigned;

   ---------------------------------------------------------------------------

   -- Result = Left * Right (mod N)
   function Mult(Left, Right, N : Big_Unsigned) return Big_Unsigned is
      T : DWord;
      Carry : Word := 0;
      R : D_Big_Unsigned;
   begin
      for I in 0..Left.Last_Index loop
         for J in 0..Right.Last_Index loop
            T := DWord(Left.Number(I)) *  DWord(Right.Number(J))
               + DWord(R.Number(I+J)) + DWord(Carry);

            R.Number(I+J) := Word(T and DWord(Word'Last));

            Carry:= Word(Shift_Right(T,Word'Size));
         end loop;
         R.Number(I+Right.Last_Index+1) := Carry +
          R.Number(I+Right.Last_Index+1);
         Carry := 0;
      end loop;

     for I in reverse 0..D_Max_Length loop
        if R.Number(I) /= 0 then
            R.Last_Index := I;
            exit;
         end if;
      end loop;
     return R mod N;
   end Mult;

 ---------------------------------------------------------------------------


   -- Returns a probability  N-bit prime (Result).
   function Get_N_Bit_Prime(N : Positive) return Big_Unsigned  is
      J : Big_Unsigned := Get_Random(Shift_Left(Big_Unsigned_One,N-2));
      Index : constant Natural := (N-1)/Word'Size;
      Amount : constant Natural := (N-1) mod Word'Size;
      Result : Big_Unsigned := Shift_Left(J,1);

   begin
      if N = 1 or N > Size then
         raise Constraint_Error;
      end if;

      loop
         -- Make sure that Result is an odd
         Set_Least_Significant_Bit (Result);

         -- Make sure that Result is a N-Bit-Number;
         Result.Number (Index) := Result.Number (Index) or
           Shift_Left (Word (1), Amount);

         if Amount = 0 then
	    Result.Last_Index := Index;
         end if;

         if Is_Prime(Result) then
	    return Result;
         else
            Result:=Result-2;
            if Is_Prime(Result) then 
	       return Result;
            end if;
         end if;

         J := Get_Random (Shift_Left (Big_Unsigned_One, N - 2));
         Result := Shift_Left (J, 1);
      end loop;

   end Get_N_Bit_Prime;

   ---------------------------------------------------------------------------

   -- computes the jacobi-symbol
   -- return value:
   --   0 : if X mod N = 0
   --   1 : if X is a quadratic resuide mod N
   --  -1 : if X is a quadratic non-resuide mod N

   function Jacobi(X, N : Big_Unsigned) return Integer is
      A : Big_Unsigned :=  X mod N;
   begin

      if Is_Even(N) then
         raise Constraint_Error;
      end if;

      if N = Big_Unsigned_One then return 1;
      elsif A = Big_Unsigned_Zero then return 0;
      elsif A =  Big_Unsigned_One then return 1;
      end if;

      while (A mod 4) = Big_Unsigned_Zero loop
         exit when (A mod 4) = Big_Unsigned_Zero;
         A := Shift_Right(A,2);
      end loop;

      if Is_Even(A) then
         if (N mod 8 = 1) or (N mod 8 = 7) then
            return Jacobi(Shift_Right(A,1),N);
         else return -1*Jacobi(Shift_Right(A,1),N);
         end if;
      else
         if (A mod 4 = 1)  or (N mod 4 = 1) then
            return Jacobi(N mod A, A);
         else  return -1*Jacobi(N mod A, A);
         end if;
      end if;
   end Jacobi;

  ----------------------------------------------------------------------------
  -----------------------------DOUBLE_SIZE------------------------------------
  ----------------------------------------------------------------------------

   --only needed for multiplication mod N
   --here we need 2*Size-bit numbers to avoid an overflow because
   --if one of our provisional result t > BIG_Unsigned_Last
   --then there ist no well known algortihm to compute the
   -- result of an multiplication mod m

   -- same algorithm for D_Big_Unsigned as for Big_Unsigned

   function "="(Left, Right : D_Big_Unsigned) return Boolean is
   begin
      if Left.Last_Index = Right.Last_Index then
         for I in 0..Left.Last_Index loop
            if  Left.Number(I) /= Right.Number(I) then return False;
            end if;
         end loop;
      else return False;
      end if;
      return True;
   end"=";

   ----------------------------------------------------------------------------

   function Shift_Left(Value : D_Big_Unsigned; Amount : Natural)
                      return D_Big_Unsigned is
   begin
      if Amount >= (D_Max_Length+1)*Word'Size or
        Value = D_Big_Unsigned_Zero
      then  return D_Big_Unsigned_Zero;
      elsif Amount = 0 then return Value;
      end if;

      declare
         Result : D_Big_Unsigned;
         Temp : DLimbs :=(others => 0);
         L : constant Natural := Amount mod Word'Size;
         R : constant Natural := Word'Size-L;
         M : constant Natural := Amount/Word'Size;
      begin
         Temp(0) := Shift_Left(Value.Number(0), L);

         for I in 1..Value.Last_Index loop
            Temp(I) := Shift_Right(Value.Number(I-1), R) +
              Shift_Left(Value.Number(I), L);
         end loop;

         if Value.Last_Index /= D_Max_Length then
            Temp(Value.Last_Index+1):=
              Shift_Right(Value.Number(Value.Last_Index), R);
         end if;

         for I in Temp'Range loop
            if (I+M) > D_Max_Length then
               exit;
            end if;
            Result.Number(I+M):= Temp(I);
         end loop;

         for I in reverse 0..D_Max_Length loop
            if Result.Number(I) /=0 then
               Result.Last_Index:=I;
               exit;
            end if;
         end loop;
         return Result;
      end;
   end Shift_Left; pragma Inline (Shift_Left);

   ---------------------------------------------------------------------------


   function Bit_Length(X : D_Big_Unsigned) return Natural is
   begin
      if X = D_Big_Unsigned_Zero then
         return 0;
      end if;

      for I in reverse 0..Word'Size-1 loop
         if Shift_Left(1,I) <= X.Number(X.Last_Index) then
            return  Word'Size * X.Last_Index + I + 1 ;
         end if;
      end loop;
      return X.Last_Index * Word'Size;
   end Bit_Length; pragma Inline (Bit_Length);


   ---------------------------------------------------------------------------

   function "<"(Left, Right : D_Big_Unsigned) return Boolean is
   begin
      if Left.Last_Index < Right.Last_Index then return True;
      elsif Left.Last_Index > Right.Last_Index then  return False;
      else
         for I in reverse  0..Left.Last_Index loop
            if  Left.Number(I) < Right.Number(I) then return True;
            elsif Left.Number(I) > Right.Number(I) then return False;
            end if;
         end loop;
      end if;
      return False;
   end "<"; pragma Inline ("<");

   ---------------------------------------------------------------------------

   function ">"(Left, Right : D_Big_Unsigned) return Boolean is
   begin
      return Right < Left;
   end ">"; pragma Inline (">");



   ---------------------------------------------------------------------------

   function ">="(Left, Right : D_Big_Unsigned) return Boolean is
   begin
      return not(Left < Right);
   end ">="; pragma Inline (">=");

   ---------------------------------------------------------------------------

   function "+"(Left, Right : D_Big_Unsigned) return D_Big_Unsigned;

   function "-"(Left, Right : D_Big_Unsigned) return D_Big_Unsigned is
   begin
      if Left = Right then  return D_Big_Unsigned_Zero;
      elsif Left = Right+D_Big_Unsigned_One then return D_Big_Unsigned_One;
      elsif Left+D_Big_Unsigned_One = Right then return D_Big_Unsigned_Last;

      -- add the modulus if Right > Left
      elsif Right > Left then
         return D_Big_Unsigned_Last - Right + Left + D_Big_Unsigned_One;
      else
         declare
            Result : D_Big_Unsigned;
            Carry : Word:=0;
         begin
            -- Remember Left > Right
            for I in 0..Left.Last_Index loop
               Result.Number(I) := Left.Number(I) - Right.Number(I) - Carry;
               if (Right.Number(I) >  Left.Number(I)) or
                 (Carry= 1 and Right.Number(I) =  Left.Number(I))
               then Carry :=1;
               else Carry :=0;
               end if;
               if Result.Number(I) /= 0 then
                  Result.Last_Index := I;
               end if;
            end loop;
            return Result;
         end;
      end if;
   end "-";


   ---------------------------------------------------------------------------

   function "mod"(Left :  D_Big_Unsigned; Right : Big_Unsigned)
                 return Big_Unsigned is
   begin
      if Left.Last_Index <= Max_Length then
         declare
            L : Big_Unsigned;
         begin
            L.Last_Index := Left.Last_Index;
            L.Number(0..Left.Last_Index) := Left.Number(0..Left.Last_Index);
         return L mod Right;
         end;
      end if;

      if Right = Big_Unsigned_Zero then raise Division_By_Zero;
      --elsif Right = Big_Unsigned_One then return Big_Unsigned_Zero;
      end if;

      -- Now, there is only the case where (Left > Right), (Right /= 0)
      -- and |Left|>|Right|.

      declare
         Remainder : D_Big_Unsigned:=Left;
         Temp_Right, R : D_Big_Unsigned;
         Result : Big_Unsigned;
         Diff: Natural;

      begin
         Temp_Right.Last_Index := Right.Last_Index;
         Temp_Right.Number(0..Right.Last_Index) :=
           Right.Number(0..Right.Last_Index);
         R:=Temp_Right;

         while(Remainder >= R) loop
            Diff := Bit_Length(Remainder) - Bit_Length(R);
            if Diff = 0 then
               Remainder := Remainder-R;
               exit;
            else Diff:=Diff-1;
            end if;
            Temp_Right := Shift_Left(R, Diff);
            Remainder :=  Remainder-Temp_Right;
         end loop;

         Result.Last_Index := Remainder.Last_Index;
         Result.Number(0..Result.Last_Index) :=
           Remainder.Number(0..Result.Last_Index);
         return Result;
      end;
   end "mod";

   ---------------------------------------------------------------------------

   function "+"(Left, Right : D_Big_Unsigned) return D_Big_Unsigned is
      Result : D_Big_Unsigned;
      M : constant Natural  := Natural'Max(Left.Last_Index, Right.Last_Index);
      Temp : Word;
      Carry : Word :=0;
   begin

        for I in 0..M loop
         Temp :=Carry;
         Result.Number(I) := Left.Number(I) + Right.Number(I) +Temp;
         if Result.Number(I) < Word'Max(Left.Number(I), Right.Number(I))
         then  Carry := 1;
         else Carry := 0;
         end if;
      end loop;

      if Carry  =1 and M < Max_Length then
         Result.Number(M+1) := 1;
         Result.Last_Index := M+1;
      else
         -- Set Result.Last_Index
         for I in reverse 0..M loop
            if Result.Number(I) /= 0 then
               Result.Last_Index := I;
               return Result;
            end if;
         end loop;
      end if;
      return Result;
   end "+";

   ---------------------------------------------------------------------------

end Mod_Utils;