let rec zsolve (a, b) =    
  if Term.eq a b then [] else
    if is_var a && is_var b then
      [Term.orient(a, b)]
    else 
      let (q, ml) = poly_of (mk_sub a b) in   (* [q + ml = 0] *)
        if ml = [] then
          if Q.is_zero q then [] else raise(Exc.Inconsistent)
        else
          let (cl, xl) = vectorize ml in     (* [cl * xl = ml] in vector notation *)
            match Euclid.solve cl (Q.minus q) with
              | None -> raise Exc.Inconsistent
              | Some(d, pl) -> 
                  let (kl, gl) = general cl (d, pl) in
                    List.combine xl gl
             
and vectorize ml =
  let rec loop (ql, xl) = function
    | [] -> 
        (List.rev ql, List.rev xl) 
    | m :: ml ->
        let (q, x) = mono_of m in
          loop (q :: ql, x :: xl) ml
  in
    loop ([], []) ml




(** Compute the general solution of a linear Diophantine equation with coefficients al, the gcd d of al and a particular solution pl. In the case of four coeffients, compute, for example, (p0 p1 p2 p3) + k0/d * (a1 -a0 0 0) + k1/d * (0 a2 -a1 0) + k2/d * (0 0 a3 -a2) Here, k0, k1, and k2 are fresh variables. Note that any basis of the vector space of solutions xl of the equation al * xl = 0 would be appropriate. *)


and general al (d, pl) =
  let fl = ref [] in
  let rec loop al zl =
    match al, zl with
      | [_], [_] -> zl
      | a0 :: ((a1 :: al'') as al'),  z0 :: z1 :: zl'' ->
          let k = mk_fresh () in
            fl := k :: !fl;
            let e0 = mk_add z0 (mk_multq (Q.div a1 d) k) in
            let e1 = mk_add z1 (mk_multq (Q.div (Q.minus a0) d) k) in
              e0 :: loop al' (e1 :: zl'')
      | _ -> assert false
  in
    (!fl, loop al (List.map mk_num pl))