Fast path for Z.divisible on small arguments (#147)

Closes: #140

tests/zq.ml

@@ -657,6 +657,16 @@ let test_Z() =
   Printf.printf "hamdist 2^120 (2^120-1)\n = %i\n" (I.hamdist p120 (I.pred p120));
   Printf.printf "hamdist 2^120 2^300\n = %i\n" (I.hamdist p120 p300);
   Printf.printf "hamdist (2^120-1) (2^300-1)\n = %i\n" (I.hamdist (I.pred p120) (I.pred p300));
+  Printf.printf "divisible 42 7\n = %B\n" (I.divisible (I.of_int 42) (I.of_int 7));
+  Printf.printf "divisible 43 7\n = %B\n" (I.divisible (I.of_int 43) (I.of_int 7));
+  Printf.printf "divisible 0 0\n = %B\n" (I.divisible I.zero I.zero);
+  Printf.printf "divisible 0 2^120\n = %B\n" (I.divisible I.zero p120);
+  Printf.printf "divisible 2 2^120\n = %B\n" (I.divisible (I.of_int 2) p120);
+  Printf.printf "divisible 2^300 2^120\n = %B\n" (I.divisible p300 p120);
+  Printf.printf "divisible (2^300-1) 32\n = %B\n" (I.divisible (I.pred p300) (I.of_int 32));
+  Printf.printf "divisible min_int (max_int+1)\n = %B\n" (I.divisible (I.of_int min_int) (I.succ (I.of_int max_int)));
+  Printf.printf "divisible (max_int+1) min_int\n = %B\n" (I.divisible (I.succ (I.of_int max_int)) (I.of_int min_int));
+
   (* always 0 when not using custom blocks *)
   Printf.printf "hash(2^120)\n = %i\n" (Hashtbl.hash p120);
   Printf.printf "hash(2^121)\n = %i\n" (Hashtbl.hash p121);

tests/zq.output32

@@ -903,6 +903,24 @@ hamdist 2^120 2^300
  = 2
 hamdist (2^120-1) (2^300-1)
  = 180
+divisible 42 7
+ = true
+divisible 43 7
+ = false
+divisible 0 0
+ = true
+divisible 0 2^120
+ = true
+divisible 2 2^120
+ = false
+divisible 2^300 2^120
+ = true
+divisible (2^300-1) 32
+ = false
+divisible min_int (max_int+1)
+ = true
+divisible (max_int+1) min_int
+ = true
 hash(2^120)
  = 691199303
 hash(2^121)

tests/zq.output64

@@ -903,6 +903,24 @@ hamdist 2^120 2^300
  = 2
 hamdist (2^120-1) (2^300-1)
  = 180
+divisible 42 7
+ = true
+divisible 43 7
+ = false
+divisible 0 0
+ = true
+divisible 0 2^120
+ = true
+divisible 2 2^120
+ = false
+divisible 2^300 2^120
+ = true
+divisible (2^300-1) 32
+ = false
+divisible min_int (max_int+1)
+ = true
+divisible (max_int+1) min_int
+ = true
 hash(2^120)
  = 691199303
 hash(2^121)

z.ml

@@ -261,7 +261,27 @@ external nextprime: t -> t = "ml_z_nextprime"
 external hash: t -> int = "ml_z_hash" [@@noalloc]
 external to_bits: t -> string = "ml_z_to_bits"
 external of_bits: string -> t = "ml_z_of_bits"
-external divisible: t -> t -> bool = "ml_z_divisible"
+
+external c_divisible: t -> t -> bool = "ml_z_divisible"
+
+let divisible x y =
+  if is_small_int x then
+    if is_small_int y then
+      if unsafe_to_int y = 0
+      then unsafe_to_int x = 0
+      else (unsafe_to_int x) mod (unsafe_to_int y) = 0
+    else
+      (* If y divides x, we have |y| <= |x| or x = 0.
+         Here, x is small:  min_int <= x <= max_int
+         and y is not small:  y < min_int \/ y > max_int.
+         |y| <= |x| is possible only if
+         x = min_int and y = -min_int = max_int+1 .
+         So, the only two cases where y divides x are
+         x = 0 or x = min_int /\ y = -min_int. *)
+      unsafe_to_int x = 0 || (unsafe_to_int x = min_int && y = c_neg x)
+  else
+    c_divisible x y
+
 external congruent: t -> t -> t -> bool = "ml_z_congruent"
 external jacobi: t -> t -> int = "ml_z_jacobi"
 external legendre: t -> t -> int = "ml_z_legendre"

z.mli

@@ -209,7 +209,7 @@ val divexact: t -> t -> t
     Can raise a [Division_by_zero].
 *)
 
-external divisible: t -> t -> bool = "ml_z_divisible"
+val divisible: t -> t -> bool
 (** [divisible a b] returns [true] if [a] is exactly divisible by [b].
     Unlike the other division functions, [b = 0] is accepted 
     (only 0 is considered divisible by 0).