From 41e06ddd55a5b4ea2c00eeadb3d3074889c384e8 Mon Sep 17 00:00:00 2001
From: Emilio Jesus Gallego Arias <e+git@x80.org>
Date: Wed, 17 Oct 2018 18:40:03 +0200
Subject: [PATCH] [micromega] Workaround bug / discordance in gcd calculation.

We workaround a bug in Zarith's `gcd`. Indeed, `gcd n 0` when `n != 0`
should be `n`. The problem can be witnessed here:
```
utop # Big_int_Z.(gcd_big_int zero_big_int unit_big_int);;
Exception: Division_by_zero.

utop # Big_int.(gcd_big_int zero_big_int unit_big_int);;
- : Big_int.big_int = <big_int 1>
```

c.f: https://github.com/ocaml/Zarith/pull/39
---
 plugins/micromega/sos.ml  | 11 +++++++++--
 plugins/micromega/vect.ml | 15 +++++++++++----
 2 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/plugins/micromega/sos.ml b/plugins/micromega/sos.ml
index 162171d3e4af8..5ca349d1506ef 100644
--- a/plugins/micromega/sos.ml
+++ b/plugins/micromega/sos.ml
@@ -634,19 +634,26 @@ let diag m =
 (* Adjust a diagonalization to collect rationals at the start.               *)
 (* ------------------------------------------------------------------------- *)
 
+(* Workaround to removed once https://github.com/ocaml/Zarith/pull/39
+   is fixed *)
+let fix_gcd c n =
+  if Z.(equal c zero) then
+    if Z.(equal n zero) then Z.(div zero zero) else n
+  else Z.gcd c n
+
 let deration d =
   if d = [] then Q.zero, d else
   let adj(c,l) =
     let a =
       Q.make
         (foldl (fun a i c -> Z.lcm a (Q.den c)) Z.one (snd l))
-        (foldl (fun a i c -> Z.gcd a (Q.num c)) Z.zero (snd l)) in
+        (foldl (fun a i c -> fix_gcd a (Q.num c)) Z.zero (snd l)) in
     (Q.div c (Q.mul a a)), mapa (Q.mul a) l in
   let d' = List.map adj d in
   let a =
     Q.make
       (List.fold_right ((o) Z.lcm ((o) Q.den fst)) d' Z.one)
-      (List.fold_right ((o) Z.gcd ((o) Q.num fst)) d' Z.zero) in
+      (List.fold_right ((o) fix_gcd ((o) Q.num fst)) d' Z.zero) in
   Q.inv a, List.map (fun (c,l) -> (Q.mul a c, l)) d'
 
 (* ------------------------------------------------------------------------- *)
diff --git a/plugins/micromega/vect.ml b/plugins/micromega/vect.ml
index 43c096ea96fea..9100fb5e02372 100644
--- a/plugins/micromega/vect.ml
+++ b/plugins/micromega/vect.ml
@@ -250,17 +250,24 @@ let incr_var i v = List.map (fun (v,n) -> (v+i,n)) v
 
 open Big_int_Z
 
+(* Workaround to removed once https://github.com/ocaml/Zarith/pull/39
+   is fixed *)
+let fix_gcd c n =
+  if Z.(equal c zero) then
+    if Z.(equal n zero) then Z.(div zero zero) else n
+  else gcd_big_int c n
+
 let gcd v =
   let res = fold (fun c _ n  ->
-                assert (Int.equal (compare_big_int (Q.den n)  unit_big_int) 0);
-                gcd_big_int c (Q.num n)) zero_big_int v in
+                assert (Int.equal (compare_big_int (Q.den n) unit_big_int) 0);
+                fix_gcd c (Q.num n)) zero_big_int v in
   if Int.equal (compare_big_int res zero_big_int) 0
   then unit_big_int else res
 
 let normalise v =
   let ppcm = fold (fun c _ n -> ppcm c (Q.den n)) unit_big_int v in
-  let gcd  =
-    let gcd = fold (fun c _ n -> gcd_big_int c (Q.num n)) zero_big_int v in
+  let gcd =
+    let gcd = fold (fun c _ n -> fix_gcd c (Q.num n)) unit_big_int v in
     if Int.equal (compare_big_int gcd zero_big_int) 0 then unit_big_int else gcd in
   List.map (fun (x,v) -> (x, Q.(mul v (div (of_bigint ppcm) (of_bigint gcd))))) v