summaryrefslogtreecommitdiff
path: root/plugins/micromega/polynomial.ml
blob: 14d312a5ccff466f79ac2521bf1fbfad02df5d5f (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
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
(************************************************************************)
(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
(*   \VV/  **************************************************************)
(*    //   *      This file is distributed under the terms of the       *)
(*         *       GNU Lesser General Public License Version 2.1        *)
(************************************************************************)
(*                                                                      *)
(* Micromega: A reflexive tactic using the Positivstellensatz           *)
(*                                                                      *)
(*  Frédéric Besson (Irisa/Inria) 2006-20011                            *)
(*                                                                      *)
(************************************************************************)

open Num
module Utils = Mutils
open Utils

type var = int


let (<+>) = add_num
let (<->) = minus_num
let (<*>) = mult_num


module Monomial :
sig
 type t
 val const : t
 val is_const : t -> bool
 val var : var -> t
 val is_var : t -> bool
 val find : var -> t -> int
 val mult : var -> t -> t
 val prod : t -> t -> t
 val exp  : t -> int -> t
 val div  : t -> t -> t *  int
 val compare : t -> t -> int
 val pp : out_channel -> t -> unit
 val fold : (var -> int -> 'a -> 'a) -> t -> 'a -> 'a
 val sqrt : t -> t option 
end
 =
struct
 (* A monomial is represented by a multiset of variables  *)
 module Map = Map.Make(struct type t = var let compare = Pervasives.compare end)
 open Map
 
 type t = int Map.t

 let pp o m = Map.iter 
   (fun k v -> 
      if v = 1 then Printf.fprintf o "x%i." k
      else     Printf.fprintf o "x%i^%i." k v) m


 (* The monomial that corresponds to a constant *)
 let const = Map.empty

 let sum_degree m = Map.fold (fun _ n s -> s + n) m 0

 (* Total ordering of monomials *)
 let compare: t -> t -> int = 
   fun m1 m2 -> 
     let s1 = sum_degree m1 
     and s2 = sum_degree m2 in
       if s1 = s2 then Map.compare Pervasives.compare m1 m2
       else Pervasives.compare s1 s2

 let is_const m = (m = Map.empty)
  
 (* The monomial 'x' *)
 let var x = Map.add x 1 Map.empty

 let is_var m = 
   try
     not (Map.fold (fun _ i fk -> 
		 if fk = true (* first key *)
		 then 
		   if i = 1 then false
		   else raise Not_found
		 else raise Not_found) m true)
   with Not_found -> false

 let sqrt m = 
   if is_const m then None
   else
     try 
       Some (Map.fold (fun v i acc ->
			 let i' = i / 2 in
			   if i mod 2 = 0
			   then add v i' m
			   else raise Not_found) m const)
     with Not_found -> None

 (* Get the degre of a variable in a monomial *)
 let find x m = try find x m with Not_found -> 0
  
 (* Multiply a monomial by a variable *)
 let mult x m = add x (  (find x m) + 1) m
   
 (* Product of monomials *)
 let prod m1 m2 = Map.fold (fun k d m -> add k ((find k m) + d) m) m1 m2


 let exp m n =
   let rec exp acc n = 
     if n = 0 then acc
     else  exp (prod acc m) (n - 1) in

     exp const n


 (*  [div m1 m2 = mr,n] such that mr * (m2)^n = m1 *)
 let div m1 m2 =
   let n = fold (fun x i n -> let i' = find x m1 in 
		     let nx = i' / i in
		       min n nx) m2 max_int in

   let mr = fold (fun x i' m -> 
		    let i = find x m2 in
		    let ir = i' - i * n in
		      if ir = 0 then m
		      else add x ir m) m1 empty in
     (mr,n)


 let fold = fold

end

module Poly :
 (* A polynomial is a map of monomials *)
 (* 
    This is probably a naive implementation 
    (expected to be fast enough - Coq is probably the bottleneck)
    *The new ring contribution is using a sparse Horner representation.
 *)
sig
 type t
 val get : Monomial.t -> t -> num
 val variable : var -> t
 val add : Monomial.t -> num -> t -> t
 val constant : num -> t
 val mult : Monomial.t -> num -> t -> t
 val product : t -> t -> t
 val addition : t -> t -> t
 val uminus : t -> t
 val fold : (Monomial.t -> num -> 'a -> 'a) -> t -> 'a -> 'a
 val pp : out_channel -> t -> unit
 val compare : t -> t -> int
 val is_null : t -> bool
 val is_linear : t -> bool
end =
struct
 (*normalisation bug : 0*x ... *)
 module P = Map.Make(Monomial)
 open P

 type t = num P.t

 let pp o p = P.iter 
   (fun k v -> 
	if Monomial.compare Monomial.const k = 0
	then         Printf.fprintf o "%s " (string_of_num v)
	else Printf.fprintf o "%s*%a " (string_of_num v) Monomial.pp k) p  

 (* Get the coefficient of monomial mn *)
 let get : Monomial.t -> t -> num = 
  fun mn p -> try find mn p with Not_found -> (Int 0)


 (* The polynomial 1.x *)
 let variable : var -> t =
  fun  x ->  add (Monomial.var x) (Int 1) empty
   
 (*The constant polynomial *)
 let constant : num -> t =
  fun c -> add (Monomial.const) c empty

 (* The addition of a monomial *)

 let add : Monomial.t -> num -> t -> t =
  fun mn v p -> 
    if sign_num v = 0 then p
    else
      let vl = (get mn p) <+> v in
	if sign_num vl = 0 then
	  remove mn p
	else add mn vl p


 (** Design choice: empty is not a polynomial 
     I do not remember why .... 
 **)

 (* The product by a monomial *)
 let mult : Monomial.t -> num -> t -> t =
  fun mn v p -> 
    if sign_num v = 0 
    then constant (Int 0)
    else
      fold (fun mn' v' res -> P.add (Monomial.prod mn mn') (v<*>v') res) p empty


 let  addition : t -> t -> t =
  fun p1 p2 -> fold (fun mn v p -> add mn v p) p1 p2
   

 let product : t -> t -> t =
  fun p1 p2 -> 
   fold (fun mn v res -> addition (mult mn v p2) res ) p1 empty


 let uminus : t -> t =
  fun p -> map (fun v -> minus_num v) p

 let fold = P.fold

 let is_null p = fold (fun mn vl b -> b & sign_num vl = 0) p  true

 let compare = compare compare_num

 let is_linear p = P.fold (fun m _ acc -> acc && (Monomial.is_const m || Monomial.is_var m)) p true

(* let is_linear p = 
   let res = is_linear p in
     Printf.printf "is_linear %a = %b\n" pp p res ; res
*)
end


module Vect =
  struct
    (** [t] is the type of vectors.
	A vector [(x1,v1) ; ... ; (xn,vn)] is such that:
	- variables indexes are ordered (x1 <c ... < xn
	- values are all non-zero
    *)
    type var = int
    type t = (var * num) list

(** [equal v1 v2 = true] if the vectors are syntactically equal.
    ([num] is not handled by  [Pervasives.equal] *)

    let rec equal v1 v2 =
      match v1 , v2 with
	| [] , []   -> true
	| [] , _    -> false
	| _::_ , [] -> false
	| (i1,n1)::v1 , (i2,n2)::v2 ->
	    (i1 = i2) && n1 =/ n2 && equal v1 v2

    let hash v =
      let rec hash i = function
	| [] -> i
	| (vr,vl)::l ->  hash (i + (Hashtbl.hash (vr, float_of_num vl))) l in
	Hashtbl.hash (hash 0 v )


    let null = []

    let pp_vect o vect =
      List.iter (fun (v,n) -> Printf.printf "%sx%i + " (string_of_num n) v) vect

    let from_list (l: num list) =
      let rec xfrom_list i l =
	match l with
	  | [] -> []
	  | e::l ->
	      if e <>/ Int 0
	      then (i,e)::(xfrom_list (i+1) l)
	      else xfrom_list (i+1) l in

	xfrom_list 0 l

    let zero_num = Int 0
    let unit_num = Int 1


    let to_list m =
      let rec xto_list i l =
	match l with
	  | [] -> []
	  |	(x,v)::l' ->
		  if i = x then v::(xto_list (i+1) l') else zero_num ::(xto_list (i+1) l) in
	xto_list 0 m


    let cons i v rst = if v =/ Int 0 then rst else (i,v)::rst

    let rec update i f t =
      match t with
	| [] -> cons i (f zero_num) []
	| (k,v)::l ->
	    match Pervasives.compare i k with
	  | 0 -> cons k (f v) l
	  | -1 -> cons i (f zero_num) t
	  |  1 -> (k,v) ::(update i f l)
	  |  _  -> failwith "compare_num"

    let rec set i n t =
      match t with
	| [] -> cons i n []
	| (k,v)::l ->
	    match Pervasives.compare i k with
	      | 0 -> cons k n l
	      | -1 -> cons i n t
	      |  1 -> (k,v) :: (set i n l)
	      |  _  -> failwith "compare_num"

    let gcd m =
      let res = List.fold_left (fun x (i,e) -> Big_int.gcd_big_int  x (Utils.numerator e))  Big_int.zero_big_int m in
	if Big_int.compare_big_int res Big_int.zero_big_int = 0
	then Big_int.unit_big_int else res

    let rec mul z t =
      match z with
	| Int 0 -> []
	| Int 1 -> t
	|  _    -> List.map (fun (i,n) -> (i, mult_num z n)) t


    let rec add v1 v2  =
      match v1 , v2 with
	| (x1,n1)::v1' , (x2,n2)::v2' ->
	    if x1 = x2
	    then
	      let n' = n1  +/ n2 in
		if n' =/ Int 0 then add v1' v2'
		else
		  let res = add v1' v2'  in
		    (x1,n') ::res
	    else if x1 < x2
	    then let res = add v1' v2   in
	      (x1, n1)::res
	    else let res = add v1 v2'  in
	      (x2, n2)::res
	|  [] , [] -> []
	|  [] ,  _ ->  v2
	|  _  , [] ->  v1 




    let compare : t -> t -> int = Utils.Cmp.compare_list (fun x y -> Utils.Cmp.compare_lexical
      [
	(fun () -> Pervasives.compare (fst x) (fst y));
	(fun () -> compare_num (snd x) (snd y))])

    (** [tail v vect] returns
	- [None] if [v] is not a variable of the vector [vect]
	- [Some(vl,rst)]  where [vl] is the value of [v] in vector [vect]
        and [rst] is the remaining of the vector
	We exploit that vectors are ordered lists
    *)
    let rec tail (v:var) (vect:t) =
      match vect with
	| [] -> None
	| (v',vl)::vect' ->
	    match Pervasives.compare v' v with
	      | 0 -> Some (vl,vect) (* Ok, found *)
	      | -1 -> tail v vect' (* Might be in the tail *)
	      | _  -> None (* Hopeless *)

    let get v vect =
      match tail v vect with
	| None -> None
	| Some(vl,_) -> Some vl


    let rec fresh v =
	match v with
	  | [] -> 1
	  | [v,_] -> v + 1
	  | _::v  -> fresh v

  end

type vector = Vect.t

type cstr_compat = {coeffs : vector ; op : op ; cst : num}
and op = |Eq | Ge

let string_of_op = function Eq -> "=" | Ge -> ">="

let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} = 
  Printf.fprintf o "%a %s %s" Vect.pp_vect coeffs (string_of_op op) (string_of_num cst)

let opMult o1 o2 = 
  match o1, o2 with
    | Eq , Eq -> Eq
    | Eq , Ge | Ge , Eq -> Ge
    | Ge , Ge -> Ge

let opAdd o1 o2 =
  match o1 , o2 with
    | Eq , _ | _ , Eq -> Eq
    | Ge , Ge -> Ge




open Big_int

type index = int

type prf_rule = 
  | Hyp of int 
  | Def of int
  | Cst  of big_int
  | Zero
  | Square of (Vect.t * num) 
  | MulC of (Vect.t * num) * prf_rule 
  | Gcd of big_int * prf_rule 
  | MulPrf of prf_rule * prf_rule
  | AddPrf of prf_rule * prf_rule
  | CutPrf of prf_rule
      
type proof = 
  | Done
  | Step of int * prf_rule * proof
  | Enum of int * prf_rule * Vect.t * prf_rule * proof list 
      

let rec output_prf_rule o = function
  | Hyp i -> Printf.fprintf o "Hyp %i" i
  | Def i -> Printf.fprintf o "Def %i" i
  | Cst c -> Printf.fprintf o "Cst %s" (string_of_big_int c)
  | Zero  -> Printf.fprintf o "Zero"
  | Square _ -> Printf.fprintf o "( )^2"
  | MulC(p,pr) -> Printf.fprintf o "P * %a" output_prf_rule pr
  | MulPrf(p1,p2) -> Printf.fprintf o "%a * %a" output_prf_rule p1 output_prf_rule p2
  | AddPrf(p1,p2) -> Printf.fprintf o "%a + %a" output_prf_rule p1 output_prf_rule p2
  | CutPrf(p) -> Printf.fprintf o "[%a]" output_prf_rule p
  | Gcd(c,p)  -> Printf.fprintf o "(%a)/%s" output_prf_rule p (string_of_big_int c)

let rec output_proof o = function
  | Done -> Printf.fprintf o "."
  | Step(i,p,pf) -> Printf.fprintf o "%i:= %a ; %a" i output_prf_rule p output_proof pf
  | Enum(i,p1,v,p2,pl) -> Printf.fprintf o "%i{%a<=%a<=%a}%a" i 
      output_prf_rule p1 Vect.pp_vect v output_prf_rule p2
	(pp_list output_proof) pl

let rec pr_rule_max_id = function
  | Hyp i | Def i -> i
  | Cst _ | Zero | Square _ -> -1
  | MulC(_,p) | CutPrf p | Gcd(_,p) -> pr_rule_max_id p
  | MulPrf(p1,p2)| AddPrf(p1,p2) -> max (pr_rule_max_id p1) (pr_rule_max_id p2)

let rec proof_max_id = function
  | Done -> -1
  | Step(i,pr,prf) -> max i (max (pr_rule_max_id pr) (proof_max_id prf))
  | Enum(i,p1,_,p2,l) -> 
      let m = max (pr_rule_max_id p1) (pr_rule_max_id p2) in
	List.fold_left (fun i prf -> max i (proof_max_id prf)) (max i m) l

let rec pr_rule_def_cut id = function
  | MulC(p,prf) -> 
      let (bds,id',prf') = pr_rule_def_cut id prf  in
	(bds, id', MulC(p,prf'))
  | MulPrf(p1,p2) -> 
      let (bds1,id,p1) = pr_rule_def_cut id p1 in
      let (bds2,id,p2) = pr_rule_def_cut id p2 in
	(bds2@bds1,id,MulPrf(p1,p2))
  | AddPrf(p1,p2) -> 
      let (bds1,id,p1) = pr_rule_def_cut id p1 in
      let (bds2,id,p2) = pr_rule_def_cut id p2 in
	(bds2@bds1,id,AddPrf(p1,p2))
  | CutPrf p -> 
      let (bds,id,p) = pr_rule_def_cut id p in
	((id,p)::bds,id+1,Def id)
  | Gcd(c,p) ->       
      let (bds,id,p) = pr_rule_def_cut id p in
	((id,p)::bds,id+1,Def id)
  | Square _|Cst _|Def _|Hyp _|Zero as x -> ([],id,x)


(* Do not define top-level cuts *)
let pr_rule_def_cut id = function
  | CutPrf p -> 
      let (bds,ids,p') = pr_rule_def_cut id p in
	bds,ids, CutPrf p'
  |    p     -> pr_rule_def_cut id p


let rec implicit_cut p = 
  match p with
    | CutPrf p -> implicit_cut p
    |   _      -> p
       

let rec normalise_proof id prf =
  match prf with
  | Done -> (id,Done)
  | Step(i,Gcd(c,p),Done) ->  normalise_proof id (Step(i,p,Done))
  | Step(i,p,prf) -> 
      let bds,id,p' = pr_rule_def_cut id p in
      let (id,prf) = normalise_proof id prf in
      let prf = List.fold_left (fun  acc (i,p)  -> Step(i, CutPrf p,acc)) 
	(Step(i,p',prf)) 	  bds in

	(id,prf)
  | Enum(i,p1,v,p2,pl) -> 
      (* Why do I have  top-level cuts ? *)
(*      let p1 = implicit_cut p1 in
      let p2 = implicit_cut p2 in
      let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in
	(List.fold_left max  0 ids , 
	   Enum(i,p1,v,p2,prfs))
*)

      let bds1,id,p1' = pr_rule_def_cut id (implicit_cut p1) in
      let bds2,id,p2' = pr_rule_def_cut id (implicit_cut p2) in
      let (ids,prfs) = List.split (List.map (normalise_proof id) pl) in
	(List.fold_left max  0 ids , 
	 List.fold_left (fun  acc (i,p)  -> Step(i, CutPrf p,acc))
	   (Enum(i,p1',v,p2',prfs)) (bds2@bds1))


let normalise_proof id prf = 
  let res = normalise_proof id prf in
  if debug then Printf.printf "normalise_proof %a -> %a" output_proof  prf output_proof (snd res) ; 
    res



let add_proof x y = 
  match x, y with 
   | Zero , p | p , Zero -> p
   | _    -> AddPrf(x,y)


let mul_proof c p = 
  match  sign_big_int c with
   | 0 -> Zero (* This is likely to be a bug *)
   | -1 -> MulC(([],Big_int c),p) (* [p] should represent an equality *)
   | 1 -> 
       if eq_big_int c unit_big_int
       then p
       else MulPrf(Cst c,p)
   | _ -> assert false


let mul_proof_ext (p,c) prf = 
  match p with
   | [] -> mul_proof (numerator c) prf
   |  _ -> MulC((p,c),prf)
  


(*
  let rec scale_prf_rule = function
	    | Hyp i -> (unit_big_int, Hyp i)
  | Def i -> (unit_big_int, Def i)
  | Cst c -> (unit_big_int, Cst i)
  | Zero  -> (unit_big_int, Zero)
  | Square p -> (unit_big_int,Square p)
  | Div(c,pr) -> 
  let (bi,pr') = scale_prf_rule pr in
  (mult_big_int c bi , pr')
  | MulC(p,pr) -> 
  let bi,pr' = scale_prf_rule pr in
  (bi,MulC p,pr')
  | MulPrf(p1,p2) -> 
  let b1,p1 = scale_prf_rule p1 in
  let b2,p2 = scale_prf_rule p2 in
	    
  
  | AddPrf(p1,p2) -> 
	    let b1,p1 = scale_prf_rule p1 in
  let b2,p2 = scale_prf_rule p2 in
  let g =  gcd_big_int
*)
	  




module LinPoly = 
struct
  type t = Vect.t * num

  module MonT = 
  struct
    module MonoMap = Map.Make(Monomial)
    module IntMap  = Map.Make(struct type t = int let compare = Pervasives.compare end)
	  
    (** A hash table might be preferable but requires a hash function. *)
    let (index_of_monomial : int MonoMap.t ref) = ref (MonoMap.empty)
    let (monomial_of_index : Monomial.t IntMap.t ref) = ref (IntMap.empty)
    let fresh = ref 0

    let clear () = 
      index_of_monomial := MonoMap.empty;
      monomial_of_index := IntMap.empty ; 
      fresh := 0


    let register m = 
      try
	MonoMap.find m !index_of_monomial
      with Not_found -> 
	begin
	  let res = !fresh in
	    index_of_monomial := MonoMap.add m res !index_of_monomial ; 
	    monomial_of_index := IntMap.add res m !monomial_of_index ;
	    incr fresh ; res
	end

    let retrieve i = IntMap.find i !monomial_of_index


  end 

  let normalise (v,c) = 
    (List.sort  (fun x y -> Pervasives.compare (fst x) (fst y)) v , c)


  let output_mon o (x,v) = 
    Printf.fprintf o "%s.%a +" (string_of_num v) Monomial.pp (MonT.retrieve x)



  let output_cstr o {coeffs = coeffs ; op = op ; cst = cst} = 
    Printf.fprintf o "%a %s %s" (pp_list output_mon) coeffs (string_of_op op) (string_of_num cst)



  let linpol_of_pol p = 
    let (v,c) = 
      Poly.fold 
	(fun  mon num (vct,cst)  -> 
	   if Monomial.is_const mon then (vct,num) 
	   else
	     let vr = MonT.register mon in
	       ((vr,num)::vct,cst)) p ([], Int 0) in
      normalise (v,c)

  let mult v m (vect,c) =
    if Monomial.is_const m
    then
      (Vect.mul v vect, v <*> c)
    else
      if sign_num v <> 0
      then
	let hd = 
	  if sign_num c <> 0
	  then [MonT.register m,v <*> c]
	  else [] in
	  
	let vect = hd @ (List.map (fun (x,n) -> 
					let x = MonT.retrieve x in
					let x_m = MonT.register (Monomial.prod m x) in
					  (x_m, v <*> n)) vect ) in
	  normalise (vect , Int 0)
      else ([],Int 0)

  let mult v m (vect,c) = 
    let (vect',c') = mult v m (vect,c) in
      if debug then 
      Printf.printf "mult %s %a (%a,%s) -> (%a,%s)\n" (string_of_num v) Monomial.pp m 
	(pp_list output_mon) vect (string_of_num c) 
	(pp_list output_mon) vect' (string_of_num c') ;
      (vect',c')



  let make_lin_pol v mon =
    if Monomial.is_const mon
    then [] , v
    else [MonT.register mon, v],Int 0

      




  let xpivot_eq (c,prf) x v (c',prf') = 
    if debug then Printf.printf "xpivot_eq {%a} %a %s {%a}\n" 
      output_cstr c 
      Monomial.pp (MonT.retrieve x)
      (string_of_num v) output_cstr c' ; 


    let {coeffs = coeffs ; op = op ; cst = cst} = c' in
    let m = MonT.retrieve x in

    let apply_pivot (vqn,q,n) (c',prf') = 
      (* Morally, we have (Vect.get (q*x^n) c'.coeffs) = vmn with n >=0 *)

      let cc' = abs_num v in
      let cc_num  = Int (- (sign_num v)) <*> vqn in
      let cc_mon  = Monomial.prod q (Monomial.exp m (n-1)) in

      let (c_coeff,c_cst) = mult cc_num cc_mon (c.coeffs, minus_num c.cst) in
	
      let c' = {coeffs = Vect.add (Vect.mul cc' c'.coeffs) c_coeff ; op = op ; cst = (minus_num c_cst) <+> (cc' <*> c'.cst)} in
      let prf' = add_proof 
	(mul_proof_ext (make_lin_pol cc_num cc_mon) prf)
	(mul_proof (numerator cc') prf') in

	if debug then Printf.printf "apply_pivot -> {%a}\n" output_cstr c' ; 
	(c',prf') in


    let cmp (q,n) (q',n') = 
      if n < n' then -1
      else if n = n' then  Monomial.compare q q'
      else 1 in

      
    let find_pivot (c',prf') = 
      let (v,q,n) = List.fold_left 
	(fun (v,q,n) (x,v') -> 
	   let x = MonT.retrieve x in 
	   let (q',n') = Monomial.div x m in
	     if cmp (q,n) (q',n') = -1 then (v',q',n') else (v,q,n)) (Int 0, Monomial.const,0) c'.coeffs in
	if n > 0 then Some (v,q,n) else None in

    let rec pivot (q,n) (c',prf') = 
      match find_pivot (c',prf') with
       | None -> (c',prf') 
       | Some(v,q',n') -> 
	   if cmp (q',n')  (q,n) = -1
	   then pivot (q',n') (apply_pivot (v,q',n') (c',prf'))
	   else (c',prf') in

      pivot (Monomial.const,max_int) (c',prf')


  let pivot_eq x (c,prf)   = 
    match Vect.get x c.coeffs with
     | None -> (fun x -> None)
     | Some v -> fun cp' -> Some (xpivot_eq (c,prf) x v cp')


end