Merge pull request #39 from vbgl/fix-gcd-zero

Fix Z.gcd when exactly one argument is zero
ocaml · Mar 4, 2019 · 9ab3443 · 9ab3443
2 parents 6f84992 + e54936d
commit 9ab3443
Show file tree

Hide file tree

Showing 4 changed files with 58 additions and 46 deletions.
diff --git a/caml_z.c b/caml_z.c
@@ -1740,7 +1740,6 @@ CAMLprim value ml_z_gcd(value arg1, value arg2)
     /* fast path */
     intnat a1 = Long_val(arg1);
     intnat a2 = Long_val(arg2);
-    if (!a1 || !a2) ml_z_raise_divide_by_zero();
     if (a1 < 0) a1 = -a1;
     if (a2 < 0) a2 = -a2;
     if (a1 < a2) { intnat t = a1; a1 = a2; a2 = t; }
@@ -1759,52 +1758,55 @@ CAMLprim value ml_z_gcd(value arg1, value arg2)
     mp_size_t sz, pos1, pos2, limb1, limb2, bit1, bit2, pos, limb, bit, i;
     Z_DECL(arg1); Z_DECL(arg2);
     Z_ARG(arg1);  Z_ARG(arg2);
-    if (!size_arg1 || !size_arg2) ml_z_raise_divide_by_zero();
-    /* copy args to tmp storage & remove lower 0 bits */
-    pos1 = mpn_scan1(ptr_arg1, 0);
-    pos2 = mpn_scan1(ptr_arg2, 0);
-    limb1 = pos1 / Z_LIMB_BITS;
-    limb2 = pos2 / Z_LIMB_BITS;
-    bit1 = pos1 % Z_LIMB_BITS;
-    bit2 = pos2 % Z_LIMB_BITS;
-    size_arg1 -= limb1;
-    size_arg2 -= limb2;
-    tmp1 = ml_z_alloc(size_arg1 + 1);
-    tmp2 = ml_z_alloc(size_arg2 + 1);
-    Z_REFRESH(arg1);
-    Z_REFRESH(arg2);
-    if (bit1) {
-      mpn_rshift(Z_LIMB(tmp1), ptr_arg1 + limb1, size_arg1, bit1);
-      if (!Z_LIMB(tmp1)[size_arg1-1]) size_arg1--;
-    }
-    else ml_z_cpy_limb(Z_LIMB(tmp1), ptr_arg1 + limb1, size_arg1);
-    if (bit2) { 
-      mpn_rshift(Z_LIMB(tmp2), ptr_arg2 + limb2, size_arg2, bit2);
-      if (!Z_LIMB(tmp2)[size_arg2-1]) size_arg2--;
-    }
-    else ml_z_cpy_limb(Z_LIMB(tmp2), ptr_arg2 + limb2, size_arg2);
-    /* compute gcd of 2^pos1 & 2^pos2 */
-    pos = (pos1 <= pos2) ? pos1 : pos2;
-    limb = pos / Z_LIMB_BITS;
-    bit = pos % Z_LIMB_BITS;
-    /* compute gcd of arg1 & arg2 without lower 0 bits */
-    /* second argument must have less bits than first  */
-    if ((size_arg1 > size_arg2) ||
-        ((size_arg1 == size_arg2) && 
-         (Z_LIMB(tmp1)[size_arg1 - 1] >= Z_LIMB(tmp2)[size_arg1 - 1]))) {
-      r = ml_z_alloc(size_arg2 + limb + 1);
-      sz = mpn_gcd(Z_LIMB(r) + limb, Z_LIMB(tmp1), size_arg1, Z_LIMB(tmp2), size_arg2);
-    }
+    if (!size_arg1) r = arg2;
+    else if (!size_arg2) r = arg1;
     else {
-      r = ml_z_alloc(size_arg1 + limb + 1);
-      sz = mpn_gcd(Z_LIMB(r) + limb, Z_LIMB(tmp2), size_arg2, Z_LIMB(tmp1), size_arg1);
-    } 
-    /* glue the two results */
-    for (i = 0; i < limb; i++) 
-      Z_LIMB(r)[i] = 0;
-    Z_LIMB(r)[sz + limb] = 0;
-    if (bit) mpn_lshift(Z_LIMB(r) + limb, Z_LIMB(r) + limb, sz + 1, bit);
-    r = ml_z_reduce(r, limb + sz + 1, 0);
+      /* copy args to tmp storage & remove lower 0 bits */
+      pos1 = mpn_scan1(ptr_arg1, 0);
+      pos2 = mpn_scan1(ptr_arg2, 0);
+      limb1 = pos1 / Z_LIMB_BITS;
+      limb2 = pos2 / Z_LIMB_BITS;
+      bit1 = pos1 % Z_LIMB_BITS;
+      bit2 = pos2 % Z_LIMB_BITS;
+      size_arg1 -= limb1;
+      size_arg2 -= limb2;
+      tmp1 = ml_z_alloc(size_arg1 + 1);
+      tmp2 = ml_z_alloc(size_arg2 + 1);
+      Z_REFRESH(arg1);
+      Z_REFRESH(arg2);
+      if (bit1) {
+        mpn_rshift(Z_LIMB(tmp1), ptr_arg1 + limb1, size_arg1, bit1);
+        if (!Z_LIMB(tmp1)[size_arg1-1]) size_arg1--;
+      }
+      else ml_z_cpy_limb(Z_LIMB(tmp1), ptr_arg1 + limb1, size_arg1);
+      if (bit2) { 
+        mpn_rshift(Z_LIMB(tmp2), ptr_arg2 + limb2, size_arg2, bit2);
+        if (!Z_LIMB(tmp2)[size_arg2-1]) size_arg2--;
+      }
+      else ml_z_cpy_limb(Z_LIMB(tmp2), ptr_arg2 + limb2, size_arg2);
+      /* compute gcd of 2^pos1 & 2^pos2 */
+      pos = (pos1 <= pos2) ? pos1 : pos2;
+      limb = pos / Z_LIMB_BITS;
+      bit = pos % Z_LIMB_BITS;
+      /* compute gcd of arg1 & arg2 without lower 0 bits */
+      /* second argument must have less bits than first  */
+      if ((size_arg1 > size_arg2) ||
+          ((size_arg1 == size_arg2) && 
+           (Z_LIMB(tmp1)[size_arg1 - 1] >= Z_LIMB(tmp2)[size_arg1 - 1]))) {
+        r = ml_z_alloc(size_arg2 + limb + 1);
+        sz = mpn_gcd(Z_LIMB(r) + limb, Z_LIMB(tmp1), size_arg1, Z_LIMB(tmp2), size_arg2);
+      }
+      else {
+        r = ml_z_alloc(size_arg1 + limb + 1);
+        sz = mpn_gcd(Z_LIMB(r) + limb, Z_LIMB(tmp2), size_arg2, Z_LIMB(tmp1), size_arg1);
+      } 
+      /* glue the two results */
+      for (i = 0; i < limb; i++) 
+        Z_LIMB(r)[i] = 0;
+      Z_LIMB(r)[sz + limb] = 0;
+      if (bit) mpn_lshift(Z_LIMB(r) + limb, Z_LIMB(r) + limb, sz + 1, bit);
+      r = ml_z_reduce(r, limb + sz + 1, 0);
+    }
     Z_CHECK(r);
     CAMLreturn(r);
   }

diff --git a/tests/zq.ml b/tests/zq.ml
@@ -505,12 +505,14 @@ let test_Z() =
   Printf.printf "sign -1\n = %i\n" (I.sign I.minus_one);
   Printf.printf "sign 2^300\n = %i\n" (I.sign p300);
   Printf.printf "sign -2^300\n = %i\n" (I.sign (I.neg p300));
+  Printf.printf "gcd 0 -137\n = %a\n" pr (I.gcd (I.of_int 0) (I.of_int (-137)));
   Printf.printf "gcd 12 27\n = %a\n" pr (I.gcd (I.of_int 12) (I.of_int 27));
   Printf.printf "gcd 27 12\n = %a\n" pr (I.gcd (I.of_int 27) (I.of_int 12));
   Printf.printf "gcd 27 27\n = %a\n" pr (I.gcd (I.of_int 27) (I.of_int 27));
   Printf.printf "gcd -12 27\n = %a\n" pr (I.gcd (I.of_int (-12)) (I.of_int 27));
   Printf.printf "gcd 12 -27\n = %a\n" pr (I.gcd (I.of_int 12) (I.of_int (-27)));
   Printf.printf "gcd -12 -27\n = %a\n" pr (I.gcd (I.of_int (-12)) (I.of_int (-27)));
+  Printf.printf "gcd 0 2^300\n = %a\n" pr (I.gcd (I.of_int 0) p300);
   Printf.printf "gcd 2^120 2^300\n = %a\n" pr (I.gcd p120 p300);
   Printf.printf "gcd 2^300 2^120\n = %a\n" pr (I.gcd p300 p120);
   Printf.printf "gcdext 12 27\n = %a\n" pr3 (I.gcdext (I.of_int 12) (I.of_int 27));

diff --git a/tests/zq.output32 b/tests/zq.output32
@@ -700,6 +700,8 @@ sign 2^300
  = 1
 sign -2^300
  = -1
+gcd 0 -137
+ = 137
 gcd 12 27
  = 3
 gcd 27 12
@@ -712,6 +714,8 @@ gcd 12 -27
  = 3
 gcd -12 -27
  = 3
+gcd 0 2^300
+ = 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397376
 gcd 2^120 2^300
  = 1329227995784915872903807060280344576
 gcd 2^300 2^120

diff --git a/tests/zq.output64 b/tests/zq.output64
@@ -700,6 +700,8 @@ sign 2^300
  = 1
 sign -2^300
  = -1
+gcd 0 -137
+ = 137
 gcd 12 27
  = 3
 gcd 27 12
@@ -712,6 +714,8 @@ gcd 12 -27
  = 3
 gcd -12 -27
  = 3
+gcd 0 2^300
+ = 2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397376
 gcd 2^120 2^300
  = 1329227995784915872903807060280344576
 gcd 2^300 2^120