Skip to content

Commit

Permalink
[micromega] Workaround bug / discordance in gcd calculation.
Browse files Browse the repository at this point in the history
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
  • Loading branch information
ejgallego committed Nov 17, 2018
1 parent c747a7d commit 41e06dd
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 6 deletions.
11 changes: 9 additions & 2 deletions plugins/micromega/sos.ml
Original file line number Diff line number Diff line change
Expand Up @@ -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'

(* ------------------------------------------------------------------------- *)
Expand Down
15 changes: 11 additions & 4 deletions plugins/micromega/vect.ml
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down

0 comments on commit 41e06dd

Please sign in to comment.