[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: ocaml/Zarith#39

plugins/micromega/certificate.ml

@@ -527,7 +527,7 @@ let rec scale_term t =
  (s, Opp t)
  | Add (t1, t2) ->
  let s1, y1 = scale_term t1 and s2, y2 = scale_term t2 in
- let g = gcd_big_int s1 s2 in
+ let g = fix_gcd s1 s2 in
  let s1' = div_big_int s1 g in
  let s2' = div_big_int s2 g in
  let e = mult_big_int g (mult_big_int s1' s2') in
@@ -567,7 +567,7 @@ let rec scale_certificate pos =
  (mult_big_int s1 s2, Eqmul (y1, y2))
  | Sum (y, z) ->
  let s1, y1 = scale_certificate y and s2, y2 = scale_certificate z in
- let g = gcd_big_int s1 s2 in
+ let g = fix_gcd s1 s2 in
  let s1' = div_big_int s1 g in
  let s2' = div_big_int s2 g in
  ( mult_big_int g (mult_big_int s1' s2')
@@ -717,7 +717,7 @@ let extract_coprime (c1, p1) (c2, p2) =
  Vect.exists2
  (fun n1 n2 ->
  Int.equal
- (compare_big_int (gcd_big_int (Q.num n1) (Q.den n2)) unit_big_int)
+ (compare_big_int (fix_gcd (Q.num n1) (Q.den n2)) unit_big_int)
  0)
  c1.coeffs c2.coeffs
  else None
@@ -779,7 +779,7 @@ let reduce_var_change psys =
  Vect.find
  (fun x' v' ->
  let v' = Q.num v' in
- if eq_big_int (gcd_big_int v v') unit_big_int then Some (x', v')
+ if eq_big_int (fix_gcd v v') unit_big_int then Some (x', v')
  else None)
  vect
  with

plugins/micromega/mutils.ml

@@ -161,8 +161,12 @@ let saturate_bin (f : 'a -> 'a -> 'a option) (l : 'a list) =
 
 open Big_int_Z
 
+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 ppcm x y =
- let g = gcd_big_int x y in
+ let g = fix_gcd x y in
  let x' = div_big_int x g in
  let y' = div_big_int y g in
  mult_big_int g (mult_big_int x' y')
@@ -456,6 +460,7 @@ module NumNotation = struct
  let z10 = Z.of_int 10
  let pow2 x = Z.pow z2 x |> Q.of_bigint
  let pow10 x = Z.pow z10 x |> Q.of_bigint
+ let fix_gcd c n = fix_gcd c n
 end
 
 (* Local Variables: *)

plugins/micromega/mutils.mli

@@ -135,4 +135,8 @@ module NumNotation : sig
  val round : Q.t -> Q.t
  val pow2 : int -> Q.t
  val pow10 : int -> Q.t
+
+ (* Workaround to removed once https://github.com/ocaml/Zarith/pull/39
+ is fixed *)
+ val fix_gcd : Z.t -> Z.t -> Z.t
 end

plugins/micromega/sos.ml

@@ -669,15 +669,15 @@ let deration d =
  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))
+ (foldl (fun a i c -> fix_gcd a (Q.num c)) Z.zero (snd l))
  in
  (c // (a */ a), mapa (fun x -> a */ x) 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)
+ (List.fold_right (o fix_gcd (o Q.num fst)) d' Z.zero)
  in
  (Q.one // a, List.map (fun (c, l) -> (a */ c, l)) d')
 

plugins/micromega/vect.ml

@@ -238,15 +238,15 @@ let gcd v =
  fold
  (fun c _ n ->
  assert (Int.equal (compare_big_int (Q.den n) unit_big_int) 0);
- gcd_big_int c (Q.num n))
+ NumNotation.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 = fold (fun c _ n -> fix_gcd c (Q.num n)) zero_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, v */ Q.of_bigint ppcm // Q.of_bigint gcd)) v