diff --git a/libpoly/ocaml_poly.ml b/libpoly/ocaml_poly.ml
index 5539f21f05536260426a9d2b13769266111a95c0..51c49af0b016040c1625071aa200d2503115c864 100644
--- a/libpoly/ocaml_poly.ml
+++ b/libpoly/ocaml_poly.ml
@@ -133,6 +133,9 @@ module A = struct
     binop C.Function.lp_algebraic_number_positive_root a
       (Unsigned.UInt.of_int n)
 
+  let positive_pow a (q : Q.t) =
+    positive_root (pow a (Z.to_int q.num)) (Z.to_int q.den)
+
   let is_integer a = C.Function.lp_algebraic_number_is_integer a
   let is_rational a = C.Function.lp_algebraic_number_is_rational a
 
diff --git a/libpoly/ocaml_poly.mli b/libpoly/ocaml_poly.mli
index afe21a80f1b287b100cdeeaed8eef62968c80518..89d5984195bdc5f8ab2a14b317b080f2880037ed 100644
--- a/libpoly/ocaml_poly.mli
+++ b/libpoly/ocaml_poly.mli
@@ -46,6 +46,9 @@ module A : sig
   val mod_f : t -> t -> t
   val pow : t -> int -> t
   val positive_root : t -> int -> t
+
+  val positive_pow : t -> Q.t -> t
+  (** the numerator and denominator should fit in an int *)
 end
 
 (** Polynomial in one variable *)
diff --git a/src_colibri2/stdlib/std.ml b/src_colibri2/stdlib/std.ml
index c666f604005085ca4485229c16a69abb79517766..5d9e610f10b55c5f3afc7395db4d3f09aa99f88b 100644
--- a/src_colibri2/stdlib/std.ml
+++ b/src_colibri2/stdlib/std.ml
@@ -277,4 +277,6 @@ module A = struct
         empty
 
   let pow q n = match q with Q q -> Q (Q.pow q n) | A a -> !!(A.pow a n)
+  let positive_root q n = !!(A.positive_root (to_a q) n)
+  let positive_pow q n = !!(A.positive_pow (to_a q) n)
 end
diff --git a/src_colibri2/stdlib/std.mli b/src_colibri2/stdlib/std.mli
index 0a461204a0c75779b24988bbdc2cdd4dc80145ae..07e0c10974bd02a9a078d768fa7804af9541372f 100644
--- a/src_colibri2/stdlib/std.mli
+++ b/src_colibri2/stdlib/std.mli
@@ -154,4 +154,8 @@ module A : sig
   val ( / ) : t -> t -> t
   val min : t -> t -> t
   val max : t -> t -> t
+  val positive_root : t -> int -> t
+
+  val positive_pow : t -> Q.t -> t
+  (** the numerator and denominator should fit in an int *)
 end
diff --git a/src_colibri2/theories/LRA/dom_product.ml b/src_colibri2/theories/LRA/dom_product.ml
index 1b61a57f466445fcf24be0c77459427e8527d60d..bb2bfced681703ebd413ff953c62dabfcd8eb40d 100644
--- a/src_colibri2/theories/LRA/dom_product.ml
+++ b/src_colibri2/theories/LRA/dom_product.ml
@@ -297,7 +297,7 @@ let converter d (f : Ground.t) =
       reg a;
       reg b;
       SolveAbs.assume_equality d res
-        (Product.of_list [ (a, Q.one); (b, Q.one) ]);
+        (Product.of_list A.one [ (a, Q.one); (b, Q.one) ]);
       SolveSign.assume_equality d res
         (Sign_product.mul
            (Sign_product.of_one_node a)
@@ -310,7 +310,7 @@ let converter d (f : Ground.t) =
           if Dom_interval.is_not_zero d den then (
             (* num = res * den *)
             SolveAbs.assume_equality d num
-              (Product.of_list [ (res, Q.one); (den, Q.one) ]);
+              (Product.of_list A.one [ (res, Q.one); (den, Q.one) ]);
             SolveSign.assume_equality d num
               (Sign_product.mul
                  (Sign_product.of_one_node res)
@@ -328,7 +328,7 @@ let converter d (f : Ground.t) =
   | { app = { builtin = Expr.Abs }; tyargs = []; args; _ } ->
       let a = IArray.extract1_exn args in
       reg a;
-      SolveAbs.assume_equality d res (Product.of_list [ (a, Q.one) ]);
+      SolveAbs.assume_equality d res (Product.of_list A.one [ (a, Q.one) ]);
       SolveSign.assume_equality d res Sign_product.one
   | _ -> ()
 
@@ -345,7 +345,7 @@ let init env =
           if A.equal p.cst A.zero && Node.M.is_num_elt 1 p.poly then
             let cl, q = Node.M.choose p.poly in
             if A.equal q A.one then (
-              let p1 = Product.of_list [ (cl, Q.one) ] in
+              let p1 = Product.of_list A.one [ (cl, Q.one) ] in
               Debug.dprintf4 debug "[Polynome->Product] propagate %a is %a"
                 Node.pp n Node.pp cl;
               SolveAbs.assume_equality d n p1;
diff --git a/src_colibri2/theories/LRA/fourier.ml b/src_colibri2/theories/LRA/fourier.ml
index 57a14ad7851ba03f1643aba69097a58cceac6270..afc23ecee50033504762700ddaa4eb8102d88194 100644
--- a/src_colibri2/theories/LRA/fourier.ml
+++ b/src_colibri2/theories/LRA/fourier.ml
@@ -75,7 +75,7 @@ let divide d (p : Polynome.t) =
             Node.M.empty common
         in
         if Node.M.is_empty common then raise Exit;
-        let common = Product.of_map common in
+        let common = Product.of_map A.one common in
         Debug.dprintf4 debug "[Fourier] found possible division: %a / |%a|"
           Polynome.pp p Product.pp common;
         let cst, l =
diff --git a/src_colibri2/theories/LRA/product.ml b/src_colibri2/theories/LRA/product.ml
index 2652156e8e8dae84b6fe49b0c37793c221461bfa..2139dd1a235600ee11297ccbee197ade713cc68a 100644
--- a/src_colibri2/theories/LRA/product.ml
+++ b/src_colibri2/theories/LRA/product.ml
@@ -22,7 +22,7 @@ open Colibri2_popop_lib
 open Colibri2_core
 
 module T = struct
-  type t = { (* cst : Q.t; *) poly : Q.t Node.M.t } [@@deriving eq, ord, hash]
+  type t = { cst : A.t; poly : Q.t Node.M.t } [@@deriving eq, ord, hash]
 
   let pp_exposant fmt q =
     if Z.equal q.Q.den Z.one && Z.fits_int q.Q.num then
@@ -57,14 +57,13 @@ include Popop_stdlib.MkDatatype (T)
 (** different invariant *)
 
 let invariant p =
-  (* not (Q.equal p.cst Q.zero && Node.M.is_empty p.poly) && *)
-  Node.M.for_all (fun _ q -> not (Q.equal q Q.zero)) p.poly
+  A.sign p.cst > 0 && Node.M.for_all (fun _ q -> not (Q.equal q Q.zero)) p.poly
 
 (* let cst q = assert (not (Q.equal q Q.zero)); {cst = q; poly = Node.M.empty}
  * (\* let zero = cst Q.zero *\) *)
 
 (** constructor *)
-let one = { poly = Node.M.empty }
+let one = { cst = A.one; poly = Node.M.empty }
 (*
  * let is_cst p = if Node.M.is_empty p.poly then Some p.cst else None
  * (\* let is_zero p = Q.equal p.cst Q.zero && Node.M.is_empty p.poly *\)
@@ -77,10 +76,10 @@ let extract p =
     (* else Cst p.cst *)
   else
     let x, q = Shuffle.chooseb Node.M.choose Node.M.choose_rnd p.poly in
-    let p' = { (* p with  *) poly = Node.M.remove x p.poly } in
+    let p' = { p with poly = Node.M.remove x p.poly } in
     Var (q, x, p')
 
-let remove n p = { poly = Node.M.remove n p.poly }
+let remove n p = { p with poly = Node.M.remove n p.poly }
 
 type kind = ONE | NODE of Node.t | PRODUCT
 
@@ -92,7 +91,7 @@ let classify p =
   else PRODUCT
 
 let monome x c =
-  if Q.equal Q.zero c then one else { poly = Node.M.singleton x c }
+  if Q.equal Q.zero c then one else { cst = A.one; poly = Node.M.singleton x c }
 
 let is_one_node p =
   (* cst = 0 and one empty monome *)
@@ -101,9 +100,9 @@ let is_one_node p =
     if Q.equal Q.one k then Some node else None
   else None
 
-(* let mult_cst c p1 =
- *   assert (not (Q.equal c Q.zero));
- *   { p1 with cst = Q.mul p1.cst c } *)
+let mult_cst c p1 =
+  assert (A.is_not_zero c);
+  { p1 with cst = A.mul p1.cst c }
 
 let power_cst c p1 =
   if Q.equal Q.zero c then one
@@ -111,7 +110,7 @@ let power_cst c p1 =
   else
     let poly_mult c m = Node.M.map (fun c1 -> Q.mul c c1) m in
     if Q.equal Q.zero c then one
-    else { (* cst = Q.mul c p1.cst; *) poly = poly_mult c p1.poly }
+    else { cst = A.positive_pow p1.cst c; poly = poly_mult c p1.poly }
 
 (** Warning Node.M.union can be used only for defining an operation [op]
     that verifies [op 0 p = p] and [op p 0 = p] *)
@@ -119,7 +118,7 @@ let mul p1 p2 =
   let poly_add m1 m2 =
     Node.M.union (fun _ c1 c2 -> Q.none_zero (Q.add c1 c2)) m1 m2
   in
-  { (* cst = Q.mul p1.cst p2.cst; *) poly = poly_add p1.poly p2.poly }
+  { cst = A.mul p1.cst p2.cst; poly = poly_add p1.poly p2.poly }
 
 let div p1 p2 =
   let poly_sub m1 m2 =
@@ -130,7 +129,7 @@ let div p1 p2 =
         | Some c1 -> Q.none_zero (Q.sub c1 c2))
       m1 m2
   in
-  { (* cst = Q.div p1.cst p2.cst; *) poly = poly_sub p1.poly p2.poly }
+  { cst = A.div p1.cst p2.cst; poly = poly_sub p1.poly p2.poly }
 
 let x_m_yc p1 p2 c =
   assert (not (Q.equal c Q.zero));
@@ -143,7 +142,7 @@ let x_m_yc p1 p2 c =
         | Some c1 -> Q.none_zero (f c1 c2))
       m1 m2
   in
-  { (* cst = Q.mul p1.cst (Q. p2.cst; *) poly = poly p1.poly p2.poly }
+  { cst = A.mul p1.cst (A.positive_pow p2.cst c); poly = poly p1.poly p2.poly }
 
 let xc_m_yc p1 c1 p2 c2 =
   let c1_iszero = Q.equal c1 Q.zero in
@@ -163,7 +162,10 @@ let xc_m_yc p1 c1 p2 c2 =
           | Some e1, Some e2 -> Q.none_zero (f e1 e2))
         m1 m2
     in
-    { (* cst = f p1.cst p2.cst; *) poly = poly p1.poly p2.poly }
+    {
+      cst = A.mul (A.positive_pow p1.cst c1) (A.positive_pow p2.cst c2);
+      poly = poly p1.poly p2.poly;
+    }
 
 let subst_node p x y =
   let poly, qo = Node.M.find_remove x p.poly in
@@ -175,29 +177,30 @@ let subst_node p x y =
           (function None -> qo | Some q' -> Q.none_zero (Q.add q q'))
           y poly
       in
-      ({ poly }, q)
+      ({ p with poly }, q)
 
 let subst p x s =
   let poly, q = Node.M.find_remove x p.poly in
-  match q with None -> (p, Q.zero) | Some q -> (x_m_yc { poly } s q, q)
+  match q with None -> (p, Q.zero) | Some q -> (x_m_yc { p with poly } s q, q)
 
 let fold f acc p = Node.M.fold_left f acc p.poly
 let iter f p = Node.M.iter f p.poly
 
-let of_list (* cst  *) l =
+let of_list cst l =
   let fold acc (node, q) =
     Node.M.change
       (function None -> Some q | Some q' -> Q.none_zero (Q.add q q'))
       node acc
   in
-  { (* cst; *) poly = List.fold_left fold Node.M.empty l }
+  { cst; poly = List.fold_left fold Node.M.empty l }
 
-let of_map m =
+let of_map cst m =
   assert (Node.M.for_all (fun _ v -> Q.sign v <> 0) m);
-  { poly = m }
+  { cst; poly = m }
 
 let common pa pb =
   {
+    cst = A.min pa.cst pb.cst;
     poly =
       Node.M.inter
         (fun _ a b -> Q.none_zero (Q.min (Q.floor a) (Q.floor b)))
diff --git a/src_colibri2/theories/LRA/product.mli b/src_colibri2/theories/LRA/product.mli
index 042123e4a42132a429177c31b97f5ba353e523e5..07d094a9cc38e0ad2673c42fe44ea6cf2bf889fa 100644
--- a/src_colibri2/theories/LRA/product.mli
+++ b/src_colibri2/theories/LRA/product.mli
@@ -23,7 +23,7 @@ open Colibri2_popop_lib
 open Colibri2_core
 (** Polynome *)
 
-type t = private { (* cst : Q.t;  *) poly : Q.t Node.M.t }
+type t = private { cst : A.t; poly : Q.t Node.M.t }
 
 include Popop_stdlib.Datatype with type t := t
 
@@ -43,7 +43,7 @@ type kind = ONE | NODE of Node.t | PRODUCT
 
 val classify : t -> kind
 val power_cst : Q.t -> t -> t
-val of_list : (Node.t * Q.t) list -> t
+val of_list : A.t -> (Node.t * Q.t) list -> t
 val mul : t -> t -> t
 val div : t -> t -> t
 val x_m_yc : t -> t -> Q.t -> t
@@ -52,7 +52,7 @@ val subst : t -> Node.t -> t -> t * Q.t
 val subst_node : t -> Node.t -> Node.t -> t * Q.t
 val fold : ('a -> Node.t -> Q.t -> 'a) -> 'a -> t -> 'a
 val iter : (Node.t -> Q.t -> unit) -> t -> unit
-val of_map : Q.t Node.M.t -> t
+val of_map : A.t -> Q.t Node.M.t -> t
 
 val common : t -> t -> t
 (** Return the common part *)