-- 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;