let rec zsolve (a, b) = 
  fresh := [];           (* reinitialize fresh variable generation. *)
  if Term.eq a b then [] else
    if Term.is_var a && Term.is_var b then
      [Term.orient(a, b)]
    else 
      let a_sub_b = mk_sub a b in
      let q = constant_of a_sub_b 
      and ml = nonconstant_monomials_of a_sub_b in     (* [q + ml = 0] *)
        match ml with
          | [] -> 
              if Q.is_zero q then [] else raise(Exc.Inconsistent)
          | [m] ->                           (* [q + p*x = 0] *)
              let p = coefficient_of_mono m
              and x = variable_of_mono m in
              let q_div_p = Q.div q p in
                if Q.is_integer q_div_p then 
                  [(x, mk_num (Q.minus q_div_p))]
                else 
                  raise Exc.Inconsistent
          | _ ->
              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 gl = general cl (d, pl) in
                         combine xl gl)
             
and vectorize ml =
  let rec loop (ql, xl) = function
    | [] -> 
        (ql, xl)
    | m :: ml -> 
        let q = coefficient_of_mono m
        and x = variable_of_mono m in
          loop (q :: ql, x :: xl) ml
  in
    loop ([], []) ml

and combine xl bl =
  List.combine xl bl




(** 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 rec loop al zl =
    match al, zl with
      | [_], [_] -> zl
      | a0 :: ((a1 :: al'') as al'),  z0 :: z1 :: zl'' ->
          let k = mk_fresh () in
          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
          let sl' =  loop al' (e1 :: zl'') in
            e0 :: sl'
      | _ -> assert false
  in
    loop al (List.map mk_num pl)