aboutsummaryrefslogtreecommitdiffhomepage
path: root/plugins/omega/omega.ml
blob: 2510c1693401ca697f4736868cabe67816362c09 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
(************************************************************************)
(*         *   The Coq Proof Assistant / The Coq Development Team       *)
(*  v      *   INRIA, CNRS and contributors - Copyright 1999-2018       *)
(* <O___,, *       (see CREDITS file for the list of authors)           *)
(*   \VV/  **************************************************************)
(*    //   *    This file is distributed under the terms of the         *)
(*         *     GNU Lesser General Public License Version 2.1          *)
(*         *     (see LICENSE file for the text of the license)         *)
(************************************************************************)
(**************************************************************************)
(*                                                                        *)
(* Omega: a solver of quantifier-free problems in Presburger Arithmetic   *)
(*                                                                        *)
(* Pierre Crégut (CNET, Lannion, France)                                  *)
(*                                                                        *)
(* 13/10/2002 : modified to cope with an external numbering of equations  *)
(*   and hypothesis. Its use for Omega is not more complex and it makes   *)
(*   things much simpler for the reflexive version where we should limit  *)
(*   the number of source of numbering.                                   *)
(**************************************************************************)

module type INT = sig
  type bigint
  val equal : bigint -> bigint -> bool
  val less_than : bigint -> bigint -> bool
  val add : bigint -> bigint -> bigint
  val sub : bigint -> bigint -> bigint
  val mult : bigint -> bigint -> bigint
  val euclid : bigint -> bigint -> bigint * bigint
  val neg : bigint -> bigint
  val zero : bigint
  val one : bigint
  val to_string : bigint -> string
end

let debug = ref false

module MakeOmegaSolver (I:INT) = struct

type bigint = I.bigint
let (=?) = I.equal
let (<?) = I.less_than
let (<=?) x y = I.less_than x y || x = y
let (>?) x y = I.less_than y x
let (>=?) x y = I.less_than y x || x = y
let (+) = I.add
let (-) = I.sub
let ( * ) = I.mult
let (/) x y = fst (I.euclid x y)
let (mod) x y = snd (I.euclid x y)
let zero = I.zero
let one = I.one
let two = one + one
let negone = I.neg one
let abs x = if I.less_than x zero then I.neg x else x
let string_of_bigint = I.to_string
let neg = I.neg

(* To ensure that polymorphic (<) is not used mistakenly on big integers *)
(* Warning: do not use (=) either on big int *)
let (<) = ((<) : int -> int -> bool)
let (>) = ((>) : int -> int -> bool)
let (<=) = ((<=) : int -> int -> bool)
let (>=) = ((>=) : int -> int -> bool)

let pp i = print_int i; print_newline (); flush stdout

let push v l = l := v :: !l

let rec pgcd x y = if y =? zero then x else pgcd y (x mod y)

let pgcd_l = function
  | [] -> failwith "pgcd_l"
  | x :: l -> List.fold_left pgcd x l

let floor_div a b =
  match a >=? zero , b >? zero with
    | true,true -> a / b
    | false,false -> a / b
    | true, false -> (a-one) / b - one
    | false,true  -> (a+one) / b - one

type coeff = {c: bigint ; v: int}

type linear = coeff list

type eqn_kind = EQUA | INEQ | DISE

type afine = {
  (* a number uniquely identifying the equation *)
  id: int ;
  (* a boolean true for an eq, false for an ineq (Sigma a_i x_i >= 0) *)
  kind: eqn_kind;
  (* the variables and their coefficient *)
  body: coeff list;
  (* a constant *)
  constant: bigint }

type state_action = {
  st_new_eq : afine;
  st_def    : afine; (* /!\ this represents [st_def = st_var] *)
  st_orig   : afine;
  st_coef   : bigint;
  st_var    : int }

type action =
  | DIVIDE_AND_APPROX of afine * afine * bigint * bigint
  | NOT_EXACT_DIVIDE of afine * bigint
  | FORGET_C of int
  | EXACT_DIVIDE of afine * bigint
  | SUM of int * (bigint * afine) * (bigint * afine)
  | STATE of state_action
  | HYP of afine
  | FORGET   of int * int
  | FORGET_I of int * int
  | CONTRADICTION of afine * afine
  | NEGATE_CONTRADICT of afine * afine * bool
  | MERGE_EQ of int * afine * int
  | CONSTANT_NOT_NUL of int * bigint
  | CONSTANT_NUL of int
  | CONSTANT_NEG of int * bigint
  | SPLIT_INEQ of afine * (int * action list) * (int * action list)
  | WEAKEN of int * bigint

exception UNSOLVABLE

exception NO_CONTRADICTION

let display_eq print_var (l,e) =
  let _ =
    List.fold_left
      (fun not_first f ->
	 print_string
	   (if f.c <? zero then "- " else if not_first then "+ " else "");
	 let c = abs f.c in
	 if c =? one then
	   Printf.printf "%s " (print_var f.v)
	 else
	   Printf.printf "%s %s " (string_of_bigint c) (print_var f.v);
	 true)
      false l
  in
  if e >? zero then
    Printf.printf "+ %s " (string_of_bigint e)
  else if e <? zero then
    Printf.printf "- %s " (string_of_bigint (abs e))

let rec trace_length l =
  let action_length accu = function
    | SPLIT_INEQ (_,(_,l1),(_,l2)) ->
	accu + one + trace_length l1 + trace_length l2
    | _ -> accu + one in
  List.fold_left action_length zero l

let operator_of_eq = function
  | EQUA -> "=" | DISE -> "!=" | INEQ -> ">="

let kind_of = function
  | EQUA -> "equation" | DISE -> "disequation" | INEQ -> "inequation"

let display_system print_var l =
  List.iter
    (fun { kind=b; body=e; constant=c; id=id} ->
       Printf.printf "E%d: " id;
       display_eq print_var (e,c);
       Printf.printf "%s 0\n" (operator_of_eq b))
    l;
  print_string "------------------------\n\n"

let display_inequations print_var l =
  List.iter (fun e -> display_eq print_var e;print_string ">= 0\n") l;
  print_string "------------------------\n\n"

let sbi = string_of_bigint

let rec display_action print_var = function
  | act :: l -> begin match act with
      | DIVIDE_AND_APPROX (e1,e2,k,d) ->
          Printf.printf
            "Inequation E%d is divided by %s and the constant coefficient is \
            rounded by substracting %s.\n" e1.id (sbi k) (sbi d)
      | NOT_EXACT_DIVIDE (e,k) ->
          Printf.printf
            "Constant in equation E%d is not divisible by the pgcd \
            %s of its other coefficients.\n" e.id (sbi k)
      | EXACT_DIVIDE (e,k) ->
          Printf.printf
            "Equation E%d is divided by the pgcd \
            %s of its coefficients.\n" e.id (sbi k)
      | WEAKEN (e,k) ->
          Printf.printf
            "To ensure a solution in the dark shadow \
            the equation E%d is weakened by %s.\n" e (sbi k)
      | SUM (e,(c1,e1),(c2,e2)) ->
          Printf.printf
            "We state %s E%d = %s %s E%d + %s %s E%d.\n"
            (kind_of e1.kind) e (sbi c1) (kind_of e1.kind) e1.id (sbi c2)
            (kind_of e2.kind) e2.id
      | STATE { st_new_eq = e } ->
          Printf.printf "We define a new equation E%d: " e.id;
          display_eq print_var (e.body,e.constant);
          print_string (operator_of_eq e.kind); print_string " 0"
      | HYP e ->
          Printf.printf "We define E%d: " e.id;
          display_eq print_var (e.body,e.constant);
          print_string (operator_of_eq e.kind); print_string " 0\n"
      | FORGET_C e ->  Printf.printf "E%d is trivially satisfiable.\n" e
      | FORGET (e1,e2) -> Printf.printf "E%d subsumes E%d.\n" e1 e2
      | FORGET_I (e1,e2) -> Printf.printf "E%d subsumes E%d.\n" e1 e2
      | MERGE_EQ (e,e1,e2) ->
         Printf.printf "E%d and E%d can be merged into E%d.\n" e1.id e2 e
      | CONTRADICTION (e1,e2) ->
          Printf.printf
            "Equations E%d and E%d imply a contradiction on their \
            constant factors.\n" e1.id e2.id
      | NEGATE_CONTRADICT(e1,e2,b) ->
          Printf.printf
            "Equations E%d and E%d state that their body is at the same time \
            equal and different\n" e1.id e2.id
      | CONSTANT_NOT_NUL (e,k) ->
          Printf.printf "Equation E%d states %s = 0.\n" e (sbi k)
      | CONSTANT_NEG(e,k) ->
          Printf.printf "Equation E%d states %s >= 0.\n" e (sbi k)
      | CONSTANT_NUL e ->
          Printf.printf "Inequation E%d states 0 != 0.\n" e
      | SPLIT_INEQ (e,(e1,l1),(e2,l2)) ->
          Printf.printf "Equation E%d is split in E%d and E%d\n\n" e.id e1 e2;
          display_action print_var l1;
          print_newline ();
          display_action print_var l2;
          print_newline ()
    end; display_action print_var l
  | [] ->
      flush stdout

let default_print_var v = Printf.sprintf "X%d" v (* For debugging *)

(*""*)
let add_event, history, clear_history =
  let accu = ref [] in
  (fun (v:action) -> if !debug then display_action default_print_var [v]; push v accu),
  (fun () -> !accu),
  (fun () -> accu := [])

let nf_linear = List.sort (fun x y -> Pervasives.(-) y.v x.v)

let nf ((b : bool),(e,(x : int))) = (b,(nf_linear e,x))

let map_eq_linear f =
  let rec loop = function
    | x :: l -> let c = f x.c in if c=?zero then loop l else {v=x.v; c=c} :: loop l
    | [] -> []
  in
  loop

let map_eq_afine f e =
  { id = e.id; kind = e.kind; body = map_eq_linear f e.body;
    constant = f e.constant }

let negate_eq = map_eq_afine (fun x -> neg x)

let rec sum p0 p1 = match (p0,p1) with
  | ([], l) -> l | (l, []) -> l
  | (((x1::l1) as l1'), ((x2::l2) as l2')) ->
      if x1.v = x2.v then
	let c = x1.c + x2.c in
	if c =? zero then sum l1 l2 else {v=x1.v;c=c} :: sum l1 l2
      else if x1.v > x2.v then
	x1 :: sum l1 l2'
      else
	x2 :: sum l1' l2

let sum_afine new_eq_id eq1 eq2 =
  { kind = eq1.kind; id = new_eq_id ();
    body = sum eq1.body eq2.body; constant = eq1.constant + eq2.constant }

exception FACTOR1

let rec chop_factor_1 = function
  | x :: l ->
    if abs x.c =? one then x,l else let (c',l') = chop_factor_1 l in (c',x::l')
  | [] -> raise FACTOR1

exception CHOPVAR

let rec chop_var v = function
  | f :: l -> if f.v = v then f,l else let (f',l') = chop_var v l in (f',f::l')
  | [] -> raise CHOPVAR

let normalize ({id=id; kind=eq_flag; body=e; constant =x} as eq) =
  if e = [] then begin
    match eq_flag with
      | EQUA ->
          if x =? zero then [] else begin
            add_event (CONSTANT_NOT_NUL(id,x)); raise UNSOLVABLE
         end
      | DISE ->
          if x <> zero then [] else begin
            add_event (CONSTANT_NUL id); raise UNSOLVABLE
          end
      | INEQ ->
          if x >=? zero then [] else begin
            add_event (CONSTANT_NEG(id,x)); raise UNSOLVABLE
          end
  end else
    let gcd = pgcd_l (List.map (fun f -> abs f.c) e) in
    if eq_flag=EQUA && x mod gcd <> zero then begin
      add_event (NOT_EXACT_DIVIDE (eq,gcd)); raise UNSOLVABLE
    end else if eq_flag=DISE && x mod gcd <> zero then begin
      add_event (FORGET_C eq.id); []
    end else if gcd <> one then begin
      let c = floor_div x gcd in
      let d = x - c * gcd in
      let new_eq = {id=id; kind=eq_flag; constant=c;
                    body=map_eq_linear (fun c -> c / gcd) e} in
      add_event (if eq_flag=EQUA || eq_flag = DISE then EXACT_DIVIDE(eq,gcd)
                 else DIVIDE_AND_APPROX(eq,new_eq,gcd,d));
      [new_eq]
    end else [eq]

let eliminate_with_in new_eq_id {v=v;c=c_unite} eq2
                        ({body=e1; constant=c1} as eq1) =
  try
    let (f,_) = chop_var v e1 in
    let coeff = if c_unite=?one then neg f.c else if c_unite=? negone then f.c
                else failwith "eliminate_with_in" in
    let res = sum_afine new_eq_id eq1 (map_eq_afine (fun c -> c * coeff) eq2) in
    add_event (SUM (res.id,(one,eq1),(coeff,eq2))); res
  with CHOPVAR -> eq1

let omega_mod a b = a - b * floor_div (two * a + b) (two * b)
let banerjee_step (new_eq_id,new_var_id,print_var) original l1 l2 =
  let e = original.body in
  let sigma = new_var_id () in
  if e == [] then begin
    display_system print_var [original] ; failwith "TL"
    end;
  let smallest,var =
    List.fold_left (fun (v,p) c ->  if v >? (abs c.c) then abs c.c,c.v else (v,p))
      (abs (List.hd e).c, (List.hd e).v) (List.tl e)
  in
  let m = smallest + one in
  let new_eq =
    { constant = omega_mod original.constant m;
      body = {c= neg m;v=sigma} ::
             map_eq_linear (fun a -> omega_mod a m) original.body;
      id = new_eq_id (); kind = EQUA } in
  let definition =
    { constant = neg (floor_div (two * original.constant + m) (two * m));
      body = map_eq_linear (fun a -> neg (floor_div (two * a + m) (two * m)))
                             original.body;
      id = new_eq_id (); kind = EQUA } in
  add_event (STATE {st_new_eq = new_eq; st_def = definition;
		    st_orig = original; st_coef = m; st_var = sigma});
  let new_eq = List.hd (normalize new_eq) in
  let eliminated_var, def = chop_var var new_eq.body in
  let other_equations =
    Util.List.map_append
      (fun e ->
        normalize (eliminate_with_in new_eq_id eliminated_var new_eq e)) l1 in
  let inequations =
    Util.List.map_append
      (fun e ->
        normalize (eliminate_with_in new_eq_id eliminated_var new_eq e)) l2 in
  let original' = eliminate_with_in new_eq_id eliminated_var new_eq original in
  let mod_original = map_eq_afine (fun c -> c / m) original' in
  add_event (EXACT_DIVIDE (original',m));
  List.hd (normalize mod_original),other_equations,inequations

let rec eliminate_one_equation ((new_eq_id,new_var_id,print_var) as new_ids) (e,other,ineqs) =
  if !debug then display_system print_var (e::other);
  try
    let v,def = chop_factor_1 e.body in
    (Util.List.map_append
      (fun e' -> normalize (eliminate_with_in new_eq_id v e e')) other,
     Util.List.map_append
       (fun e' -> normalize (eliminate_with_in new_eq_id v e e')) ineqs)
  with FACTOR1 ->
    eliminate_one_equation new_ids (banerjee_step new_ids e other ineqs)

let rec banerjee ((_,_,print_var) as new_ids) (sys_eq,sys_ineq) =
  let rec fst_eq_1 = function
     (eq::l) ->
        if List.exists (fun x -> abs x.c =? one) eq.body then eq,l
        else let (eq',l') = fst_eq_1 l in (eq',eq::l')
   | [] -> raise Not_found in
  match sys_eq with
     [] -> if !debug then display_system print_var sys_ineq; sys_ineq
   | (e1::rest) ->
       let eq,other = try fst_eq_1 sys_eq with Not_found -> (e1,rest) in
       if eq.body = [] then
         if eq.constant =? zero then begin
           add_event (FORGET_C eq.id); banerjee new_ids (other,sys_ineq)
         end else begin
           add_event (CONSTANT_NOT_NUL(eq.id,eq.constant)); raise UNSOLVABLE
         end
       else
	 banerjee new_ids
	   (eliminate_one_equation new_ids (eq,other,sys_ineq))

type kind = INVERTED | NORMAL

let redundancy_elimination new_eq_id system =
  let normal = function
     ({body=f::_} as e) when f.c <? zero ->  negate_eq e, INVERTED
   | e -> e,NORMAL in
  let table = Hashtbl.create 7 in
  List.iter
    (fun e ->
       let ({body=ne} as nx) ,kind = normal e in
       if ne = [] then
         if nx.constant <? zero then begin
           add_event (CONSTANT_NEG(nx.id,nx.constant)); raise UNSOLVABLE
         end else add_event (FORGET_C nx.id)
       else
       try
         let (optnormal,optinvert) = Hashtbl.find table ne in
         let final =
           if kind = NORMAL then begin
             match optnormal with
                Some v ->
                  let kept =
                    if v.constant <? nx.constant
                    then begin add_event (FORGET (v.id,nx.id));v end
                    else begin add_event (FORGET (nx.id,v.id));nx end in
                  (Some(kept),optinvert)
              | None -> Some nx,optinvert
           end else begin
             match optinvert with
                Some v ->
                  let _kept =
                    if v.constant >? nx.constant
                    then begin add_event (FORGET_I (v.id,nx.id));v end
                    else begin add_event (FORGET_I (nx.id,v.id));nx end in
                  (optnormal,Some(if v.constant >? nx.constant then v else nx))
              | None -> optnormal,Some nx
           end in
         begin match final with
            (Some high, Some low) ->
              if high.constant <? low.constant then begin
                add_event(CONTRADICTION (high,negate_eq low));
                raise UNSOLVABLE
              end
          | _ -> () end;
         Hashtbl.remove table ne;
         Hashtbl.add table ne final
       with Not_found ->
         Hashtbl.add table ne
           (if kind = NORMAL then (Some nx,None) else (None,Some nx)))
    system;
  let accu_eq = ref [] in
  let accu_ineq = ref [] in
  Hashtbl.iter
    (fun p0 p1 -> match (p0,p1) with
       | (e, (Some x, Some y)) when x.constant =? y.constant ->
           let id=new_eq_id () in
           add_event (MERGE_EQ(id,x,y.id));
           push {id=id; kind=EQUA; body=x.body; constant=x.constant} accu_eq
       | (e, (optnorm,optinvert)) ->
           begin match optnorm with
               Some x -> push x accu_ineq | _ -> () end;
           begin match optinvert with
               Some x -> push (negate_eq x) accu_ineq | _ -> () end)
    table;
  !accu_eq,!accu_ineq

exception SOLVED_SYSTEM

let select_variable system =
  let table = Hashtbl.create 7 in
  let push v c=
    try let r = Hashtbl.find table v in r := max !r (abs c)
    with Not_found -> Hashtbl.add table v (ref (abs c)) in
  List.iter (fun {body=l} -> List.iter (fun f -> push f.v f.c) l) system;
  let vmin,cmin = ref (-1), ref zero in
  let var_cpt = ref 0 in
  Hashtbl.iter
    (fun v ({contents =  c}) ->
       incr var_cpt;
       if c <? !cmin || !vmin = (-1) then begin vmin := v; cmin := c end)
    table;
  if !var_cpt < 1 then raise SOLVED_SYSTEM;
  !vmin

let classify v system =
  List.fold_left
    (fun (not_occ,below,over) eq ->
       try let f,eq' = chop_var v eq.body in
       if f.c >=? zero then (not_occ,((f.c,eq) :: below),over)
       else (not_occ,below,((neg f.c,eq) :: over))
       with CHOPVAR -> (eq::not_occ,below,over))
    ([],[],[]) system

let product new_eq_id dark_shadow low high =
  List.fold_left
    (fun accu (a,eq1) ->
       List.fold_left
         (fun accu (b,eq2) ->
            let eq =
              sum_afine new_eq_id (map_eq_afine (fun c -> c * b) eq1)
                (map_eq_afine (fun c -> c * a) eq2) in
            add_event(SUM(eq.id,(b,eq1),(a,eq2)));
            match normalize eq with
	      | [eq] ->
                  let final_eq =
                    if dark_shadow then
                      let delta = (a - one) * (b - one) in
                      add_event(WEAKEN(eq.id,delta));
                      {id = eq.id; kind=INEQ; body = eq.body;
                       constant = eq.constant - delta}
                    else eq
                  in final_eq :: accu
              | (e::_) -> failwith "Product dardk"
              | [] -> accu)
         accu high)
    [] low

let fourier_motzkin (new_eq_id,_,print_var) dark_shadow system =
  let v = select_variable system in
  let (ineq_out, ineq_low,ineq_high) = classify v system in
  let expanded = ineq_out @ product new_eq_id dark_shadow ineq_low ineq_high in
  if !debug then display_system print_var expanded; expanded

let simplify ((new_eq_id,new_var_id,print_var) as new_ids) dark_shadow system =
  if List.exists (fun e -> e.kind = DISE) system then
    failwith "disequation in simplify";
  clear_history ();
  List.iter (fun e -> add_event (HYP e)) system;
  let system = Util.List.map_append normalize system in
  let eqs,ineqs = List.partition (fun e -> e.kind=EQUA) system in
  let simp_eq,simp_ineq = redundancy_elimination new_eq_id ineqs in
  let system = (eqs @ simp_eq,simp_ineq) in
  let rec loop1a system =
    let sys_ineq = banerjee new_ids system in
    loop1b sys_ineq
  and loop1b sys_ineq =
    let simp_eq,simp_ineq = redundancy_elimination new_eq_id sys_ineq in
    if simp_eq = [] then simp_ineq else loop1a (simp_eq,simp_ineq)
  in
  let rec loop2 system =
    try
      let expanded = fourier_motzkin new_ids dark_shadow system in
      loop2 (loop1b expanded)
    with SOLVED_SYSTEM ->
      if !debug then display_system print_var system; system
  in
  loop2 (loop1a system)

let rec depend relie_on accu = function
  | act :: l ->
      begin match act with
	| DIVIDE_AND_APPROX (e,_,_,_) ->
            if Int.List.mem e.id relie_on then depend relie_on (act::accu) l
            else depend relie_on accu l
	| EXACT_DIVIDE (e,_) ->
            if Int.List.mem e.id relie_on then depend relie_on (act::accu) l
            else depend relie_on accu l
	| WEAKEN (e,_) ->
            if Int.List.mem e relie_on then depend relie_on (act::accu) l
            else depend relie_on accu l
	| SUM (e,(_,e1),(_,e2)) ->
            if Int.List.mem e relie_on then
	      depend (e1.id::e2.id::relie_on) (act::accu) l
            else
	      depend relie_on accu l
	| STATE {st_new_eq=e;st_orig=o} ->
            if Int.List.mem e.id relie_on then depend (o.id::relie_on) (act::accu) l
            else depend relie_on accu l
	| HYP e ->
            if Int.List.mem e.id relie_on then depend relie_on (act::accu) l
            else depend relie_on accu l
	| FORGET_C _ -> depend relie_on accu l
	| FORGET _ -> depend relie_on accu l
	| FORGET_I _ -> depend relie_on accu l
	| MERGE_EQ (e,e1,e2) ->
            if Int.List.mem e relie_on then
	      depend (e1.id::e2::relie_on) (act::accu) l
            else
	      depend relie_on accu l
	| NOT_EXACT_DIVIDE (e,_) -> depend (e.id::relie_on) (act::accu) l
	| CONTRADICTION (e1,e2) ->
	    depend (e1.id::e2.id::relie_on) (act::accu) l
	| CONSTANT_NOT_NUL (e,_) -> depend (e::relie_on) (act::accu) l
	| CONSTANT_NEG (e,_) -> depend (e::relie_on) (act::accu) l
	| CONSTANT_NUL e -> depend (e::relie_on) (act::accu) l
	| NEGATE_CONTRADICT (e1,e2,_) ->
            depend (e1.id::e2.id::relie_on) (act::accu) l
	| SPLIT_INEQ _ -> failwith "depend"
      end
  | [] -> relie_on, accu

let negation (eqs,ineqs) =
  let diseq,_ = List.partition (fun e -> e.kind = DISE) ineqs in
  let normal = function
    | ({body=f::_} as e) when f.c <? zero ->  negate_eq e, INVERTED
    | e -> e,NORMAL in
  let table = Hashtbl.create 7 in
  List.iter (fun e ->
               let {body=ne;constant=c} ,kind = normal e in
               Hashtbl.add table (ne,c) (kind,e)) diseq;
  List.iter (fun e ->
               assert (e.kind = EQUA);
               let {body=ne;constant=c},kind = normal e in
               try
		 let (kind',e') = Hashtbl.find table (ne,c) in
		 add_event (NEGATE_CONTRADICT (e,e',kind=kind'));
		 raise UNSOLVABLE
               with Not_found -> ()) eqs

exception FULL_SOLUTION of action list * int list

let simplify_strong ((new_eq_id,new_var_id,print_var) as new_ids) system =
  clear_history ();
  List.iter (fun e -> add_event (HYP e)) system;
  (* Initial simplification phase *)
  let rec loop1a system =
    negation system;
    let sys_ineq = banerjee new_ids system in
    loop1b sys_ineq
  and loop1b sys_ineq =
    let dise,ine = List.partition (fun e -> e.kind = DISE) sys_ineq in
    let simp_eq,simp_ineq = redundancy_elimination new_eq_id ine in
    if simp_eq = [] then dise @ simp_ineq
    else loop1a (simp_eq,dise @ simp_ineq)
  in
  let rec loop2 system =
    try
      let expanded = fourier_motzkin new_ids false system in
      loop2 (loop1b expanded)
    with SOLVED_SYSTEM -> if !debug then display_system print_var system; system
  in
  let rec explode_diseq = function
    | (de::diseq,ineqs,expl_map) ->
	let id1 = new_eq_id ()
	and id2 = new_eq_id () in
	let e1 =
          {id = id1; kind=INEQ; body = de.body; constant = de.constant -one} in
	let e2 =
	  {id = id2; kind=INEQ; body = map_eq_linear neg de.body;
           constant = neg de.constant - one} in
	let new_sys =
          List.map (fun (what,sys) -> ((de.id,id1,true)::what, e1::sys))
	    ineqs @
          List.map (fun (what,sys) -> ((de.id,id2,false)::what,e2::sys))
	    ineqs
	in
	explode_diseq (diseq,new_sys,(de.id,(de,id1,id2))::expl_map)
    | ([],ineqs,expl_map) -> ineqs,expl_map
  in
  try
    let system = Util.List.map_append normalize system in
    let eqs,ineqs = List.partition (fun e -> e.kind=EQUA) system in
    let dise,ine = List.partition (fun e -> e.kind = DISE) ineqs in
    let simp_eq,simp_ineq = redundancy_elimination new_eq_id ine in
    let system = (eqs @ simp_eq,simp_ineq @ dise) in
    let system' = loop1a system in
    let diseq,ineq = List.partition (fun e -> e.kind = DISE) system' in
    let first_segment = history () in
    let sys_exploded,explode_map = explode_diseq (diseq,[[],ineq],[]) in
    let all_solutions =
      List.map
        (fun (decomp,sys) ->
           clear_history ();
           try let _ = loop2 sys in raise NO_CONTRADICTION
           with UNSOLVABLE ->
             let relie_on,path = depend [] [] (history ()) in
             let dc,_ = List.partition (fun (_,id,_) -> Int.List.mem id relie_on) decomp in
             let red = List.map (fun (x,_,_) -> x) dc in
             (red,relie_on,decomp,path))
        sys_exploded
    in
    let max_count sys =
      let tbl = Hashtbl.create 7 in
      let augment x =
        try incr (Hashtbl.find tbl x)
	with Not_found -> Hashtbl.add tbl x (ref 1) in
      let eq = ref (-1) and c = ref 0 in
      List.iter (function
		   | ([],r_on,_,path) -> raise (FULL_SOLUTION (path,r_on))
                   | (l,_,_,_) -> List.iter augment l) sys;
      Hashtbl.iter (fun x v -> if !v > !c then begin eq := x; c := !v end) tbl;
      !eq
    in
    let rec solve systems =
      try
	let id = max_count systems in
        let rec sign = function
	  | ((id',_,b)::l) -> if id=id' then b else sign l
          | [] -> failwith "solve" in
        let s1,s2 =
          List.partition (fun (_,_,decomp,_) -> sign decomp) systems in
        let remove_int (dep,ro,dc,pa) =
          (Util.List.except Int.equal id dep,ro,dc,pa)
        in
        let s1' = List.map remove_int s1 in
        let s2' = List.map remove_int s2 in
        let (r1,relie1) = solve s1'
	and (r2,relie2) = solve s2' in
	let (eq,id1,id2) = Int.List.assoc id explode_map in
	[SPLIT_INEQ(eq,(id1,r1),(id2, r2))],
        eq.id :: Util.List.union Int.equal relie1 relie2
      with FULL_SOLUTION (x0,x1) -> (x0,x1)
    in
    let act,relie_on = solve all_solutions in
    snd(depend relie_on act first_segment)
  with UNSOLVABLE -> snd (depend [] [] (history ()))

end