Ajout de l'exponentiation modulaire

git-svn-id: http://caml.inria.fr/svn/ocaml/trunk@4492 f963ae5c-01c2-4b8c-9fe0-0dff7051ff02

5 files changed, 72 insertions, 4 deletions

diff --git a/otherlibs/num/big_int.ml b/otherlibs/num/big_int.ml
index 7542f1f89..7b1091096 100644
--- a/otherlibs/num/big_int.ml
+++ b/otherlibs/num/big_int.ml
@@ -531,6 +531,19 @@ let base_power_big_int base n bi =
                then zero_big_int
                else create_big_int (bi.sign) res
 
+(* Modular exponentiation *)
+
+let modexp_big_int a b c =
+  if b.sign < 0 || c.sign < 0 then invalid_arg "modexp_big_int";
+  if c.sign = 0 then raise Division_by_zero;
+  let na =
+    if a.sign < 0 || ge_big_int a c
+    then mod_big_int a c
+    else a in
+  let res_nat = modexp_nat na.abs_value b.abs_value c.abs_value in
+  { sign = if is_zero_nat res_nat 0 (length_nat res_nat) then 0 else 1;
+    abs_value = res_nat }
+
 (* Coercion with float type *)
 
 let float_of_big_int bi = 
diff --git a/otherlibs/num/big_int.mli b/otherlibs/num/big_int.mli
index 9b140abf2..a92354c47 100644
--- a/otherlibs/num/big_int.mli
+++ b/otherlibs/num/big_int.mli
@@ -74,10 +74,15 @@ val power_big_int_positive_int: big_int -> int -> big_int
 val power_int_positive_big_int: int -> big_int -> big_int
 val power_big_int_positive_big_int: big_int -> big_int -> big_int
         (** Exponentiation functions.  Return the big integer
-           representing the first argument [a] raised to the power [b]
-           (the second argument).  Depending
-           on the function, [a] and [b] can be either small integers
-           or big integers.  Raise [Invalid_argument] if [b] is negative. *)
+            representing the first argument [a] raised to the power [b]
+            (the second argument).  Depending
+            on the function, [a] and [b] can be either small integers
+            or big integers.  Raise [Invalid_argument] if [b] is negative. *)
+val modexp_big_int: big_int -> big_int -> big_int -> big_int
+        (** Modular exponentiation.  [modexp_big_int a b c] returns
+            [a] to the [b]-th power, modulo [c].  Raise [Invalid_argument]
+            if [b] or [c] are negative, and [Division_by_zero] if [c]
+            is zero. *)
 
 (** {6 Comparisons and tests} *)
 
diff --git a/otherlibs/num/nat.ml b/otherlibs/num/nat.ml
index 02f1b1b7a..3cc0fca4b 100644
--- a/otherlibs/num/nat.ml
+++ b/otherlibs/num/nat.ml
@@ -234,6 +234,34 @@ let sqrt_nat rad off len =
   loop ()
  end;;
 
+(* Modular exponentiation.  Return a fresh nat equal to (a ^ b) mod c.
+   We assume c > 0 and a < c. *)
+
+let modexp_nat a b c =
+  let len_c = length_nat c in
+  let res = make_nat len_c in
+  set_digit_nat res 0 1;
+  let prod = create_nat (2 * len_c) in
+  let modmult_res x =
+    (* res <- res * x mod c *)
+    let len_x = length_nat x in
+    set_to_zero_nat prod 0 (2 * len_c);
+    ignore(mult_nat prod 0 (len_c + len_x) res 0 len_c x 0 len_x);
+    div_nat prod 0 (len_c + len_x) c 0 len_c;
+    blit_nat res 0 prod 0 len_c in
+  let digit = make_nat 1
+  and carry = make_nat 1 in
+  (* Iterate over each bit of b, from most significant to least significant *)
+  for i = length_nat b - 1 downto 0 do
+    blit_nat digit 0 b i 1;
+    for i = 1 to length_of_digit do
+      modmult_res res;                  (* res <- res * res *)
+      shift_left_nat digit 0 1 carry 0 1;
+      if is_digit_odd carry 0 then modmult_res a (* res <- res * a *)
+    done
+  done;
+  res
+
 let power_base_max = make_nat 2;;
 
 match length_of_digit with
diff --git a/otherlibs/num/nat.mli b/otherlibs/num/nat.mli
index df52473c0..a024e927a 100644
--- a/otherlibs/num/nat.mli
+++ b/otherlibs/num/nat.mli
@@ -63,6 +63,7 @@ external lxor_digit_nat: nat -> int -> nat -> int -> unit = "lxor_digit_nat"
 val square_nat : nat -> int -> int -> nat -> int -> int -> int
 val gcd_nat : nat -> int -> int -> nat -> int -> int -> int
 val sqrt_nat : nat -> int -> int -> nat
+val modexp_nat : nat -> nat -> nat -> nat
 val string_of_nat : nat -> string
 val nat_of_string : string -> nat
 val sys_nat_of_string : int -> string -> int -> int -> nat
diff --git a/otherlibs/num/test/test_big_ints.ml b/otherlibs/num/test/test_big_ints.ml
index 61e9ae4df..2335d2edd 100644
--- a/otherlibs/num/test/test_big_ints.ml
+++ b/otherlibs/num/test/test_big_ints.ml
@@ -466,3 +466,24 @@ test 3 eq_big_int
  (square_big_int (big_int_of_string "-1"), big_int_of_string "1");;
 test 4 eq_big_int
  (square_big_int (big_int_of_string "-7"), big_int_of_string "49");;
+
+testing_function "modexp_big_int";;
+
+test 1 eq_big_int
+  (modexp_big_int (big_int_of_int 0) (big_int_of_int 0) (big_int_of_int 42),
+   big_int_of_int 1);;
+test 2 eq_big_int
+  (modexp_big_int (big_int_of_int 0) (big_int_of_int 1) (big_int_of_int 42),
+   big_int_of_int 0);;
+test 3 eq_big_int
+  (modexp_big_int (big_int_of_int 12345) (big_int_of_int 6789)
+                  (big_int_of_string "4782546123567"),
+   big_int_of_string "3726846414024");;
+test 4 eq_big_int
+  (modexp_big_int (big_int_of_int (-12345)) (big_int_of_int 6789)
+                  (big_int_of_string "4782546123567"),
+   big_int_of_string "1055699709543");;
+test 5 eq_big_int
+  (modexp_big_int (big_int_of_string "200866937202159") (big_int_of_int 6789)
+                  (big_int_of_string "4782546123567"),
+   big_int_of_string "3726846414024");;