package coq-core

  1. Overview
  2. Docs
package coq-core
Legend:
Page
Library
Module
Module type
Parameter
Class
Class type
Source

Source file linsolve.ml

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
(************************************************************************)
(*         *   The Coq Proof Assistant / The Coq Development Team       *)
(*  v      *         Copyright INRIA, CNRS and contributors             *)
(* <O___,, * (see version control and CREDITS file for authors & dates) *)
(*   \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)         *)
(************************************************************************)


(** An equation is of the following form a1.x1 + a2.x2 = c
*)

type var = int
type id  = int

module Itv =
struct
  let debug = false

  type t = int * int (* We only consider closed intervals *)

  exception Empty

  let output o (lb,ub) =
    Printf.fprintf o "[%i,%i]" lb ub

  let wf o (lb,ub) =
    if lb <= ub then ()
    else Printf.fprintf o "error %a\n" output (lb,ub)

  (** [mul_cst c (lb,ub)] requires c > 0 *)
  let mul_cst c (lb,ub) = (c * lb, c * ub)

  (** [opp (lb,ub)] is multplication by -1 *)
  let opp (lb,ub) = (-ub,-lb)

  let opp i1  =
    let i = opp i1 in
    if debug then Printf.printf  "opp %a -> %a\n" output i1  output i;
      i

  (** [div (lb,ub) c] requires c > 0, lb >= 0 *)
  let div (lb,ub) c =
    let lb = lb /c + (if lb  mod c = 0 then 0 else 1) in
    let ub = ub / c in
    if lb <= ub then (lb,ub) else raise Empty

  let div i c =
    try let r = div i c in
      if debug then Printf.printf  "%a div %i -> %a\n" output i c output r;
    r
    with Empty ->
      if debug then Printf.printf  "%a div %i -> Empty \n" output i c;
      raise Empty

  let add (lb1,ub1) (lb2,ub2) =
    (lb1+lb2,ub1+ub2)

  let add i1 i2 =
      let i = add i1 i2 in
      if debug then Printf.printf  "%a add %a -> %a\n" output i1 output i2 output i;
      i

  let inter : t -> t -> t =
    fun (lb1,ub1) (lb2,ub2) ->
    let ub = max lb1 lb2 in
    let lb = min ub1 ub2 in
    if ub <= lb
    then (ub,lb)
    else raise Empty

  let inter i1 i2 =
    try
      let i = inter i1 i2 in
      if debug then Printf.printf  "%a inter %a -> %a\n" output i1 output i2 output i;
      i
    with Empty ->
      if debug then Printf.printf  "%a inter %a -> Empty\n" output i1 output i2 ;
      raise Empty

  (* [enum (lb,ub)] is only defined for finite intervals *)
  let enum (lb,ub) =
    match Int.compare lb ub with
    | 0 -> (lb,None)
    | _ -> (lb,Some (lb+1,ub))

  let range (lb,ub) = ub - lb + 1

  let top = (min_int,max_int)

  let lt i1 i2 =
    range i1 < range i2

end

module ItvMap =
struct

  module M = Map.Make(Int)

  include M

  let refine_with v i m =
    try
      let i0 = M.find v m in
      let i' = Itv.inter i i0 in
      (Itv.lt i' i0,i',M.add v i' m)
    with Not_found ->
      (true, i, M.add v i m)

  let pick m =
    let (x,i,r) =
      fold (fun v i (v',i',r') ->
          let r = Itv.range i in
          if r < r' then (v,i,r) else (v',i',r')) m (min_int,Itv.top,max_int) in
    if x = min_int then raise Not_found else (x,i)

  let output o m =
    Printf.fprintf o "[";
    iter (fun k (lb,ub) -> Printf.fprintf o "x%i -> [%i,%i] " k lb ub) m;
    Printf.fprintf o "]";


end

exception Unsat

module Eqn =
struct
  type t   = (var * int) list * int

  let empty = ([],0)

  let rec output_lin o l =
    match l with
    | [] -> Printf.fprintf o "0"
    | [x,v] -> Printf.fprintf o "%i.x%i" v x
    | (x,v)::l' -> Printf.fprintf o "%i.x%i + %a" v x output_lin l'

  let normalise (l,c) =
    match l with
    | [] -> if c = 0 then None else raise Unsat
    | _  -> Some(l,c)


  let rec no_dup l =
    match l with
    | [] -> true
    | (x,v)::l -> try
        let _ = List.assoc x l in
        false
      with Not_found -> no_dup l


  let add (l1,c1) (l2,c2) =
    (l1@l2,c1+c2)

  let add e1 e2 =
    let r = add e1 e2 in
    if no_dup (fst r) then ()
    else Printf.printf "add(duplicate)%a %a" output_lin (fst e1) output_lin (fst e2) ;
    r


  let itv_of_ax m (var,coe) =
    Itv.mul_cst coe (ItvMap.find var m)

  let itv_list m l =
    List.fold_left (fun i (var,coe) ->
        Itv.add i (itv_of_ax m (var,coe))) (0,0) l

  let get_remove x l =
    let l' = List.remove_assoc x l in
    let c  = try List.assoc x l with Not_found -> 0 in
    (c,l')



end
type eqn = Eqn.t

open Eqn

let debug = false

(** Given an equation a1.x1 + ... an.xn = c,
    bound all the variables xi in [0; c/ai] *)

let init_bound m (v,c) =
  match v with
  | [] -> if c = 0 then m else raise Unsat
  | [x,v] -> let (_,_,m) = ItvMap.refine_with x (Itv.div (c,c) v) m in m
  |   _   ->
    List.fold_left (fun m (var,coe)  ->
        let (_,_,m) = ItvMap.refine_with var (0,c / coe) m in
        m)  m v

let init_bounds sys =
  List.fold_left init_bound ItvMap.empty sys

let init_bounds sys =
  let m = init_bounds sys in
  if debug then Printf.printf "init_bound : %a\n" ItvMap.output m;
  m


(* [refine_bound p m acc (v,c)]
    improves the bounds of the equation v + acc = c
*)


let rec refine_bound p m acc (v,c) =
  Itv.wf stdout acc;
  match v with
  | [] -> (m,p)
  | (var,coe)::v' ->
    if debug then Printf.printf "Refining %i.x%i + %a + %a = %i\n" coe var
      Itv.output acc output_lin v' c;
    let itv_acc_l = Itv.inter (0,c) (Itv.add acc (itv_list m v')) in
    let itv_coe_var =  Itv.add (c,c) (Itv.opp itv_acc_l) in
    let i      =  Itv.div itv_coe_var coe in
    let (b,i',m) = ItvMap.refine_with var i m in
    refine_bound (p || b)
      m (Itv.add (Itv.mul_cst coe i') acc) (v',c)

let refine_bounds p m l =
  List.fold_left (fun (m,p) eqn ->
      refine_bound p m (0,0) eqn) (m,p) l

let refine_until_fix m l =
  let rec iter_refine m =
    let (m',b) = refine_bounds false m l in
    if b then iter_refine m' else m' in
  iter_refine m




let subst x a l =

  let subst_eqn acc (v,c) =
    let (coe,v') = Eqn.get_remove x v in
    let (v',c') = (v', c - coe * a) in
    match v' with
    | [] -> if c' = 0 then acc
      else raise Unsat
    | _ -> (v',c')::acc in

  List.fold_left subst_eqn [] l

let output_list elt o l =
  Printf.fprintf o "[";
  List.iter (fun e -> Printf.fprintf o "%a; " elt e) l;
  Printf.fprintf o "]"

let output_equations o l =
    let output_equation o (l,c) =
        Printf.fprintf o "%a = %i" output_lin l c in
    output_list output_equation o l

let output_intervals o m =
  ItvMap.iter (fun k v -> Printf.fprintf o "x%i:%a " k Itv.output v) m


type solution = (var * int) list

let solve_system l =

  let rec solve m l =
    if debug then Printf.printf "Solve %a\n" output_equations l;

    match l with
    | [] -> [m] (* we have a solution *)
    | _  ->
      try
        let m' = refine_until_fix m l in
        try
        if debug then Printf.printf "Refined %a\n" ItvMap.output m' ;
        let (k,i) = ItvMap.pick m' in
        let (v,itv') = Itv.enum i in
        (* We recursively solve using k = v *)
        let sol1 =
          List.map (ItvMap.add k (v,v))
            (solve (ItvMap.remove k m) (subst k v l)) in
        let sol2 =
          match itv' with
          | None -> []
          | Some itv' ->
            (* We recursively solve with a smaller interval *)
            solve (ItvMap.add k itv' m) l in
        sol1 @ sol2
        with    | Not_found -> Printf.printf "NOT FOUND %a %a\n" output_equations l output_intervals m'; raise Not_found
      with (Unsat | Itv.Empty) as e ->
        begin
          if debug then Printf.printf "Unsat detected %s\n" (Printexc.to_string e);
          [] end
  in

  try
    let l = CList.map_filter Eqn.normalise l in
    solve (init_bounds l) l
  with Itv.Empty | Unsat -> []



let enum_sol m =
  let rec augment_sols x (lb,ub) s =
    let slb = if lb = 0 then s
      else List.rev_map (fun s -> (x,lb)::s) s in
    if lb = ub then slb
    else let sl = augment_sols x (lb+1,ub) s in
      List.rev_append slb sl in

  ItvMap.fold augment_sols m [[]]

let enum_sols l =
  List.fold_left (fun s m -> List.rev_append (enum_sol m) s) [] l


let solve_and_enum l = enum_sols (solve_system l)


let output_solution o s =
  let output_var_coef o (x,v) =
    Printf.fprintf o "x%i:%i" x v in
  output_list output_var_coef o s ;
  Printf.fprintf o "\n"

let output_solutions o l = output_list output_solution o l

(** Incremental construction of systems of equations *)
open Mutils

type system = Eqn.t IMap.t

let empty : system = IMap.empty

let set_constant (idx:int) (c:int) (s:system) : Eqn.t =
  let e = try IMap.find idx s with
    |Not_found -> Eqn.empty in
  (fst e,c)

let make_mon (idx:int) (v:var) (c:int) (s:system) : system  =
  IMap.add idx ([v,c],0) s

let merge (s1:system) (s2:system) : system =
  IMap.merge (fun k e1 e2 ->
      match e1 , e2 with
      | None , None -> None
      | None , Some e | Some e , None -> Some e
      | Some e1, Some e2 -> Some (Eqn.add e1 e2)) s1 s2
OCaml

Innovation. Community. Security.

GitHub Discord Twitter Peertube RSS
About
Changelog Releases Industrial Users Academic Users Why OCaml
Ecosystem
Packages Community Events OCaml Planet Jobs
Resources
Install OCaml Get Started Platform Tools Language Manual Standard Library API Books Exercises Papers OCaml Playground Logo
Policies
Carbon Footprint Governance Privacy Code of Conduct