(*
 * Cooper's Quantifier Elimination Method for Presburger Arithmetic
 * ----------------------------------------------------------------
 * D. C. Cooper, "Theorem Proving in Arithmetic without Multiplication"
 *
 * Author: Aaron R. Bradley
 * 
 * Last modified:
 * -------------
 * 042406 -- final optimization & cleaned up
 * 042506 -- removed usemem, as unnecessary
 * 042606 -- improved divopt for alternation;
 *           fixed overflow bug in divfloor/ceil;
 *           modifed left/right heuristic to take #vars into account
 * 050306 -- changed band, bandl, bor, borl to sets;
 *           replaced int with Big_int.big_int (many changes)
 * 050406 -- changed qe to run depth-first; added iprog;
 *           improved performance of bnot
 * 031307 -- fixed to_unit: conjoin 1 | s when lcm = 1 in Step 4
 * 112707 -- fixed bug in disj;
 *           fixed case when triangulated constraints have True, False,
 *             or no constraint over one variable (only in alternations);
 *           fixed conflict in parsr.mly
 *
 * To do:
 * -----
 *
 * Copyright (c) 2006, Aaron R. Bradley
 *
 * <GNU public license>
 *)

module H = Hashtbl
module L = List

open Big_int

(* ---------- Options ---------- *)

(* true: use division optimization described in appendix of paper *)
let appendix = ref true
let verbose = ref false

(* for convenience *)
let soi x = string_of_big_int x
let pr x = if !verbose then print_endline x
let prs x = if !verbose then print_string x

let bi x = big_int_of_int x
let ne_big_int x y = not (eq_big_int x y)
let biz = bi 0
let bio = bi 1
let eqz x = eq_big_int x biz
let nez x = not (eqz x)

(* ---------- Expression type and methods ---------- *)

type expr =
  | Int of big_int                       (* constants *)
  | Var of big_int * string              (* variable w/ coefficient *)
  | Plus of expr * expr
  | GT of expr * expr                    (* predicate: greater-than *)
  | Div of big_int * expr                (* predicate: divides *)
  | NDiv of big_int * expr               (* predicate: not divides *)
  | True
  | False
  | And of expr * expr
  | Or of expr * expr
  | BigAnd of string * big_int * big_int * expr  
                                         (* v, l, u, e: A v \in [l, u]. e *)
  | BigOr of string * big_int * big_int * expr   
                                         (* v, l, u, e: E v \in [l, u]. e *)
  | Forall of string * expr
  | Exists of string * expr

(* in: e : expr
 * out: string
 * Converts e to string.
 *)
let rec to_string e = if !verbose then ts e else ""
and string_of_expr e = ts e
and ts = function
  | Int n -> soi n
  | Var (n, s) -> (if ne_big_int n bio then (soi n) ^ " " else "") ^ s
  | Plus (e1, e2) -> (ts e1) ^ " + " ^ (ts e2)
  | GT (e1, e2) -> (ts e1) ^ " > " ^ (ts e2)
  | Div (n, e) -> (soi n) ^ " | " ^ (ts e)
  | NDiv (n, e) -> "~(" ^ (soi n) ^ " | " ^ (ts e) ^ ")"
  | True -> "true"
  | False -> "false"
  | And (e1, e2) -> "(" ^ (ts e1) ^ " && " ^ (ts e2) ^ ")"
  | Or (e1, e2) -> "(" ^ (ts e1) ^ " || " ^ (ts e2) ^ ")"
  | BigAnd (s, l, u, e) -> 
      "(AND(" ^ s ^ ", " ^ (soi l) ^ ", " ^ (soi u) ^ "). " ^ (ts e) ^")"
  | BigOr (s, l, u, e) -> 
      "(OR(" ^ s ^ ", " ^ (soi l) ^ ", " ^ (soi u) ^ "). " ^ (ts e) ^ ")"
  | Forall (s, e) -> "A " ^ s ^ ". " ^ (ts e)
  | Exists (s, e) -> "E " ^ s ^ ". " ^ (ts e)

(* expr set:
 * Unfortunately, compare cannot be used on big_ints, so I need to
 * implement my own comparision.  Another option would have been to
 * convert to strings and compare.
 *)
module ES = Set.Make (
  struct
    type t = expr
    let compare e1 e2 = 
      let rec cmp e1 e2 =
        match e1, e2 with
        | Int n1, Int n2 -> compare_big_int n1 n2
        | Var (n1, s1), Var (n2, s2) -> 
            let c = compare s1 s2 in
            if c = 0 then  compare_big_int n1 n2 else c
        | Div (n1, e1), Div (n2, e2)
        | NDiv (n1, e1), NDiv (n2, e2) ->
            let c = compare_big_int n1 n2 in
            if c = 0 then cmp e1 e2 else c
        | BigAnd (s1, l1, u1, e1), BigAnd (s2, l2, u2, e2)
        | BigOr (s1, l1, u1, e1), BigOr (s2, l2, u2, e2) ->
            let c = compare s1 s2 in
            if c = 0 then
              let c = compare_big_int l1 l2 in
              if c = 0 then
                let c = compare_big_int u1 u2 in
                if c = 0 then cmp e1 e2
                else c
              else c
            else c
        | Plus (e1, e2), Plus (f1, f2)
        | GT (e1, e2), GT (f1, f2)
        | And (e1, e2), And (f1, f2)
        | Or (e1, e2), Or (f1, f2) ->
            let c = cmp e1 f1 in
            if c = 0 then cmp e2 f2 else c
        | Forall (s1, e1), Forall (s2, e2)
        | Exists (s1, e1), Exists (s2, e2) ->
            let c = compare s1 s2 in
            if c = 0 then cmp e1 e2 else c
        | True, True
        | False, False -> 0
        | x, y -> (* different heads *)
            let hdv = function
              | Int _ -> 0
              | Var _ -> 1
              | Plus _ -> 2
              | GT _ -> 3
              | Div _ -> 4
              | NDiv _ -> 5
              | True -> 6
              | False -> 7
              | And _ -> 8
              | Or _ -> 9
              | BigAnd _ -> 10
              | BigOr _ -> 11
              | Forall _ -> 12
              | Exists _ -> 13
	    in
	    compare (hdv x) (hdv y)
      in
      cmp e1 e2
  end)

(* ---------- Constructors and expr utils ---------- *)
(* Methods for constructing canonical expressions. *)

let num n = Int n

let var n s = if eqz n then num biz else Var (n, s)

(* string set *)
module SS = Set.Make (
  struct 
    type t = string
    let compare s1 s2 = compare s1 s2
  end)

(* in: expr
 * out: expr set
 * Returns variables of e.
 *)
let rec vars = function
  | Var (_, s) -> SS.singleton s
  | Plus (e1, e2)
  | GT (e1, e2)
  | And (e1, e2)
  | Or (e1, e2) -> SS.union (vars e1) (vars e2)
  | Div (_, e)
  | NDiv (_, e) -> vars e
  | BigAnd (s, _, _, e)
  | BigOr (s, _, _, e)
  | Forall (s, e)
  | Exists (s, e) -> SS.remove s (vars e)
  | _ -> SS.empty

let plus e1 e2 = Plus (e1, e2)

(* factor out plusl, borl, andl *)
let fnl b fn = function
  | [] -> b
  | h :: [] -> h
  | h :: t -> L.fold_left fn h t

(* apply plus to list *)
let plusl es = fnl (num biz) plus es

(* return terms of sum *)
let rec terms = function
  | Plus (e1, e2) -> (terms e1) @ (terms e2)
  | x -> [x]

(* scale expr by c *)
let rec scale c = function
  | Int n -> num (mult_big_int c n)
  | Var (n, s) -> var (mult_big_int c n) s
  | Plus (e1, e2) -> plus (scale c e1) (scale c e2)
  | x -> pr (to_string x); assert false

(* negate expr *)
let minus e = scale (bi (-1)) e

(* greatest common divisor *)
let gcd m n = gcd_big_int m n

(* gcd applied to list *)
let gcdl = function
  | [] -> bi 1
  | n :: [] -> if eqz n then bi 1 else abs_big_int n
  | h :: t -> 
      L.fold_left gcd (if eqz h then bi 1 else abs_big_int h) t

(* remove 0s before applying actual gcdl *)
let gcdl ns = gcdl (L.filter (fun n -> nez n) ns)

(* least common multiple *)
let lcm m n = abs_big_int (div_big_int (mult_big_int m n) (gcd m n))

(* lcm applied to list *)
let lcml = function
  | [] -> bi 1
  | n :: [] -> abs_big_int n
  | h :: t -> L.fold_left lcm (abs_big_int h) t

(* filter: remove 0s before applying actual lcml *)
let lcml ns = lcml (L.filter (fun n -> nez n) ns)

(* in: e : expr, scalar : bool
 * out: expr
 * Required: e is term (summation, constant, or variable)
 * Result: at most one occurrence per variable, one constant.
 * scalar canonical: divides by gcd of coeff if scalar is true.
 *)
let canonical e scalar =
  let vs = H.create 13 in
  let add v n =
    let curr = try H.find vs v with
    | Not_found -> bi 0
    in
    H.replace vs v (add_big_int curr n)
  in
  let const = ref biz in
  let addc n = const := add_big_int !const n in
  L.iter
    (function
      | Int n -> addc n
      | Var (n, s) -> add s n
      | x -> pr (to_string x); assert false)
    (terms e);
  let div = 
    if scalar then begin
      (* scalar canonical, so divide by gcd of coeffs *)
      let coeffs = ref [!const] in
      H.iter (fun _ n -> coeffs := n :: !coeffs) vs;
      gcdl !coeffs
    end
    else bi 1
  in
  let can = 
    ref (if eqz !const then ES.empty else  
    ES.singleton (num (div_big_int !const div)))
  in
  H.iter
    (fun v n -> if nez n then can := ES.add (var (div_big_int n div) v) !can)
    vs;
  plusl (ES.elements !can)

(* in: s' : string (symbolic var), expr
 * out: expr
 * Assumes expr is linear and canonical; returns coefficient of s.
 * Raises Not_found if s' not in expr.
 *)
let rec coeff s' = function
  | Int _ -> raise Not_found
  | Var (n, s) -> if s = s' then n else raise Not_found
  | Plus (e1, e2) -> begin
      try coeff s' e1 with
      | Not_found -> coeff s' e2
  end
  | x -> pr (to_string x); assert false

(* Predicate constructors *)

(* simplifies *)
let gt e1 e2 = 
  match canonical (plus e1 (minus e2)) true with
  | Int n -> if gt_big_int n biz then True else False
  | e -> GT (e, num biz)

(* simplifies *)
let div n e = 
  assert (gt_big_int n biz); 
  if eq_big_int n bio then True else
  match canonical e false with
  | Int m -> if eqz (mod_big_int m n) then True else False
  | e -> Div (n, e)

(* simplifies *)
let ndiv n e = 
  assert (gt_big_int n biz); 
  if eq_big_int n bio then False else
  match canonical e false with
  | Int m -> if eqz (mod_big_int m n) then False else True
  | e -> NDiv (n, e)

(* Boolean constructors *)

(* returns conjuncts of expr *)
let rec conj = function
  | And (e1, e2) -> ES.union (conj e1) (conj e2)
  | x -> ES.singleton x

let fbandl es = fnl True (fun e1 e2 -> And (e1, e2)) es

let bandl es =
  let conj = 
    L.fold_left (fun s1 s2 -> ES.union s1 s2) ES.empty (L.map conj es) 
  in
  if ES.exists (function False -> true | _ -> false) conj then False else
  let es = ES.elements (ES.filter (function True -> false | _ -> true) conj) in
  fbandl es

let band e1 e2 = bandl [e1; e2]

(* returns disjuncts of expr *)
let rec disj = function
  | Or (e1, e2) -> ES.union (disj e1) (disj e2)
  | x -> ES.singleton x

let fborl es = fnl False (fun e1 e2 -> Or (e1, e2)) es

let borl es =
  let disj = 
    L.fold_left (fun s1 s2 -> ES.union s1 s2) ES.empty (L.map disj es) 
  in
  if ES.exists (function True -> true | _ -> false) disj then True else
  let es = ES.elements (ES.filter (function False -> false | _ -> true) disj) in
  fborl es

let bor e1 e2 = borl [e1; e2]

(* in: sub : (string, expr) H.t, expr
 * out: expr
 * Variable substitution.
 *)
let rec vsub sub = function
  | Var (n, s) as x -> begin
      try scale n (H.find sub s) with Not_found -> x
  end
  | Plus (e1, e2) -> plus (vsub sub e1) (vsub sub e2)
  | GT (e1, e2) -> gt (vsub sub e1) (vsub sub e2)
  | Div (n, e) -> div n (vsub sub e)
  | NDiv (n, e) -> ndiv n (vsub sub e)
  | And (e1, e2) -> band (vsub sub e1) (vsub sub e2)
  | Or (e1, e2) -> bor (vsub sub e1) (vsub sub e2)
  (* respect scoping through scope function *)
  | BigAnd (s, l, u, e) ->
      scope sub s (fun sub -> BigAnd (s, l, u, vsub sub e))
  | BigOr (s, l, u, e) ->
      scope sub s (fun sub -> BigOr (s, l, u, vsub sub e))
  | Forall (s, e) ->
      scope sub s (fun sub -> Forall (s, vsub sub e))
  | Exists (s, e) ->
      scope sub s (fun sub -> Exists (s, vsub sub e))
  | x -> x
(* removes mapping for s; applies fn; restores mapping *)
and scope sub s fn =
  try
    let v = H.find sub s in
    H.remove sub s;
    let rv = fn sub in
    H.add sub s v;
    rv
  with
  | Not_found -> fn sub

(* Quantifier constructors *)

(* in: v : string, l, u : int, e : expr, 
 *     b : expr, fn : string -> int -> int -> expr -> expr
 * out: expr
 * Constructs BigAnd or BigOr.  b is base value if [l, u] is empty.
 * fn constructs final expression.
 *)
let bigX v l u e b fn =
  match e with
  | True -> True
  | False -> False
  | e ->
      let l, u =
        if ge_big_int l u then (l, u) else begin
          (* use simple atoms to reduce l, u, if possible *)
          let l = ref l in
          let u = ref u in
          ES.iter
            (function
              | GT (Plus (Var (m, v'), Int n), Int z)
              | GT (Plus (Int n, Var (m, v')), Int z) ->
                  if v = v' && eqz z then
                    if gt_big_int m biz && gt_big_int biz n then begin
                      let k = div_big_int (minus_big_int n) m in
                      if ge_big_int k !l then l := succ_big_int k
                    end
                    else if gt_big_int biz m && gt_big_int n biz then begin
                      let k = 
                        if eqz (mod_big_int n m) 
                        then div_big_int n (minus_big_int m) 
                        else succ_big_int (div_big_int n (minus_big_int m))
                      in
                      if le_big_int k !u then u := pred_big_int k
                    end
              | GT (Var (n, v'), Int z) when 
                  eq_big_int n (minus_big_int bio) && eqz z && v = v' -> 
                    u := biz
              | _ -> ())
            (conj e);
          (!l, !u)
        end
      in
      (* construct the BigX expr *)
      if lt_big_int l u then fn v l u e
      else if eq_big_int l u then begin
        let sub = H.create 1 in
        H.add sub v (num l);
        vsub sub e
      end
      else b

let bigand v l u e = bigX v l u e True (fun v l u e -> BigAnd (v, l, u, e))

let bigor v l u e = bigX v l u e False (fun v l u e -> BigOr (v, l, u, e))

let forall v e = Forall (v, e)

let exists v e = Exists (v, e)

(* Derived constructors *)

let lt e1 e2 = gt e2 e1

let ge e1 e2 = gt (plus e1 (num bio)) e2

let le e1 e2 = ge e2 e1

let eq e1 e2 = band (le e1 e2) (ge e1 e2)

let ne e1 e2 = bor (lt e1 e2) (gt e1 e2)

(* needs to be really efficient, so rely on canonical form when possible *)
let rec bnot = function
  | GT (e, Int z) ->
      assert (eqz z);
      let hasc = ref false in
      let terms =
	L.map
	  (function
	    | Int n -> hasc := true; Int (sub_big_int bio n)
	    | Var (n, s) -> Var (minus_big_int n, s)
	    | _ -> assert false)
	  (terms e)
      in
      GT (plusl (if !hasc then terms else (num bio) :: terms), Int z)
  | Div (n, e) -> NDiv (n, e)
  | NDiv (n, e) -> Div (n, e)
  | True -> False
  | False -> True
  | And (e1, e2) -> Or (bnot e1, bnot e2)
  | Or (e1, e2) -> And (bnot e1, bnot e2)
  | BigAnd (v, l, u, e) -> BigOr (v, l, u, (bnot e))
  | BigOr (v, l, u, e) -> BigAnd (v, l, u, (bnot e))
  | Forall (v, e) -> Exists (v, bnot e)
  | Exists (v, e) -> Forall (v, bnot e)
  | x -> pr (to_string x); assert false

(* ---------- Cooper functions ---------- *)

(* set of integers *)
module IS = Set.Make (
  struct
    type t = big_int
    let compare e1 e2 = compare_big_int e1 e2
  end)

(* in: s : string (var), e : expr
 * out: int
 * Computes lcm of coefficients of s, for Step 4 
 *)
let lcm_coeff s e =
  let rec coeff = function
    | Var (n, s') -> if s = s' then IS.singleton n else IS.empty
    | Plus (e1, e2)
    | GT (e1, e2)
    | And (e1, e2)
    | Or (e1, e2) -> IS.union (coeff e1) (coeff e2)
    | Div (_, e)
    | NDiv (_, e)
    | BigAnd (_, _, _, e)
    | BigOr (_, _, _, e) -> coeff e
    | Forall _ | Exists _ -> assert false
    | _ -> IS.empty
  in
  lcml (IS.elements (coeff e))

(* in: s : string (var), e : expr
 * out: int
 * Computes lcm of n in n | x + ... of e, for Step 5 
 *)
let lcm_div s e =
  let rec div = function
    | True | False | GT _ -> IS.empty
    | Div (n, e)
    | NDiv (n, e) -> begin
        try ignore (coeff s e); IS.singleton n with Not_found -> IS.empty
    end
    | And (e1, e2)
    | Or (e1, e2) -> IS.union (div e1) (div e2)
    | BigAnd (_, _, _, e)
    | BigOr (_, _, _, e) -> div e
    | _ -> assert false
  in
  lcml (IS.elements (div e))

(* in: s : string, e : expr
 * out: expr
 * Step 4: Convert occurrences of s to have unit coefficients.
 * Assumes e is quantifer-free formula, linear terms are canonical.
 *)
let to_unit s e =
  (* delta' in notes *)
  let lcm = lcm_coeff s e in
  pr ("      [lcm: " ^ (soi lcm) ^ "]");
  (* check if already unit (including no occurrences of s) *)
  if eq_big_int lcm bio then
      (* Bug fix 3/13/2007: conjoin 1 | s when lcm = 1 *)
      band e (div lcm (var bio s))
  else
  (* absolute value of coeff of s *)
  let acoeff e = abs_big_int (coeff s e) in
  (* s in linear expr to unit *)
  let rec ltu = function
    | Int _ as x -> x
    | Var (n, s') as x -> 
        if s = s' then begin
          assert (eq_big_int (abs_big_int n) lcm);
          Var (div_big_int n lcm, s)
        end
        else x
    | Plus (e1, e2) -> Plus (ltu e1, ltu e2)
    | x -> pr (to_string x); assert false
  in
  (* s in expr to unit *)
  let rec tu = function
    | GT (e, Int z) as x when eqz z -> begin
        try 
          let c = div_big_int lcm (acoeff e) in
          GT (ltu (scale c e), Int biz) 
        with
        | Not_found -> x
    end
    | Div (n, e) as x -> begin
        try
          let c = div_big_int lcm (acoeff e) in
          Div (mult_big_int c n, ltu (scale c e))
        with
        | Not_found -> x
    end
    | NDiv (n, e) as x -> begin
        try
          let c = div_big_int lcm (acoeff e) in
          NDiv (mult_big_int c n, ltu (scale c e))
        with
        | Not_found -> x
    end
    | True -> True
    | False -> False
    | And (e1, e2) -> And (tu e1, tu e2)
    | Or (e1, e2) -> Or (tu e1, tu e2)
    | BigAnd (s', l, u, e) as x -> 
        if s <> s' then BigAnd (s', l, u, tu e) else x
    | BigOr (s', l, u, e) as x -> 
        if s <> s' then BigOr (s', l, u, tu e) else x
    | x -> pr (to_string x); assert false
  in
  band (tu e) (div lcm (var bio s))

(* named variable constructor *)
let vi = ref 0
let fresh () = 
  incr vi;
  "j" ^ (string_of_int !vi)

(* for Step 4 & 5: left version or right version *)
type dir =
  | Left
  | Right

(* in: s : string (var), e : expr, dir : dir
 * out: expr
 * Step 5: Computes unbounded disjunct; dir indicates type.
 * Assumes e is quantifer-free, linear terms are canonical.
 *)
let rec inf s e dir =
  let ntrue, nfalse = 
    match dir with
    | Left -> True, False
    | Right -> False, True
  in
  let rec inf = function
    | GT (e, Int z) as x when eqz z -> begin
        try
          let c = coeff s e in
          assert (nez c);
          if lt_big_int c biz then ntrue else nfalse
        with
        | Not_found -> x
    end
    | And (e1, e2) -> band (inf e1) (inf e2)
    | Or (e1, e2) -> bor (inf e1) (inf e2)
    | BigAnd (s', l, u, e) as x -> 
        if s <> s' then BigAnd (s', l, u, inf e) else x
    | BigOr (s', l, u, e) as x -> 
        if s <> s' then BigOr (s', l, u, inf e) else x
    | Forall _ | Exists _ -> assert false
    | x -> x
  in
  let lcm = lcm_div s e in
  let e' = inf e in
  let j = fresh () in
  let sub = H.create 1 in
  let jv = var bio j in
  H.add sub s (if dir = Left then jv else minus jv);
  bigor j bio lcm (vsub sub e')

(* in: s : string (var), e : expr, dir : dir
 * out: expr set
 * Computes set of base terms:
 *   c s + t > 0  ==> -t if left
 *   -c s + t > 0 ==> t if right
 *)
let bases s e dir =
  let rec bases = function
    | GT (e, Int z) when eqz z -> begin
        try
          let c = coeff s e in
          assert (nez c);
          (* if c * s > b then {b} OR if c * s < a then {a} *)
          if gt_big_int c biz && dir = Left
          then ES.singleton (canonical (plus (var c s) (minus e)) false)
          else if lt_big_int c biz && dir = Right
          then ES.singleton (canonical (plus (var (minus_big_int c) s) e) false)
          else ES.empty
        with
        | Not_found -> ES.empty
    end
    | True | False | Div _ | NDiv _ -> ES.empty
    | And (e1, e2)
    | Or (e1, e2) -> ES.union (bases e1) (bases e2)
    | BigAnd (s', _, _, e)
    | BigOr (s', _, _, e) -> if s <> s' then bases e else ES.empty
    | x -> pr (to_string x); assert false
  in
  bases e

(* in: s : string, e : expr, b : expr, dir : dir
 * out: expr
 * Step 5: Computes bounded disjunct.
 * Assumes e is quantifer-free formula, linear terms are canonical,
 * each occurrence of s has unit coefficient.
 *)
let base s e b dir =
  let j = fresh () in
  (* delta in notes *)
  let lcm = lcm_div s e in
  pr ("      [lcm: " ^ (soi lcm) ^ "]");
  let sub = H.create 1 in
  pr ("      [base: " ^ (to_string b) ^ "]");
  let jv = var bio j in
  H.add sub s (plus b (if dir = Left then jv else minus jv));
  bigor j bio lcm (vsub sub e)

(* constructs one-time iterator *)
let once e =
  let once = ref 0 in
  let fn () =
    incr once;
    if !once = 1 then Some e else None
  in
  fn

(* in: s : string, e : expr
 * out: expr
 * The Big Event.  Eliminate s from e.
 * Assumes e is quantifer-free formula, linear terms are canonical.
 *)
let eid = ref 0  (* unique id *)
let eliminate s e =
  if not (SS.mem s (vars e)) then once e else begin
    incr eid;
    let sid = "(" ^ (string_of_int !eid) ^ ")" in
    pr (sid ^ " Eliminating " ^ s ^ " from: " ^ (to_string e));
    let ue = to_unit s e in
    pr (sid ^ "  Step 4: " ^ (to_string ue));
    (* optimization: choose direction based on #vars in base sets *)
    let lbases = bases s ue Left in
    let rbases = bases s ue Right in
    let dir, bases =
      let bvars bases =
        let bvars = ref SS.empty in
        ES.iter (fun e -> bvars := SS.union (vars e) !bvars) bases;
        !bvars
      in
      let lvars, rvars = bvars lbases, bvars rbases in
      if ((SS.cardinal lvars < SS.cardinal rvars) 
            or (SS.cardinal lvars = SS.cardinal rvars 
                  && ES.cardinal lbases < ES.cardinal rbases))
      then (Left, lbases) 
      else (Right, rbases)
    in
    pr (sid ^ "    Direction: " ^ (if dir = Left then "left" else "right"));
    let state = ref 0 in
    let bases = ref bases in
    let fn () =
      incr state;
      match !state with
      | 1 ->
	  let inf = inf s ue dir in
	  pr (sid ^ "  Step 5 (inf): " ^ (to_string inf));
	  Some inf
      | _ when ES.is_empty !bases -> None
      | _ ->
	  let bf = ES.choose !bases in
	  bases := ES.remove bf !bases;
	  let base = base s ue bf dir in
	  pr (sid ^ "  Step 5 (base): " ^ (to_string base));
	  Some base
    in
    fn
  end

exception AllTrue

(* in: e : expr
 * out: expr
 * Instantiate all bounded quantifiers (hence the name "blowup").  This 
 * implementation is the naive one.  Enumerate possibilities, stopping if 
 * an expression simplifies to True.
 *)
let blowup e =
  pr ("Expanding: " ^ (to_string e));
  let count = ref 0 in
  let incr () = 
    incr count;
    if !count > 100000 then begin
      (* too many subexpressions, so abort *)
      pr "Are you crazy?";
      assert false
    end
  in
  (* should not contain BigAnd or BigOr within conjunction *)
  let rec bu = function
    | Or (e1, e2) -> bor (bu e1) (bu e2)
    | BigOr (v, l, u, e) -> 
        let e' = ref False in
        let sub = H.create 1 in
        let i = ref l in
        while le_big_int l u do
          H.replace sub v (num !i);
          match vsub sub e with
          | True -> raise AllTrue
          | False -> i := succ_big_int !i
          | x -> 
              let e = bu x in
              incr (); 
              e' := bor !e' e;
              i := succ_big_int !i
        done;
        !e'
    | BigAnd _ -> assert false
    | x -> x
  in
  try bu e with
  | AllTrue -> True

(* Functions for divides optimizations *)

(* O'Caml's division and modulo semantics leave a bit to be desired. *)

let mdiv a b =
  let a' = abs_big_int a in
  let b' = abs_big_int b in
  let s = 
    if ((ge_big_int a biz && gt_big_int b biz)
      || (le_big_int a biz && lt_big_int b biz))
    then 1 else -1
  in
  let d = div_big_int a' b' in
  if eqz (mod_big_int a' b') then mult_big_int (bi s) d 
  else if s > 0 then d 
  else (pred_big_int (minus_big_int d))

let divceil a b = 
  if eqz (mod_big_int a b) then mdiv a b else succ_big_int (mdiv a b)

let divfloor a b = mdiv a b

let extended_gcd a b =
  let rec mmod a b = 
    if gt_big_int a biz && gt_big_int b biz 
    then mod_big_int a b
    else if lt_big_int a biz && lt_big_int b biz 
    then minus_big_int (mod_big_int (abs_big_int a) (abs_big_int b)) 
    else mmod (add_big_int a b) b
  in
  let rec egcd a b =
    if eqz (mmod a b) then 
      (if gt_big_int b biz then (biz, bio) 
      else (biz, minus_big_int bio)) 
    else
    let x, y = egcd b (mmod a b) in
    (y, sub_big_int x (mult_big_int y (mdiv a b)))
  in
  if eqz a || eqz b then (bio, biz, biz) else
  let p, q = egcd a b in
  let d = add_big_int (mult_big_int a p) (mult_big_int b q) in
  (d, p, q)  

exception NoSolution

(* in: domain : (string, (int, int)) list, divs : expr set, 
       cb : (string, expr) H.t -> expr
 * out: unit
 * Solve divs constraints within domain.  domain is a list of variables and
 * their interval domains.  cb is a callback function to be called with 
 * solutions.
 *)
let solve domain divs cb =
  (* choose subset divs' of divs s.t. |divs'| > 1 and all of divs' contain
   * "largest" var in order defined by given domain list
   *)
  let rec subdivs divs = function
    | [] -> None
    | (v, _) :: t as x ->
        (* set of divs with only v, t vars *)
        let rvars = ref SS.empty in
        L.iter (fun (v, _) -> rvars := SS.add v !rvars) x;
        let fdivs = 
          ES.filter (fun d -> SS.is_empty (SS.diff (vars d) !rvars)) divs
        in
        (* set of remaining divs with v, and rest *)
        let dv, _ = ES.partition (fun d -> SS.mem v (vars d)) fdivs in
        if ES.cardinal dv > 1 then Some (v, dv, ES.diff divs dv)
        else subdivs divs t
  in
  (* Theorem 1 of appendix *)
  let appThI x d1 d2 =
    match d1, d2 with
    | Div (m, e1), Div (n, e2) -> 
        (* extract m | ax + b, n | alpha x + beta *)
        let a = coeff x e1 in
        let b = canonical (plus e1 (minus (var a x))) false in
        let alpha = coeff x e2 in
        let beta = canonical (plus e2 (minus (var alpha x))) false in
        (* apply Theorem 1 of appendix of Cooper's paper *)
        let d, p, q = extended_gcd (mult_big_int a n) (mult_big_int alpha m) in
        (div (mult_big_int m n) 
           (plusl [var d x; scale (mult_big_int p n) b; 
                   scale (mult_big_int q m) beta]),
         div d (plus (scale alpha b) (minus (scale a beta))))
    | True, x | x, True -> (x, True)
    | False, _ | _, False -> raise NoSolution
    | _ -> assert false
  in
  (* return equivalent divs' s.t. only one div contains v *)
  let rec reduce v divs =
    if ES.cardinal divs <= 1 then divs else
    let d1 = ES.choose divs in
    let d2 = ES.choose (ES.remove d1 divs) in
    let r1, r2 = appThI v d1 d2 in
    let r = reduce v (ES.add r1 (ES.remove d2 (ES.remove d1 divs))) in
    match r2 with
    | True -> r
    | r2 -> ES.add r2 r
  in
  (* return equivalent divs' in triangular form, where simplest div
   * contains only last variable of domain *)
  let rec triangulate divs =
    match subdivs divs domain with
    | None -> divs
    | Some (v, dv, dr) -> triangulate (ES.union (reduce v dv) dr)
  in
  (* Theorem 2 of appendix *)
  let appThII = function
    | Div (m, e) ->
        (* extract m | ax + b, for a, b constants *)
        let vars = vars e in
        assert (SS.cardinal vars = 1);
        let x = SS.choose vars in
        let a = coeff x e in
        let b = 
          match canonical (plus e (minus (var a x))) false with
          | Int n -> n
          | _ -> assert false
        in
        (* apply Theorem II of appendix of Cooper's paper *)
        let d, p, q = extended_gcd a m in
        if eqz (mod_big_int b d) then begin
            let l, u = L.assoc x domain in
            (* divceil (p * (b / d) + l) (m / d) *)
            let l' = divceil 
                (succ_big_int (mult_big_int p (div_big_int b d))) 
                (div_big_int m d) 
            in
            (* divfloor (p * (b / d) + u) (m / d) *)
            let u' = divfloor
                (add_big_int (mult_big_int p (div_big_int b d)) u)
                (div_big_int m d)
            in
            pr ("t (" ^ x ^ ") in [" ^ (soi l') ^ ", " ^ (soi u') ^ "]");
            let sol = ref [] in
            let t = ref l' in
            while le_big_int !t u' do
              (* (-p * (b / d) + t * (m / d)) *)
              let s = add_big_int
                  (mult_big_int (minus_big_int p) (div_big_int b d))
                  (mult_big_int !t (div_big_int m d))
              in
              sol := s :: !sol;
              t := succ_big_int !t
            done;
            (x, L.rev !sol)
        end
        else (x, [])
    | _ -> assert false
  in
  (* Given divs in triangular form, back-substitute to solve.  Call cb
   * with all generated solutions. *)
  let rec backsub soln divs =
    if ES.mem False divs then () else
    let divs = ES.filter (function True -> false | _ -> true) divs in
    (* Generalize technique for quantifier-alternation case.  If not
     * all (bounded quantified) variables have an assignment, try all 
     * possible completions.  Skip dummies (representing free vars).
     *)
    let rec tryall = function
      | [] -> cb soln
      | (v, (l, u)) :: t ->
          if H.mem soln v then tryall t else
          (* skip dummies *)
          if eq_big_int l u && eqz l then tryall t else
          let i = ref l in
          while le_big_int !i u do
            H.add soln v (num !i);
            tryall t;
            H.remove soln v;
            i := succ_big_int !i
          done
    in
    if ES.cardinal divs = 0 then tryall domain else
    (* s1 is set of n | a x + b, for b a constant *)
    let s1, s2 = ES.partition 
        (function
          | Div (n, e) -> SS.cardinal (vars e) = 1
          | _ -> assert false)
        divs
    in
    if ES.cardinal s1 = 0 then tryall domain else
    let d = ES.choose s1 in
    let v, sol = appThII d in
    assert (not (H.mem soln v));
    (* recurse on all solutions to v *)
    L.iter
      (fun sol ->
        H.add soln v (Int sol);
        let divs' = ref ES.empty in
        ES.iter (fun d -> divs' := ES.add (vsub soln d) !divs') divs;
        backsub soln !divs';
        H.remove soln v)
      sol
  in
  pr ("In context: "
      ^ String.concat ", " 
          (L.map 
             (fun (v, (l, u)) -> (v ^ " in [" ^ (soi l) ^ ", " ^ (soi u) ^ "]"))
             domain));
  ES.iter (fun d -> pr (to_string d)) divs;
  pr "Triangulating...";
  let tri = triangulate divs in
  if ES.exists (function False -> true | _ -> false) tri
  then
    pr "No solution"
  else begin
    ES.iter (fun d -> pr (to_string d)) tri;
    (* now remove tri'ed divs with free vars *)
    let dummies = ref SS.empty in
    L.iter 
      (fun (v, (l, u)) -> 
         if eq_big_int l u && eqz l then dummies := SS.add v !dummies)
      domain;
    let ftri = ES.filter 
        (function
         | Div _ as x -> SS.is_empty (SS.inter (vars x) !dummies)
         | True -> false
         | _ -> assert false)
        tri
    in
    pr "Filtered:";
    ES.iter (fun d -> pr (to_string d)) ftri;
    try backsub (H.create 13) ftri with
    | NoSolution -> ()
  end

(* in: e : expr
 * out: expr
 * Instantiates bounded quantifiers the smart way, by solving divides
 * constraints and instantiating.
 *
 * For quantifier alternation, does the best it can:
 *  1. Triangulates all divs constraints, handling free vars first.
 *  2. Filters out triangularized divs constraints with free variables.
 *     Handling free vars first leaves as many constraints remaining
 *     as possible.
 *  3. Solve remaining constraints.
 *)
let opt_blowup e =
  pr ("Expanding (opt): " ^ (to_string e));
  let count = ref 0 in
  let incr () = 
    incr count;
    if !count > 1000000 then begin
      pr "Are you crazy?";
      assert false
    end
  in
  (* context : (string, (int, int)) list (var, (lower, upper)) *)
  let expand context e =
    match context with
    | [] -> e    (* base case: just return e *)
    | context -> begin
        (* Generalize to quantifier alternation: add dummy contexts
         * for free variables.  Put them at the front so that solve
         * eliminates them first.
         *)
        let allvars = ref SS.empty in
        let divs = 
          ES.filter
            (function
              | Div _ as x -> allvars := SS.union (vars x) !allvars; true
              | _ -> false)
            (conj e)
        in
        (* find free vars *)
        let cvars = ref SS.empty in
        L.iter (fun (v, _) -> cvars := SS.add v !cvars) context;
        let fvars = SS.diff !allvars !cvars in
        (* construct dummy contexts *)
        let fcontext = ref [] in
        SS.iter (fun v -> fcontext := (v, (biz, biz)) :: !fcontext) fvars;
        (* new context *)
        let context = !fcontext @ context in
        (* build formula on e' *)
        let e' = ref False in
        (* callback: instantiate e under sub *)
        let cb sub = 
          let se = vsub sub e in
          match se with
          | True -> 
              H.iter 
                (fun k v -> prs (k ^ ":" ^ (to_string v) ^ " ")) 
                sub;
              pr "";
              raise AllTrue
          | False -> ()
          | x ->
              incr ();
              e' := bor !e' se
        in
        solve context divs cb;
        !e'
    end
  in
  (* should not contain BigOr or BigAnd within disjunction *)
  let rec bu context = function
    | Or (e1, e2) as x -> 
        bor (bu context e1) (bu context e2)
    | BigOr (v, l, u, e) -> 
        bu ((v, (l, u)) :: context) e
    | BigAnd _ -> assert false
    | x -> expand context x
  in
  try bu [] e with
  | AllTrue -> True

(* ---------- Main functions ---------- *)

let blowup e = if !appendix then opt_blowup e else blowup e

(* iterate over it, building disjunction *)
let idisj it =
  let disj = ref [] in
  let fin = ref false in
  while not !fin do
    match it () with
    | None -> fin := true
    | Some True -> disj := [True]; fin := true
    | Some e -> disj := (blowup e) :: !disj
  done;
  borl !disj

(* delve into BigOr's *)
let rec elim v = function
  | BigOr (j, l, u, e) -> 
      let it = elim v e in
      let fn () = 
	match it () with
	| None -> None
	| Some e -> Some (bigor j l u e)
      in
      fn
  | x -> eliminate v x

(* for get_min *)
exception NotIProg

type bound =
  | Unbounded
  | Bound of big_int

let bmin m1 m2 =
  match m1, m2 with
  | Unbounded, _ | _, Unbounded -> Unbounded
  | Bound n1, Bound n2 -> Bound (min_big_int n1 n2)

let bmax m1 m2 =
  match m1, m2 with
  | Unbounded, _ | _, Unbounded -> Unbounded
  | Bound n1, Bound n2 -> Bound (max_big_int n1 n2)

type ival =
  | Ival of bound * bound
  | Empty

let inter i1 i2 =
  match i1, i2 with
  | Empty, _ | _, Empty -> Empty
  | Ival (min1, max1), Ival (min2, max2) ->
      let min = bmax min1 min2 in
      let max = bmin max1 max2 in
      let empty =
	match min, max with
	| Bound min, Bound max -> gt_big_int min max
	| _ -> false 
      in
      if empty then Empty else Ival (min, max)

let union i1 i2 =
  match i1, i2 with
  | Empty, x | x, Empty -> x
  | Ival (min1, max1), Ival (min2, max2) ->
      let min = bmin min1 min2 in
      let max = bmax max1 max2 in
      Ival (min, max)

(* in: e : expr
 * out: big_int * expr
 * Retrieves minimum from expression over one var.
 * Raises Not_found if unbounded.
 * TODO: 
 *  1. properly handle Div and NDiv
 *  2. SOUND???
 *)
let get_min e =
  let vars = vars e in
  if SS.cardinal vars > 1 then raise NotIProg;
  let mvar = SS.choose vars in
  let rec minmax = function
    | GT (e, Int z) when eqz z -> begin
	let c = coeff mvar e in
	match e with
	| Plus (Int n, Var _)
	| Plus (Var _, Int n) -> 
	    if lt_big_int c biz 
	    then Ival (Unbounded, 
		       Bound (pred_big_int (divceil n (minus_big_int c))))
	    else Ival (Bound (succ_big_int (divfloor (minus_big_int n) c)), 
		       Unbounded)
	| x -> pr (to_string x); assert false
    end
    | Div _ | NDiv _ -> Ival (Unbounded, Unbounded)
    | And (e1, e2) -> inter (minmax e1) (minmax e2)
    | Or (e1, e2) -> union (minmax e1) (minmax e2)
    | x -> pr (to_string x); assert false
  in
  match minmax e with
  | Empty -> (biz, False)
  | Ival (Unbounded, _) -> raise Not_found
  | Ival (Bound n, _) -> (n, ge (var bio mvar) (num n))

(* in: expr
 * out: expr iterator (unit -> e)
 * quantifier elimination: use iterators to achieve 
 * depth-first expansion 
 *)
let rec qe = function
  | Forall _ as x -> once (bnot (idisj (qe (bnot x))))
  | Exists (v, e) -> qe_exists v e
  | x -> once x
(* integer programming *)
and iprog = function
  | Forall _ as x -> begin
      let cmax = ref None in
      let fmax = ref True in
      let cmax (n, f) =
	if f = False then () else
	match !cmax with
	| None -> 
	    print_endline ("Current max: " ^ (soi n));
	    cmax := Some n;
	    fmax := f
	| Some m -> 
	    if lt_big_int m n then begin
	      cmax := Some n;
	      fmax := f;
	      print_endline ("Current max: " ^ (soi n));
	    end
      in
      let it = qe (bnot x) in
      let fin = ref false in
      while not !fin do
	match it () with
	| None -> fin := true
	| Some True -> fmax := False; fin := true
	| Some e -> begin
	    try cmax (get_min (bnot (blowup e))) with
	    | Not_found -> ()
	end
      done;
      once !fmax
  end
  | Exists (v, e) -> qe_exists v e
  | x -> once x
(* exists case *)
and qe_exists v e =
  let qiter = qe e in      (* iterator over recursive qe *)
  let eiter = ref None in  (* iterator over eliminate *)
  let fn () =
    let enext () =         (* return next eliminate disjunct *)
      match !eiter with
      | None -> assert false
      | Some it -> begin
	  match it () with
	  | None -> eiter := None; None
	  | x -> x
      end
    in
    let next () =          (* return next disjunct *)
      match !eiter with
      | None -> begin      (* move onto next qe disjunct *)
	  match qiter () with
	  | None -> None
	  | Some e -> 
	      eiter := Some (elim v e);
	      enext ()
      end
      | Some it -> enext ()
    in
    match next () with
    | None -> next ()      (* transition to next qe disjunct *)
    | x -> x
  in
  fn                       (* return iterator *)

(* in: e : expr
 * out: expr
 * Cooper's elimination procedure.  Do block elimination.
 *)
let cooper e = idisj (qe e)

let iprog e = 
  match (iprog e) () with
  | None -> assert false
  | Some e -> e