Quellcode-Bibliothek

^© Kompilation durch diese Firma

[Weder Korrektheit noch Funktionsfähigkeit der Software werden zugesichert.]

Datei: Fourier_util.v Sprache: SML

Original von: Coq^©

(************************************************************************)
(*         *   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)         *)
(************************************************************************)
(*                                                                      *)
(* Micromega: A reflexive tactic using the Positivstellensatz           *)
(*                                                                      *)
(*  Frédéric Besson (Irisa/Inria) 2006-2008                             *)
(*                                                                      *)
(************************************************************************)

(* We take as input a list of polynomials [p1...pn] and return an unfeasibility
   certificate polynomial. *)

let debug = false

open Big_int
open Num
open Polynomial

module Mc = Micromega
module Ml2C = Mutils.CamlToCoq
module C2Ml = Mutils.CoqToCaml

let use_simplex = ref true

type ('prf,'model) res =
   | Prf of 'prf
   | Model of 'model
   | Unknown

type zres = (Mc.zArithProof , (int * Mc.z  list)) res

type qres = (Mc.q Mc.psatz , (int * Mc.q  list)) res

open Mutils
type 'a number_spec = {
    bigint_to_number : big_int -> 'a;
    number_to_num : 'a -> num;
    zero : 'a;
    unit : 'a;
    mult : 'a -> 'a -> 'a;
    eqb  : 'a -> 'a -> bool
  }

let z_spec = {
    bigint_to_number = Ml2C.bigint ;
    number_to_num = (fun x -> Big_int (C2Ml.z_big_int x));
    zero = Mc.Z0;
    unit = Mc.Zpos Mc.XH;
    mult = Mc.Z.mul;
    eqb  = Mc.zeq_bool
  }

let q_spec = {
    bigint_to_number = (fun x -> {Mc.qnum = Ml2C.bigint x; Mc.qden = Mc.XH});
    number_to_num = C2Ml.q_to_num;
    zero = {Mc.qnum = Mc.Z0;Mc.qden = Mc.XH};
    unit = {Mc.qnum =  (Mc.Zpos Mc.XH) ; Mc.qden = Mc.XH};
    mult = Mc.qmult;
    eqb  = Mc.qeq_bool
  }

let dev_form n_spec  p =
  let rec dev_form p =
    match p with
    | Mc.PEc z ->  Poly.constant (n_spec.number_to_num z)
    | Mc.PEX v ->  Poly.variable (C2Ml.positive v)
    | Mc.PEmul(p1,p2) ->
       let p1 = dev_form p1 in
       let p2 = dev_form p2 in
       Poly.product p1 p2
    | Mc.PEadd(p1,p2) -> Poly.addition (dev_form p1) (dev_form p2)
    | Mc.PEopp p ->  Poly.uminus (dev_form p)
    | Mc.PEsub(p1,p2) ->  Poly.addition (dev_form p1) (Poly.uminus (dev_form p2))
    | Mc.PEpow(p,n)   ->
       let p = dev_form p in
       let n = C2Ml.n n in
       let rec pow n =
         if Int.equal n 0
         then Poly.constant (n_spec.number_to_num n_spec.unit)
         else Poly.product p (pow (n-1)) in
       pow n in
  dev_form p

let rec fixpoint f x =
  let y' = f x in
  if Pervasives.(=) y' x then y'
  else fixpoint f y'

let  rec_simpl_cone n_spec e =
  let simpl_cone =
    Mc.simpl_cone n_spec.zero n_spec.unit n_spec.mult n_spec.eqb in

  let rec rec_simpl_cone  = function
    | Mc.PsatzMulE(t1, t2) ->
       simpl_cone  (Mc.PsatzMulE (rec_simpl_cone t1, rec_simpl_cone t2))
    | Mc.PsatzAdd(t1,t2)  ->
       simpl_cone (Mc.PsatzAdd (rec_simpl_cone t1, rec_simpl_cone t2))
    |  x           -> simpl_cone x in
  rec_simpl_cone e


let simplify_cone n_spec c = fixpoint (rec_simpl_cone n_spec) c

(* The binding with Fourier might be a bit obsolete
   -- how does it handle equalities ? *)

(* Certificates are elements of the cone such that P = 0  *)

(* To begin with, we search for certificates of the form:
   a1.p1 + ... an.pn + b1.q1 +... + bn.qn + c = 0
   where pi >= 0 qi > 0
   ai >= 0
   bi >= 0
   Sum bi + c >= 1
   This is a linear problem: each monomial is considered as a variable.
   Hence, we can use fourier.

   The variable c is at index 1
*)

(* fold_left followed by a rev ! *)

let constrain_variable v l =
  let coeffs = List.fold_left (fun acc p -> (Vect.get v p.coeffs)::acc) [] l in
  { coeffs = Vect.from_list ((Big_int zero_big_int):: (Big_int zero_big_int):: (List.rev coeffs)) ;
    op = Eq ;
    cst = Big_int zero_big_int  }

let constrain_constant l =
  let coeffs = List.fold_left (fun acc p -> minus_num p.cst ::acc) [] l in
  { coeffs = Vect.from_list ((Big_int zero_big_int):: (Big_int unit_big_int):: (List.rev coeffs)) ;
    op = Eq ;
    cst = Big_int zero_big_int  }

let positivity l =
  let rec xpositivity i l =
    match l with
    | [] -> []
    | c::l -> match c.op with
              | Eq -> xpositivity (i+1) l
              |  _ ->
                  {coeffs = Vect.update (i+1) (fun _ -> Int 1) Vect.null ;
                   op = Ge ;
                   cst = Int 0 }  :: (xpositivity (i+1) l)
  in
  xpositivity 1 l

let cstr_of_poly (p,o) =
  let (c,l) = Vect.decomp_cst p in
  {coeffs = l; op = o ; cst = minus_num c}

let variables_of_cstr c = Vect.variables c.coeffs

(* If the certificate includes at least one strict inequality,
   the obtained polynomial can also be 0 *)

let build_dual_linear_system l =

  let variables =
    List.fold_left (fun acc p -> ISet.union acc (variables_of_cstr p)) ISet.empty l in
  (* For each monomial, compute a constraint *)
  let s0 =
    ISet.fold (fun mn  res -> (constrain_variable mn l)::res) variables [] in
  let c  = constrain_constant l in

  (* I need at least something strictly positive *)
  let strict = {
      coeffs = Vect.from_list ((Big_int zero_big_int) :: (Big_int unit_big_int)::
                                 (List.map (fun c ->  if is_strict c then Big_int unit_big_int else Big_int zero_big_int) l));
      op = Ge ; cst = Big_int unit_big_int } in
  (* Add the positivity constraint *)
  {coeffs = Vect.from_list ([Big_int zero_big_int ;Big_int unit_big_int]) ;
   op = Ge ;
   cst = Big_int zero_big_int}::(strict::(positivity l)@c::s0)
open Util

(** [direct_linear_prover l] does not handle strict inegalities *)
let fourier_linear_prover l =
  match Mfourier.Fourier.find_point l with
  | Inr prf ->
     if debug then Printf.printf "AProof : %a\n" Mfourier.pp_proof prf ;
     let cert = (*List.map (fun (x,n) -> x+1,n)*) (fst (List.hd (Mfourier.Proof.mk_proof l prf))) in
     if debug then Printf.printf "CProof : %a" Vect.pp cert ;
  (*Some (rats_to_ints (Vect.to_list cert))*)
     Some (Vect.normalise cert)
  | Inl _   -> None

let direct_linear_prover l =
  if !use_simplex
  then Simplex.find_unsat_certificate l
  else fourier_linear_prover l

let find_point l =
   if !use_simplex
   then Simplex.find_point l
   else match Mfourier.Fourier.find_point l with
        | Inr _ -> None
        | Inl cert -> Some cert

let optimise v l =
  if !use_simplex
  then Simplex.optimise v l
  else Mfourier.Fourier.optimise v l

let dual_raw_certificate l =
  if debug
  then begin
      Printf.printf "dual_raw_certificate\n";
      List.iter (fun c -> Printf.fprintf stdout "%a\n" output_cstr c) l
    end;

  let sys = build_dual_linear_system l in

  if debug then begin
      Printf.printf "dual_system\n";
      List.iter (fun c -> Printf.fprintf stdout "%a\n" output_cstr c) sys
    end;

  try
    match find_point sys with
    | None -> None
    | Some cert ->
       match Vect.choose cert with
       | None -> failwith "dual_raw_certificate: empty_certificate"
       | Some _ ->
          (*Some (rats_to_ints (Vect.to_list (Vect.decr_var 2 (Vect.set 1 (Int 0) cert))))*)
          Some (Vect.normalise (Vect.decr_var 2 (Vect.set 1 (Int 0) cert)))
               (* should not use rats_to_ints *)
  with x when CErrors.noncritical x ->
        if debug
        then (Printf.printf "dual raw certificate %s" (Printexc.to_string x);
              flush stdout) ;
        None

let simple_linear_prover l =
  try
    direct_linear_prover l
  with Strict ->
    (* Fourier elimination should handle > *)
    dual_raw_certificate l

let env_of_list l =
  snd (List.fold_left (fun (i,m) p -> (i+1,IMap.add i p m)) (0,IMap.empty) l)

let linear_prover_cstr sys =
  let (sysi,prfi) = List.split sys in

  match simple_linear_prover sysi with
  | None -> None
  | Some cert -> Some (ProofFormat.proof_of_farkas (env_of_list prfi) cert)

let linear_prover_cstr  =
  if debug
  then
    fun sys ->
    Printf.printf "; flush stdout ;
    let res = linear_prover_cstr sys in
    Printf.printf ">"; flush stdout ;
    res
  else linear_prover_cstr

let compute_max_nb_cstr l d =
  let len = List.length l in
  max len (max d (len * d))

let develop_constraint z_spec (e,k) =
  (dev_form z_spec e,
   match k with
   | Mc.NonStrict -> Ge
   | Mc.Equal     -> Eq
   | Mc.Strict    -> Gt
   | _     -> assert false
  )

(** A single constraint can be unsat for the following reasons:
    - 0 >= c for c a negative constant
    - 0 =  c for c a non-zero constant
    - e = c  when the coeffs of e are all integers and c is rational
*)

type checksat =
  | Tauto (* Tautology *)
  | Unsat of ProofFormat.prf_rule (* Unsatisfiable *)
  | Cut of cstr * ProofFormat.prf_rule (* Cutting plane *)
  | Normalise of cstr * ProofFormat.prf_rule (* Coefficients may be normalised i.e relatively prime *)

exception FoundProof of  ProofFormat.prf_rule

(** [check_sat]
    - detects constraints that are not satisfiable;
    - normalises constraints and generate cuts.
*)

let check_int_sat (cstr,prf) =
  let {coeffs=coeffs ; op=op ; cst=cst} = cstr in
  match Vect.choose coeffs with
  | None ->
     if eval_op op (Int 0) cst then Tauto else Unsat prf
  | _  ->
     let gcdi =  Vect.gcd coeffs in
     let gcd = Big_int gcdi in
     if eq_num gcd (Int 1)
     then Normalise(cstr,prf)
     else
       if Int.equal (sign_num (mod_num cst gcd)) 0
       then (* We can really normalise *)
         begin
           assert (sign_num gcd >=1 ) ;
           let cstr = {
               coeffs = Vect.div gcd  coeffs;
               op = op ; cst = cst // gcd
             } in
           Normalise(cstr,ProofFormat.Gcd(gcdi,prf))
                    (*     Normalise(cstr,CutPrf prf)*)
         end
       else
         match op with
         | Eq -> Unsat (ProofFormat.CutPrf prf)
         | Ge ->
            let cstr = {
                coeffs = Vect.div gcd coeffs;
                op = op ; cst = ceiling_num (cst // gcd)
              } in Cut(cstr,ProofFormat.CutPrf prf)
         | Gt -> failwith "check_sat : Unexpected operator"

let apply_and_normalise check f psys =
  List.fold_left (fun acc pc' ->
      match f pc' with
      | None -> pc'::acc
      | Some pc' ->
         match check pc' with
         | Tauto -> acc
         | Unsat prf -> raise (FoundProof prf)
         | Cut(c,p)  -> (c,p)::acc
         | Normalise (c,p) -> (c,p)::acc
    ) [] psys

let is_linear_for v pc =
  LinPoly.is_linear (fst (fst pc)) || LinPoly.is_linear_for v (fst (fst pc))

(*let non_linear_pivot sys pc v pc' =
  if LinPoly.is_linear (fst (fst pc'))
  then None (* There are other ways to deal with those *)
  else WithProof.linear_pivot sys pc v pc'
*)

let is_linear_substitution sys ((p,o),prf) =
  let pred v =  v =/ Int 1 || v =/ Int (-1)  in
  match o with
  | Eq -> begin
      match
        List.filter (fun v -> List.for_all (is_linear_for v) sys) (LinPoly.search_all_linear pred p)
      with
      | [] -> None
      | v::_  -> Some v (* make a choice *)
    end
  | _  -> None

let elim_simple_linear_equality sys0 =

  let elim sys =
    let (oeq,sys') = extract (is_linear_substitution sys) sys in
    match oeq with
    | None -> None
    | Some(v,pc) -> simplify (WithProof.linear_pivot sys0 pc v) sys' in

  iterate_until_stable elim sys0

let output_sys o sys =
  List.iter (fun s -> Printf.fprintf o "%a\n" WithProof.output s) sys

let subst sys =
  let sys' = WithProof.subst sys in
  if debug then Printf.fprintf stdout "[subst:\n%a\n==>\n%a\n]" output_sys sys output_sys sys' ;
  sys'

(** [saturate_linear_equality sys] generate new constraints
    obtained by eliminating linear equalities by pivoting.
    For integers, the obtained constraints are  sound but not complete.
*)
  let saturate_by_linear_equalities sys0 =
    WithProof.saturate_subst false sys0

let saturate_by_linear_equalities sys =
  let sys' = saturate_by_linear_equalities sys in
  if debug then Printf.fprintf stdout "[saturate_by_linear_equalities:\n%a\n==>\n%a\n]" output_sys sys output_sys sys' ;
  sys'

(* let saturate_linear_equality_non_linear sys0 =
  let (l,_) = extract_all (is_substitution false) sys0 in
  let rec elim l acc =
    match l with
    | [] -> acc
    | (v,pc)::l' ->
       let nc = saturate (non_linear_pivot sys0 pc v) (sys0@acc) in
       elim l' (nc@acc) in
  elim l []
*)

let bounded_vars (sys: WithProof.t list) =
  let l =  (fst (extract_all (fun ((p,o),prf) ->
                     LinPoly.is_variable p
                   ) sys)) in
  List.fold_left (fun acc (i,wp) -> IMap.add i wp acc) IMap.empty l

let rec power n p =
  if n = 1 then p
  else WithProof.product p (power (n-1) p)

let bound_monomial mp m =
  if  Monomial.is_var m || Monomial.is_const m
  then None
  else
     try
       Some (Monomial.fold
               (fun v i acc ->
                 let wp = IMap.find v mp in
                 WithProof.product (power i wp) acc) m (WithProof.const (Int 1))
         )
     with Not_found -> None

let bound_monomials (sys:WithProof.t list) =
  let mp = bounded_vars sys in
  let m  =
    List.fold_left (fun acc ((p,_),_) ->
        Vect.fold (fun acc v _ -> let m = LinPoly.MonT.retrieve v in
                                  match bound_monomial mp m with
                                  | None -> acc
                                  | Some r -> IMap.add v r acc) acc p) IMap.empty sys in
  IMap.fold (fun _ e acc -> e::acc) m []

let develop_constraints prfdepth n_spec sys =
  LinPoly.MonT.clear ();
  max_nb_cstr := compute_max_nb_cstr sys prfdepth ;
  let sys = List.map (develop_constraint n_spec) sys in
  List.mapi (fun i (p,o) -> ((LinPoly.linpol_of_pol p,o),ProofFormat.Hyp i)) sys

let square_of_var i =
  let x = LinPoly.var i in
  ((LinPoly.product x x,Ge),(ProofFormat.Square  x))


(** [nlinear_preprocess  sys]  augments the system [sys] by performing some limited non-linear reasoning.
    For instance, it asserts that the x² ≥0 but also that if c₁ ≥ 0 ∈ sys and c₂ ≥ 0 ∈ sys then c₁ × c₂ ≥ 0.
    The resulting system is linearised.
*)

let nlinear_preprocess  (sys:WithProof.t list) =

  let is_linear =  List.for_all (fun ((p,_),_) -> LinPoly.is_linear p) sys in

  if is_linear then sys
  else
    let collect_square =
      List.fold_left (fun acc ((p,_),_) -> MonMap.union (fun k e1 e2 -> Some e1) acc (LinPoly.collect_square p)) MonMap.empty sys in
    let sys = MonMap.fold (fun s m acc ->
                  let s = LinPoly.of_monomial s in
                  let m = LinPoly.of_monomial m in
                  ((m, Ge), (ProofFormat.Square s))::acc) collect_square  sys in

    let collect_vars = List.fold_left (fun acc p -> ISet.union acc (LinPoly.variables (fst (fst p)))) ISet.empty sys in

    let sys = ISet.fold (fun i acc -> square_of_var i :: acc) collect_vars sys  in

    let sys = sys @ (all_pairs WithProof.product sys) in

    if debug then begin
        Printf.fprintf stdout "Preprocessed\n";
        List.iter (fun s -> Printf.fprintf stdout "%a\n" WithProof.output s) sys;
      end ;

    List.map (WithProof.annot "P") sys

let nlinear_prover prfdepth sys =
  let sys = develop_constraints prfdepth q_spec sys in
  let sys1 = elim_simple_linear_equality sys in
  let sys2 = saturate_by_linear_equalities sys1 in
  let sys = nlinear_preprocess sys1@sys2 in
  let sys = List.map (fun ((p,o),prf) -> (cstr_of_poly (p,o), prf)) sys in
  let id  = (List.fold_left
              (fun acc (_,r) -> max acc (ProofFormat.pr_rule_max_id r)) 0 sys) in
  let env = CList.interval 0 id in
  match linear_prover_cstr sys with
  | None -> Unknown
  | Some cert ->
     Prf (ProofFormat.cmpl_prf_rule Mc.normQ CamlToCoq.q env cert)

let linear_prover_with_cert prfdepth sys =
  let sys = develop_constraints prfdepth q_spec sys in
  (*  let sys = nlinear_preprocess  sys in *)
  let sys = List.map (fun (c,p) -> cstr_of_poly c,p) sys in

  match linear_prover_cstr sys with
  | None -> Unknown
  | Some cert ->
     Prf (ProofFormat.cmpl_prf_rule Mc.normQ CamlToCoq.q (List.mapi (fun i e -> i) sys) cert)

(* The prover is (probably) incomplete --
   only searching for naive cutting planes *)

open Sos_types

let rec scale_term t =
  match t with
  | Zero    -> unit_big_int , Zero
  | Const n ->  (denominator n) , Const (Big_int (numerator n))
  | Var n   -> unit_big_int , Var n
  | Opp t   -> let s, t = scale_term t in s, Opp t
  | Add(t1,t2) -> let s1,y1 = scale_term t1 and s2,y2 = scale_term t2 in
                  let g = gcd_big_int s1 s2 in
                  let s1' = div_big_int s1 g in
                  let s2' = div_big_int s2 g in
                  let e = mult_big_int g (mult_big_int s1' s2') in
                  if Int.equal (compare_big_int e unit_big_int) 0
                  then (unit_big_int, Add (y1,y2))
                  else  e, Add (Mul(Const (Big_int s2'), y1),
                  Mul (Const (Big_int s1'), y2))
  | Sub _ -> failwith "scale term: not implemented"
  | Mul(y,z) ->       let s1,y1 = scale_term y and s2,y2 = scale_term z in
                      mult_big_int s1 s2 , Mul (y1, y2)
  |  Pow(t,n) -> let s,t = scale_term t in
                 power_big_int_positive_int s  n , Pow(t,n)

let scale_term t =
  let (s,t') = scale_term t in
  s,t'

let rec scale_certificate pos = match pos with
  | Axiom_eq i ->  unit_big_int , Axiom_eq i
  | Axiom_le i ->  unit_big_int , Axiom_le i
  | Axiom_lt i ->   unit_big_int , Axiom_lt i
  | Monoid l   -> unit_big_int , Monoid l
  | Rational_eq n ->  (denominator n) , Rational_eq (Big_int (numerator n))
  | Rational_le n ->  (denominator n) , Rational_le (Big_int (numerator n))
  | Rational_lt n ->  (denominator n) , Rational_lt (Big_int (numerator n))
  | Square t -> let s,t' = scale_term t in
                mult_big_int s s , Square t'
  | Eqmul (t, y) -> let s1,y1 = scale_term t and s2,y2 = scale_certificate y in
                    mult_big_int s1 s2 , Eqmul (y1,y2)
  | Sum (y, z) -> let s1,y1 = scale_certificate y
                  and s2,y2 = scale_certificate z in
                  let g = gcd_big_int s1 s2 in
                  let s1' = div_big_int s1 g in
                  let s2' = div_big_int s2 g in
                  mult_big_int g (mult_big_int s1' s2'),
                  Sum (Product(Rational_le (Big_int s2'), y1),
                       Product (Rational_le (Big_int s1'), y2))
  | Product (y, z) ->
     let s1,y1 = scale_certificate y and s2,y2 = scale_certificate z in
     mult_big_int s1 s2 , Product (y1,y2)

open Micromega
let rec term_to_q_expr = function
  | Const n ->  PEc (Ml2C.q n)
  | Zero   ->  PEc ( Ml2C.q (Int 0))
  | Var s   ->  PEX (Ml2C.index
                       (int_of_string (String.sub s 1 (String.length s - 1))))
  | Mul(p1,p2) ->  PEmul(term_to_q_expr p1, term_to_q_expr p2)
  | Add(p1,p2) ->   PEadd(term_to_q_expr p1, term_to_q_expr p2)
  | Opp p ->   PEopp (term_to_q_expr p)
  | Pow(t,n) ->  PEpow (term_to_q_expr t,Ml2C.n n)
  | Sub(t1,t2) ->  PEsub (term_to_q_expr t1,  term_to_q_expr t2)

let term_to_q_pol e = Mc.norm_aux (Ml2C.q (Int 0)) (Ml2C.q (Int 1)) Mc.qplus  Mc.qmult Mc.qminus Mc.qopp Mc.qeq_bool (term_to_q_expr e)

let rec product l =
  match l with
  | [] -> Mc.PsatzZ
  | [i] -> Mc.PsatzIn (Ml2C.nat i)
  | i ::l -> Mc.PsatzMulE(Mc.PsatzIn (Ml2C.nat i), product l)

let  q_cert_of_pos  pos =
  let rec _cert_of_pos = function
      Axiom_eq i ->  Mc.PsatzIn (Ml2C.nat i)
    | Axiom_le i ->  Mc.PsatzIn (Ml2C.nat i)
    | Axiom_lt i ->  Mc.PsatzIn (Ml2C.nat i)
    | Monoid l  -> product l
    | Rational_eq n | Rational_le n | Rational_lt n ->
       if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ else
         Mc.PsatzC (Ml2C.q   n)
    | Square t -> Mc.PsatzSquare (term_to_q_pol  t)
    | Eqmul (t, y) -> Mc.PsatzMulC(term_to_q_pol t, _cert_of_pos y)
    | Sum (y, z) -> Mc.PsatzAdd  (_cert_of_pos y, _cert_of_pos z)
    | Product (y, z) -> Mc.PsatzMulE (_cert_of_pos y, _cert_of_pos z) in
  simplify_cone q_spec (_cert_of_pos pos)

let rec term_to_z_expr = function
  | Const n ->  PEc (Ml2C.bigint (big_int_of_num n))
  | Zero   ->  PEc ( Z0)
  | Var s   ->  PEX (Ml2C.index
                       (int_of_string (String.sub s 1 (String.length s - 1))))
  | Mul(p1,p2) ->  PEmul(term_to_z_expr p1, term_to_z_expr p2)
  | Add(p1,p2) ->   PEadd(term_to_z_expr p1, term_to_z_expr p2)
  | Opp p ->   PEopp (term_to_z_expr p)
  | Pow(t,n) ->  PEpow (term_to_z_expr t,Ml2C.n n)
  | Sub(t1,t2) ->  PEsub (term_to_z_expr t1,  term_to_z_expr t2)

let term_to_z_pol e = Mc.norm_aux (Ml2C.z 0) (Ml2C.z 1) Mc.Z.add  Mc.Z.mul Mc.Z.sub Mc.Z.opp Mc.zeq_bool (term_to_z_expr e)

let  z_cert_of_pos  pos =
  let s,pos = (scale_certificate pos) in
  let rec _cert_of_pos = function
      Axiom_eq i ->  Mc.PsatzIn (Ml2C.nat i)
    | Axiom_le i ->  Mc.PsatzIn (Ml2C.nat i)
    | Axiom_lt i ->  Mc.PsatzIn (Ml2C.nat i)
    | Monoid l  -> product l
    | Rational_eq n | Rational_le n | Rational_lt n ->
       if Int.equal (compare_num n (Int 0)) 0 then Mc.PsatzZ else
         Mc.PsatzC (Ml2C.bigint (big_int_of_num  n))
    | Square t -> Mc.PsatzSquare (term_to_z_pol  t)
    | Eqmul (t, y) ->
       let is_unit =
         match t with
         | Const n -> n =/ Int 1
         |   _     -> false in
       if is_unit
       then _cert_of_pos y
       else Mc.PsatzMulC(term_to_z_pol t, _cert_of_pos y)
    | Sum (y, z) -> Mc.PsatzAdd  (_cert_of_pos y, _cert_of_pos z)
    | Product (y, z) -> Mc.PsatzMulE (_cert_of_pos y, _cert_of_pos z) in
  simplify_cone z_spec (_cert_of_pos pos)

(** All constraints (initial or derived) have an index and have a justification i.e., proof.
    Given a constraint, all the coefficients are always integers.
*)
open Mutils
open Num
open Big_int
open Polynomial

type prf_sys = (cstr * ProofFormat.prf_rule) list

(** Proof generating pivoting over variable v *)
let pivot v (c1,p1) (c2,p2) =
  let {coeffs = v1 ; op = op1 ; cst = n1} = c1
  and {coeffs = v2 ; op = op2 ; cst = n2} = c2 in

  (* Could factorise gcd... *)
  let xpivot cv1 cv2 =
    (
      {coeffs = Vect.add (Vect.mul cv1 v1) (Vect.mul cv2 v2) ;
       op = opAdd op1 op2 ;
       cst = n1 */ cv1 +/ n2 */ cv2 },

      ProofFormat.add_proof (ProofFormat.mul_cst_proof  cv1 p1) (ProofFormat.mul_cst_proof  cv2 p2)) in

  match Vect.get v v1 , Vect.get v v2 with
  | Int 0 , _ | _ , Int 0 -> None
  |  a   ,  b   ->
      if Int.equal ((sign_num a) * (sign_num b)) (-1)
      then
        let cv1 = abs_num b
        and cv2 = abs_num a  in
        Some (xpivot cv1 cv2)
      else
        if op1 == Eq
        then
          let cv1 = minus_num (b */ (Int (sign_num a)))
          and cv2 = abs_num a in
          Some (xpivot cv1 cv2)
        else if op2 == Eq
        then
          let cv1 = abs_num b
          and cv2 = minus_num (a */ (Int  (sign_num b))) in
          Some (xpivot cv1 cv2)
        else  None (* op2 could be Eq ... this might happen *)

let simpl_sys sys =
  List.fold_left (fun acc (c,p) ->
      match check_int_sat (c,p) with
      | Tauto -> acc
      | Unsat prf -> raise (FoundProof prf)
      | Cut(c,p)  -> (c,p)::acc
      | Normalise (c,p) -> (c,p)::acc) [] sys

(** [ext_gcd a b] is the extended Euclid algorithm.
    [ext_gcd a b = (x,y,g)] iff [ax+by=g]
    Source: http://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
*)
let rec ext_gcd a b =
  if Int.equal (sign_big_int b) 0
  then (unit_big_int,zero_big_int)
  else
    let (q,r) = quomod_big_int a b in
    let (s,t) = ext_gcd b r in
    (t, sub_big_int s (mult_big_int q t))

let extract_coprime (c1,p1) (c2,p2) =
  if c1.op == Eq && c2.op == Eq
  then Vect.exists2 (fun n1 n2 ->
           Int.equal (compare_big_int (gcd_big_int (numerator n1) (numerator n2)) unit_big_int) 0)
           c1.coeffs c2.coeffs
  else None

let extract2 pred l =
  let rec xextract2 rl l =
    match l with
    | [] -> (None,rl) (* Did not find *)
    | e::l ->
       match extract (pred e) l with
       | None,_ -> xextract2 (e::rl) l
       | Some (r,e'),l' -> Some (r,e,e'), List.rev_append rl l' in

  xextract2 [] l

let extract_coprime_equation psys =
  extract2 extract_coprime psys

let pivot_sys v pc psys = apply_and_normalise check_int_sat (pivot v pc) psys

let reduce_coprime psys =
  let oeq,sys = extract_coprime_equation psys in
  match oeq with
  | None -> None (* Nothing to do *)
  | Some((v,n1,n2),(c1,p1),(c2,p2) ) ->
     let (l1,l2) = ext_gcd (numerator n1) (numerator n2) in
     let l1' = Big_int l1 and l2' = Big_int l2 in
     let cstr =
       {coeffs = Vect.add (Vect.mul l1' c1.coeffs) (Vect.mul l2' c2.coeffs);
        op = Eq ;
        cst = (l1' */ c1.cst) +/ (l2' */ c2.cst)
       } in
     let prf = ProofFormat.add_proof (ProofFormat.mul_cst_proof  l1' p1) (ProofFormat.mul_cst_proof l2' p2) in

     Some (pivot_sys v (cstr,prf) ((c1,p1)::sys))

(** If there is an equation [eq] of the form 1.x + e = c, do a pivot over x with equation [eq] *)
let reduce_unary psys =
  let is_unary_equation (cstr,prf) =
    if cstr.op == Eq
    then
      Vect.find (fun v n -> if n =/ (Int 1) || n=/ (Int (-1)) then Some v else None) cstr.coeffs
    else None in

  let (oeq,sys) =  extract is_unary_equation psys in
  match oeq with
  | None -> None (* Nothing to do *)
  | Some(v,pc) ->
     Some(pivot_sys v pc sys)

let reduce_var_change psys =

  let rec rel_prime vect =
    match Vect.choose vect with
    | None -> None
    | Some(x,v,vect) ->
       let v = numerator v in
       match Vect.find (fun x' v' ->
                           let v' = numerator v' in
                           if eq_big_int (gcd_big_int  v v') unit_big_int
                           then  Some(x',v') else None) vect with
       | Some(x',v') ->  Some ((x,v),(x', v'))
       | None  -> rel_prime vect in

  let rel_prime (cstr,prf) =  if cstr.op == Eq then rel_prime cstr.coeffs else None in

  let (oeq,sys) = extract rel_prime psys in

  match oeq with
  | None -> None
  | Some(((x,v),(x',v')),(c,p)) ->
     let (l1,l2) = ext_gcd  v  v' in
     let l1,l2 = Big_int l1 , Big_int l2 in

     let pivot_eq (c',p') =
       let {coeffs = coeffs ; op = op ; cst = cst} = c' in
       let vx = Vect.get x coeffs in
       let vx' = Vect.get x' coeffs in
       let m = minus_num (vx */ l1 +/ vx' */ l2) in
       Some ({coeffs =
                Vect.add (Vect.mul m c.coeffs) coeffs ; op = op ; cst = m */ c.cst +/ cst} ,
             ProofFormat.add_proof (ProofFormat.mul_cst_proof m p) p') in

     Some (apply_and_normalise check_int_sat pivot_eq sys)

let reduction_equations psys =
  iterate_until_stable (app_funs
                          [reduce_unary ; reduce_coprime ;
                           reduce_var_change (*; reduce_pivot*)]) psys

(** [get_bound sys] returns upon success an interval (lb,e,ub) with proofs *)
let get_bound sys =
  let is_small (v,i) =
    match Itv.range i with
    | None -> false
    | Some i -> i <=/ (Int 1) in

  let select_best (x1,i1) (x2,i2) =
    if Itv.smaller_itv i1 i2
    then (x1,i1) else (x2,i2) in

  (* For lia, there are no equations => these precautions are not needed *)
  (* For nlia, there are equations => do not enumerate over equations! *)
  let all_planes sys =
    let (eq,ineq) = List.partition (fun c -> c.op == Eq) sys in
    match eq with
    | [] -> List.rev_map (fun c -> c.coeffs) ineq
    | _  ->
       List.fold_left (fun acc c ->
           if List.exists (fun c' -> Vect.equal c.coeffs c'.coeffs) eq
           then acc else c.coeffs ::acc) [] ineq in

  let smallest_interval =
    List.fold_left
      (fun acc vect ->
        if is_small acc
        then acc
        else
          match optimise vect sys with
          | None -> acc
          | Some i ->
             if debug then Printf.printf "Found a new bound %a in %a" Vect.pp vect Itv.pp i;
             select_best (vect,i) acc)  (Vect.null, (None,None)) (all_planes sys) in
  let smallest_interval =
    match smallest_interval
    with
    | (x,(Some i, Some j))  -> Some(i,x,j)
    |   x        ->   None (* This should not be possible *)
  in
  match smallest_interval with
  | Some (lb,e,ub) ->
     let (lbn,lbd) = (sub_big_int (numerator lb)  unit_big_int, denominator lb) in
     let (ubn,ubd) = (add_big_int unit_big_int (numerator ub) , denominator ub) in
     (match
        (* x <= ub ->  x  > ub *)
        direct_linear_prover   ({coeffs = Vect.mul (Big_int ubd)  e ; op = Ge ; cst = Big_int ubn} :: sys),
        (* lb <= x  -> lb > x *)
        direct_linear_prover
          ({coeffs = Vect.mul (minus_num (Big_int lbd)) e ; op = Ge ; cst = minus_num (Big_int lbn)} :: sys)
      with
      | Some cub , Some clb  -> Some(List.tl (Vect.to_list clb),(lb,e,ub), List.tl (Vect.to_list cub))
      |         _            -> failwith "Interval without proof"
     )
  | None -> None

let check_sys sys =
  List.for_all (fun (c,p) -> Vect.for_all (fun _ n -> sign_num n <> 0) c.coeffs) sys

open ProofFormat

let xlia (can_enum:bool)  reduction_equations  sys =

  let rec enum_proof (id:int) (sys:prf_sys)  =
    if debug then (Printf.printf "enum_proof\n" ; flush stdout) ;
    assert (check_sys sys) ;

    let nsys,prf = List.split sys in
    match get_bound nsys with
    | None -> Unknown (* Is the systeme really unbounded ? *)
    | Some(prf1,(lb,e,ub),prf2) ->
       if debug then Printf.printf "Found interval: %a in [%s;%s] -> " Vect.pp e (string_of_num lb) (string_of_num ub) ;
       (match start_enum  id  e  (ceiling_num lb)  (floor_num ub) sys
        with
        | Prf prfl ->
           Prf(ProofFormat.Enum(id,ProofFormat.proof_of_farkas (env_of_list prf) (Vect.from_list prf1),e,
                     ProofFormat.proof_of_farkas (env_of_list prf) (Vect.from_list prf2),prfl))
        | _ -> Unknown
       )

  and start_enum id e clb cub sys  =
    if clb >/ cub
    then Prf []
    else
      let eq = {coeffs = e ; op = Eq ; cst = clb} in
      match aux_lia (id+1) ((eq, ProofFormat.Def id) :: sys) with
      | Unknown | Model _ -> Unknown
      | Prf prf  ->
         match start_enum id e (clb +/ (Int 1)) cub sys with
         | Prf l -> Prf (prf::l)
         | _ -> Unknown

  and aux_lia (id:int)  (sys:prf_sys)   =
    assert (check_sys sys) ;
    if debug then Printf.printf "xlia: %a \n" (pp_list ";" (fun o (c,_) -> output_cstr o c)) sys ;
    try
      let sys = reduction_equations sys in
      if debug then
        Printf.printf "after reduction: %a \n" (pp_list ";" (fun o (c,_) -> output_cstr o c)) sys ;
      match linear_prover_cstr sys with
      | Some prf -> Prf (Step(id,prf,Done))
      | None ->  if can_enum then enum_proof id sys else Unknown
    with FoundProof prf ->
      (* [reduction_equations] can find a proof *)
      Prf(Step(id,prf,Done)) in

  (*  let sys' = List.map (fun (p,o) -> Mc.norm0 p , o) sys in*)
  let id  = 1 + (List.fold_left
              (fun acc (_,r) -> max acc (ProofFormat.pr_rule_max_id r)) 0 sys) in
  let orpf =
    try
      let sys = simpl_sys sys in
      aux_lia id sys
    with FoundProof pr -> Prf(Step(id,pr,Done)) in
  match orpf with
  | Unknown | Model _ -> Unknown
  | Prf prf ->
     let env = CList.interval 0 (id - 1) in
     if debug then begin
         Printf.fprintf stdout "direct proof %a\n" output_proof prf;
         flush stdout;
       end;
       let prf = compile_proof env prf in
     (*try
     if Mc.zChecker sys' prf then Some prf else
     raise Certificate.BadCertificate
     with Failure s -> (Printf.printf "%s" s ; Some prf)
      *)

let xlia_simplex env red sys =
  let compile_prf sys prf =
    let id  = 1 + (List.fold_left
                          (fun acc (_,r) -> max acc (ProofFormat.pr_rule_max_id r)) 0 sys) in
    let env = CList.interval 0 (id - 1) in
    Prf (compile_proof env prf) in

  try
    let sys = red sys in

    match Simplex.integer_solver sys with
    | None -> Unknown
    | Some prf -> compile_prf sys prf
  with FoundProof prf -> compile_prf sys (Step(0,prf,Done))

let xlia env0 en red sys =
  if !use_simplex then xlia_simplex env0 red sys
  else xlia en red sys

let dump_file = ref None

let gen_bench (tac, prover) can_enum prfdepth sys =
  let res = prover can_enum prfdepth sys in
  (match !dump_file with
  | None -> ()
  | Some file ->
     begin
       let o = open_out (Filename.temp_file ~temp_dir:(Sys.getcwd ()) file ".v") in
       let sys = develop_constraints prfdepth z_spec sys in
       Printf.fprintf o "Require Import ZArith Lia. Open Scope Z_scope.\n";
       Printf.fprintf o "Goal %a.\n" (LinPoly.pp_goal "Z") (List.map fst sys) ;
       begin
         match res with
         | Unknown | Model _ ->
            Printf.fprintf o "Proof.\n intros. Fail %s.\nAbort.\n" tac
         | Prf res ->
            Printf.fprintf o "Proof.\n intros. %s.\nQed.\n" tac
       end
       ;
         flush o ;
       close_out o ;
     end);
  res

let lia (can_enum:bool) (prfdepth:int) sys =
  let sys = develop_constraints prfdepth z_spec sys in
  if debug then begin
      Printf.fprintf stdout "Input problem\n";
      List.iter (fun s -> Printf.fprintf stdout "%a\n" WithProof.output s) sys;
      Printf.fprintf stdout "Input problem\n";
      let string_of_op = function Eq -> "=" | Ge -> ">=" | Gt -> ">" in
      List.iter (fun ((p,op),_) -> Printf.fprintf stdout "(assert (%s %a))\n" (string_of_op op) Vect.pp_smt p) sys;
    end;
  let sys = subst sys in
  let bnd = bound_monomials sys in (* To deal with non-linear monomials *)
  let sys = bnd@(saturate_by_linear_equalities sys)@sys in

  let sys' = List.map (fun ((p,o),prf) -> (cstr_of_poly (p,o), prf)) sys in
  xlia (List.map fst sys) can_enum reduction_equations sys'

let make_cstr_system sys =
  List.map (fun ((p,o),prf) -> (cstr_of_poly (p,o), prf)) sys

let nlia enum prfdepth sys =
  let sys = develop_constraints prfdepth z_spec sys in
  let is_linear =  List.for_all (fun ((p,_),_) -> LinPoly.is_linear p) sys in

  if debug then begin
      Printf.fprintf stdout "Input problem\n";
      List.iter (fun s -> Printf.fprintf stdout "%a\n" WithProof.output s) sys;
    end;

  if is_linear
  then xlia (List.map fst sys) enum reduction_equations (make_cstr_system sys)
  else
    (*
      let sys1 = elim_every_substitution sys in
      No: if a wrong equation is chosen, the proof may fail.
      It would only be safe if the variable is linear...
     *)
    let sys1 = elim_simple_linear_equality sys in
    let sys2 = saturate_by_linear_equalities sys1 in
    let sys3 = nlinear_preprocess (sys1@sys2) in

    let sys4  = make_cstr_system ((*sys2@*)sys3) in
    (* [reduction_equations] is too brutal - there should be some non-linear reasoning  *)
    xlia (List.map fst sys) enum  reduction_equations sys4

(* For regression testing, if bench = true generate a Coq goal *)

let lia can_enum prfdepth sys =
  gen_bench ("lia",lia) can_enum prfdepth sys

let nlia enum prfdepth sys =
  gen_bench ("nia",nlia) enum prfdepth sys

(* Local Variables: *)
(* coding: utf-8 *)
(* End: *)

¤ Dauer der Verarbeitung: 0.23 Sekunden (vorverarbeitet) ¤

Download des Quellennavigators
Download des sprechenden Kalenders
in der Quellcodebibliothek suchen

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.