package codex

  1. Overview
  2. Docs
package codex
Legend:
Page
Library
Module
Module type
Parameter
Class
Class type
Source

Source file float_interval.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
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
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
(**************************************************************************)
(*  This file is part of the Codex semantics library.                     *)
(*                                                                        *)
(*  Copyright (C) 2013-2025                                               *)
(*    CEA (Commissariat à l'énergie atomique et aux énergies              *)
(*         alternatives)                                                  *)
(*                                                                        *)
(*  you can redistribute it and/or modify it under the terms of the GNU   *)
(*  Lesser General Public License as published by the Free Software       *)
(*  Foundation, version 2.1.                                              *)
(*                                                                        *)
(*  It is distributed in the hope that it will be useful,                 *)
(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *)
(*  GNU Lesser General Public License for more details.                   *)
(*                                                                        *)
(*  See the GNU Lesser General Public License version 2.1                 *)
(*  for more details (enclosed in the file LICENSE).                      *)
(*                                                                        *)
(**************************************************************************)

module Comp = Abstract_interp.Comp
open Bottom.Type

type round = Float_sig.round = Up | Down | Near | Zero
type prec = Float_sig.prec = Single | Double | Long_Double | Real

module Make (F: Float_sig.S) = struct

  (* Definitions of useful floating-point constants. *)
  module Cst = struct
    let pos_zero prec = F.of_float Near prec 0.
    let neg_zero prec = F.of_float Near prec (- 0.)
    let pos_infinity prec = F.of_float Near prec infinity
    let neg_infinity prec = F.of_float Near prec neg_infinity

    (* Maximum value M such that all integers 0 < n <= M are exactly
       representable in the given precision.
       It is 2^53 for double and 2^24 for single precision. *)
    let max_precise_integer prec =
      let f = match prec with
        | Single -> 2. ** 24.
        | Double | Long_Double -> 2. ** 53.
        | Real -> infinity
      in
      F.of_float Near prec f

    (* Maximum odd integer representable in the given precision.
       It is 2^53-1 for double and 2^24-1 for single precision. *)
    let max_odd_integer prec =
      let f = match prec with
        | Single -> 2. ** 24. -. 1.
        | Double -> 2. ** 53. -. 1.
        | Real | Long_Double -> infinity
      in
      F.of_float Near prec f

    (* Successor of the maximum non-integer representable in the given precision.
       It is 2^52 in double and 2^23 in single precision. *)
    let ceil_of_max_non_integer prec =
      let f = match prec with
        | Single -> 2. ** 23.
        | Double -> 2. ** 52.
        | Real | Long_Double -> infinity
      in
      F.of_float Near prec f

    (* Returns 1.0 if [f] is positive, and -1.0 if [t] is negative.
       [f] can be infinite or a (positive or negative) zero, but not NaN. *)
    let sign prec f =
      let s = if F.is_negative f then -1. else 1. in
      F.of_float Near prec s
  end

  (* IEEE and non-IEEE comparisons. *)
  module Cmp = struct
    (* Constants in single precision: should be used only for comparisons,
       and not to create intervals of arbitrary precisions. *)
    let pos_zero = Cst.pos_zero Float_sig.Single
    let neg_zero = Cst.neg_zero Float_sig.Single
    let one = F.of_float Near Float_sig.Single 1.
    let pos_infinity = Cst.pos_infinity Float_sig.Single
    let neg_infinity = Cst.neg_infinity Float_sig.Single

    let equal_ieee f1 f2 = F.cmp_ieee f1 f2 = 0
    let le_ieee f1 f2 = F.cmp_ieee f1 f2 <= 0
    let lt_ieee f1 f2 = F.cmp_ieee f1 f2 < 0
    let ge_ieee f1 f2 = F.cmp_ieee f1 f2 >= 0
    let gt_ieee f1 f2 = F.cmp_ieee f1 f2 > 0

    let is_a_zero f = equal_ieee pos_zero f
    let is_pos_infinity f = equal_ieee pos_infinity f
    let is_neg_infinity f = equal_ieee neg_infinity f

    let equal f1 f2 = F.compare f1 f2 = 0
    let le f1 f2 = F.compare f1 f2 <= 0
    let lt f1 f2 = F.compare f1 f2 < 0
    let ge f1 f2 = F.compare f1 f2 >= 0
    let gt f1 f2 = F.compare f1 f2 > 0

    let min f1 f2 = if le f1 f2 then f1 else f2
    let max f1 f2 = if le f1 f2 then f2 else f1

    (* next_float_ieee -0. = next_float_ieee 0. *)
    let next_float_ieee prec f =
      let f = if equal f neg_zero then Cst.pos_zero prec else f in
      F.next_float prec f

    (* prev_float_ieee 0. = next_float_ieee -0. *)
    let prev_float_ieee prec f =
      let f = if equal f pos_zero then Cst.neg_zero prec else f in
      F.prev_float prec f
  end


  module FRange = struct

    (* invariants for intervals Itv (b, e, nan):
       - b and e are not NaN;
       - b <= e *)
    type t =
      | Itv of F.t * F.t * bool
      | NaN

    let nan = NaN

    let inject ?(nan=false) b e =
      if F.is_nan b || F.is_nan e || F.compare b e > 0 then
        Codex_log.fatal "Invalid bounds for float interval:@ %a .. %a@."
          F.pretty b F.pretty e;
      Itv (b, e, nan)

    let inject_after_tighten prec ~nan b e =
      let b = F.round_to_precision Up prec b
      and e = F.round_to_precision Down prec e in
      if Cmp.le b e then `Value (inject ~nan b e)
      else if nan then `Value NaN else `Bottom

    let add_nan = function
      | NaN as f -> f
      | Itv (_, _, true) as f -> f
      | Itv (b, e, false) -> Itv (b, e, true)

    let add_pos_infinity prec t =
      let pos_inf = Cst.pos_infinity prec in
      match t with
      | NaN -> Itv (pos_inf, pos_inf, true)
      | Itv (b, _e, nan) -> Itv (b, pos_inf, nan)

    let add_neg_infinity prec t =
      let neg_inf = Cst.neg_infinity prec in
      match t with
      | NaN -> Itv (neg_inf, neg_inf, true)
      | Itv (_b, e, nan) -> Itv (neg_inf, e, nan)
  end

  (* ------------------------------------------------------------------------
                                      Datatype
     ------------------------------------------------------------------------ *)

  type t = FRange.t

  (* let structural_descr =
   *   Structural_descr.t_sum [|
   *     [| F.packed_descr; F.packed_descr; Structural_descr.p_bool |];
   *     [| |]
   *   |]
   * let packed_descr = Structural_descr.pack structural_descr *)

  let compare x y =
    match x, y with
    | FRange.Itv (b1, e1, n1), FRange.Itv (b2, e2, n2) ->
      let c = Stdlib.compare n1 n2 in
      if c <> 0 then c else
        let r = F.compare b1 b2 in
        if r <> 0 then r else F.compare e1 e2
    | FRange.Itv _, FRange.NaN -> -1
    | FRange.NaN, FRange.Itv _ -> 1
    | FRange.NaN, FRange.NaN -> 0

  let equal x y =
    match x,y with
    | FRange.Itv (b1, e1, n1), FRange.Itv (b2, e2, n2) ->
      Cmp.equal b1 b2 && Cmp.equal e1 e2 && n1 = n2
    | FRange.NaN, FRange.NaN -> true
    | _ -> false

  let pretty_or_zero fmt f =
    if Cmp.(equal f pos_zero) then Format.fprintf fmt "0" else F.pretty fmt f

  let pretty fmt = function
    | FRange.Itv (b, e, nan) ->
      if Cmp.equal b e then
        Format.fprintf fmt "{%a%t}" pretty_or_zero b
          (fun fmt -> if nan then Format.pp_print_string fmt ";NaN")
      else
        let print_nan fmt nan =
          if nan then Format.fprintf fmt " %s {NaN}" ("U"(* Unicode.union_string ()*) )
        in
        Format.fprintf fmt "[%a .. %a]%a" F.pretty b F.pretty e print_nan nan
    | FRange.NaN -> Format.fprintf fmt "NaN"

  let hash = function
    | FRange.Itv (b, e, n) ->
      (2 * F.hash b) + (5 * F.hash e) + (7 * Hashtbl.hash n)
    | FRange.NaN -> 3

  (* ------------------------------------------------------------------------
                         Constants, injections, projections
     ------------------------------------------------------------------------ *)

  let min_and_max = function
    | FRange.Itv (b, e, nan) -> Some (b, e), nan
    | FRange.NaN -> None, true

  let nan = FRange.nan

  let inject = FRange.inject

  let singleton x = FRange.inject x x

  let pos_zero prec = singleton (Cst.pos_zero prec)
  let one prec = singleton (F.of_float Near prec 1.)
  let pos_infinity prec = singleton (Cst.pos_infinity prec)
  let _neg_infinity prec = singleton (Cst.neg_infinity prec)

  let minus_one_one prec ~nan =
    FRange.inject ~nan (F.of_float Near prec (-1.)) (F.of_float Near prec 1.)

  let m_pi = 3.1415929794311523 (* single-precision *)
  let m_minus_pi = -. m_pi

  let minus_pi_pi prec ~nan =
    FRange.inject ~nan (F.of_float Near prec m_minus_pi) (F.of_float Near prec m_pi)

  let top_f prec ~nan =
    FRange.inject ~nan (Cst.neg_infinity prec) (Cst.pos_infinity prec)

  let top prec = top_f prec ~nan:true

  let top_finite prec =
    if F.is_exact prec
    then
      let b = F.next_float prec (Cst.neg_infinity prec) in
      let e = F.prev_float prec (Cst.pos_infinity prec) in
      FRange.inject ~nan:false b e
    else top_f prec ~nan:false

  (* Splits the interval [b; e] into two intervals (neg, pos) such that [neg]
     (resp. [pos]) contains only negative (resp. positive) floats. *)
  let split_by_sign prec ((b, e) as x) =
    if Cmp.(le e neg_zero) then `Value x, `Bottom
    else if Cmp.(ge b pos_zero) then `Bottom, `Value x
    else `Value (b, Cst.neg_zero prec), `Value (Cst.pos_zero prec, e)

  (* Returns true iff [f] represents an odd integer. *)
  let is_odd prec f =
    let two = F.of_float Near prec 2. in
    Cmp.equal (F.abs (F.fmod Near prec f two)) Cmp.one

  (* [split_by_parity (b, e)] returns [even_bounds, odd_bounds], where:
     - [even_bounds] are the min and max even integers enclosed between b and e
       (or None if there is no such even integer).
     - [odd_bounds] are the min and max odd integers (representable as floating-
       point values) enclosed between b and e (or None if there is no such odd
       integer). These bounds cannot be infinite.  *)
  let split_by_parity prec (b, e) =
    let b = F.ceil b
    and e = F.floor e in
    if Cmp.gt_ieee b e then None, None
    else if Cmp.equal_ieee b e
    then
      let x = (b, e) in
      if is_odd prec b then None, Some x else Some x, None
    else
      let one = F.of_float Near prec 1. in
      let min_even, min_odd =
        if is_odd prec b
        (* No rounding errors may happen below because odd numbers are bounded. *)
        then F.add Near prec b one, b
        else b, Cmp.max (F.neg (Cst.max_odd_integer prec)) (F.add Near prec b one)
      and max_even, max_odd =
        if is_odd prec e
        (* No rounding errors may happen below because odd numbers are bounded. *)
        then F.sub Near prec e one, e
        else e, Cmp.min (Cst.max_odd_integer prec) (F.sub Near prec e one)
      in
      let even = Some (min_even, max_even)
      and odd =
        if Cmp.lt_ieee max_odd min_odd
        then None
        else Some (min_odd, max_odd)
      in
      even, odd

  (* -----------------------------------------------------------------------
                                      Lattice
     ----------------------------------------------------------------------- *)

  let is_included x1 x2 =
    match x1, x2 with
    | FRange.Itv (b1, e1, n1), FRange.Itv (b2, e2, n2) ->
      F.compare b2 b1 <= 0 && F.compare e1 e2 <= 0 && (not n1 || n2)
    | FRange.NaN, FRange.Itv (_, _, true) -> true
    | FRange.NaN, FRange.NaN -> true
    | _ -> false

  let join f1 f2 =
    match f1, f2 with
    | FRange.Itv (b1, e1, n1), FRange.Itv (b2, e2, n2) ->
      FRange.inject ~nan:(n1 || n2) (Cmp.min b1 b2) (Cmp.max e1 e2)
    | (FRange.Itv (b1, e1, _), FRange.NaN)
    | (FRange.NaN, FRange.Itv (b1, e1, _)) -> FRange.inject ~nan:true b1 e1
    | FRange.NaN, FRange.NaN -> FRange.nan

  let widen_down f = F.neg (F.widen_up (F.neg f))

  let widen f1 f2 =
    assert (is_included f1 f2);
    match f1, f2 with
    | FRange.Itv (b1, e1, _), FRange.Itv (b2, e2, nan) ->
      let b = if Cmp.equal b2 b1 then b2 else widen_down b2 in
      let e = if Cmp.equal e2 e1 then e2 else F.widen_up e2 in
      (** widen_up and down produce double only if the input is a double *)
      FRange.inject ~nan b e
    | FRange.NaN, f2 -> f2
    | FRange.Itv _, FRange.NaN -> assert false

  let meet f1 f2 =
    match f1, f2 with
    | FRange.Itv (b1, e1, n1), FRange.Itv (b2, e2, n2) ->
      let is_finite = F.compare b2 e1 <= 0 && F.compare b1 e2 <= 0 in
      let is_nan = n1 && n2 in
      if is_finite || is_nan
      then
        let v =
          if is_finite then
            FRange.inject ~nan:is_nan (Cmp.max b1 b2) (Cmp.min e1 e2)
          else FRange.nan
        in
        `Value v
      else `Bottom
    | (FRange.Itv (_, _, true) | FRange.NaN),
      (FRange.Itv (_, _, true) | FRange.NaN) -> `Value FRange.nan
    | _ -> `Bottom

  let narrow = meet

  (* ------------------------------------------------------------------------
                                   Tests
     ------------------------------------------------------------------------ *)

  let contains cst = function
    | FRange.NaN -> false
    | FRange.Itv (b, e, _) -> Cmp.le b cst && Cmp.ge e cst

  let contains_pos_zero t = contains Cmp.pos_zero t
  let contains_neg_zero t = contains Cmp.neg_zero t

  let contains_a_zero = function
    | FRange.Itv (b, e, _) ->
      Cmp.(le_ieee b pos_zero) && Cmp.(ge_ieee e pos_zero)
    | FRange.NaN -> false

  let contains_non_zero = function
    | FRange.Itv (b, e, nan) ->
      nan || Cmp.(lt_ieee b pos_zero) || Cmp.(gt_ieee e pos_zero)
    | FRange.NaN -> true

  let contains_strictly_pos = function
    | FRange.Itv (_, e, _) -> Cmp.(gt_ieee e pos_zero)
    | FRange.NaN -> false

  let contains_strictly_neg = function
    | FRange.Itv (b, _, _) -> Cmp.(lt_ieee b neg_zero)
    | FRange.NaN -> false

  let contains_strict_neg_finite (b, e) =
    Cmp.(gt_ieee e neg_infinity) && Cmp.(lt b neg_zero)

  let contains_finite_noninteger prec (b, e) =
    let ib = F.ceil b in
    not ((Cmp.equal_ieee ib b && Cmp.equal_ieee b e)
         || Cmp.le_ieee e (F.neg (Cst.ceil_of_max_non_integer prec))
         || Cmp.ge_ieee b (Cst.ceil_of_max_non_integer prec))

  let contains_pos_infinity = function
    | FRange.Itv (_, e, _) -> Cmp.is_pos_infinity e
    | FRange.NaN  -> false

  let contains_neg_infinity = function
    | FRange.Itv (b, _, _) -> Cmp.is_neg_infinity b
    | FRange.NaN  -> false

  let contains_infinity f =
    contains_pos_infinity f || contains_neg_infinity f

  let contains_nan = function
    | FRange.NaN -> true
    | FRange.Itv (_, _, nan) -> nan

  let is_singleton = function
    | FRange.NaN -> true
    | FRange.Itv (b, e, nan) -> Cmp.equal b e && not nan

  let is_one = function
    | FRange.NaN -> false
    | FRange.Itv (b, e, nan) ->
      Cmp.(equal b one) && Cmp.(equal e one) && not nan

  let is_a_zero = function
    | FRange.NaN -> false
    | FRange.Itv (b, e, nan) ->
      not nan && Cmp.is_a_zero b && Cmp.is_a_zero e

  let if_not_nan = function
    | FRange.NaN -> assert false
    | FRange.Itv (b, e, _) -> b, e

  let is_not_nan = function
    | FRange.NaN -> Comp.False
    | FRange.Itv (_b, _e, nan) -> if nan then Comp.Unknown else Comp.True

  let is_finite = function
    | FRange.NaN -> Comp.False
    | FRange.Itv (b, e, nan) ->
      if Cmp.is_neg_infinity e || Cmp.is_pos_infinity b
      then Comp.False
      else if nan || Cmp.is_neg_infinity b || Cmp.is_pos_infinity e
      then Comp.Unknown
      else Comp.True

  let is_negative = function
    | FRange.Itv (b, e, false) ->
      if F.is_negative e then Comp.True
      else if not (F.is_negative b) then Comp.False
      else Comp.Unknown
    | FRange.Itv (_, _, true)
    | FRange.NaN -> Comp.Unknown

  let backward_is_not_nan = function
    | FRange.NaN -> `Bottom
    | FRange.Itv (b, e, _) -> `Value (FRange.inject ~nan:false b e)

  let backward_is_finite prec = function
    | FRange.NaN -> `Bottom
    | FRange.Itv (b, e, _) as f ->
      if Cmp.equal b e && F.is_infinite b
      then `Bottom (* [f] is exactly an infinite, we can return `Bottom even
                      in the [Real] case *)
      else narrow (top_finite prec) f

  let has_greater_min_bound t1 t2 =
    match t1, t2 with
    | FRange.Itv (m1, _, _), FRange.Itv (m2, _, _) -> F.compare m1 m2
    | FRange.NaN, FRange.Itv _ -> 1
    | FRange.Itv _, FRange.NaN -> -1
    | FRange.NaN, FRange.NaN -> 0

  let has_smaller_max_bound t1 t2 =
    match t1, t2 with
    | FRange.Itv (_, m1, _), FRange.Itv (_, m2, _) -> F.compare m2 m1
    | FRange.NaN, FRange.Itv _ -> 1
    | FRange.Itv _, FRange.NaN -> -1
    | FRange.NaN, FRange.NaN -> 0

  (* ------------------------------------------------------------------------
                                Comparisons
     ------------------------------------------------------------------------ *)

  let forward_eq (b1,e1) (b2,e2) =
    let not_intersects = Cmp.lt_ieee e2 b1 || Cmp.lt_ieee e1 b2 in
    if not_intersects
    then Comp.False
    else if Cmp.equal_ieee b1 e1 && Cmp.equal_ieee b2 e2
    then Comp.True
    else Comp.Unknown

  let forward_le (b1, e1) (b2, e2) =
    if Cmp.le_ieee e1 b2 then Comp.True
    else if Cmp.lt_ieee e2 b1 then Comp.False
    else Comp.Unknown

  let forward_lt (b1, e1) (b2, e2) =
    if Cmp.lt_ieee e1 b2 then Comp.True
    else if Cmp.le_ieee e2 b1 then Comp.False
    else Comp.Unknown

  let forward_comp op f1 f2 = match f1, f2 with
    | FRange.NaN, _ | _, FRange.NaN ->
      if op = Comp.Ne then Comp.True else Comp.False
    | FRange.Itv (b1, e1, nan1), FRange.Itv (b2, e2, nan2) ->
      let r = match op with
        | Comp.Le -> forward_le (b1, e1) (b2, e2)
        | Comp.Ge -> forward_le (b2, e2) (b1, e1)
        | Comp.Lt -> forward_lt (b1, e1) (b2, e2)
        | Comp.Gt -> forward_lt (b2, e2) (b1, e1)
        | Comp.Eq -> forward_eq (b1, e1) (b2, e2)
        | Comp.Ne -> Abstract_interp.inv_truth (forward_eq (b1, e1) (b2, e2))
      in
      if nan1 || nan2
      then
        if op = Comp.Ne
        then (match r with Comp.True -> Comp.True | _ -> Comp.Unknown)
        else (match r with Comp.False -> Comp.False | _ -> Comp.Unknown)
      else r

  (* This function intentionally returns different results with
     [e2 = -0.] and [e2 = 0.] *)
  let backward_le_aux prec (b1, e1) e2 =
    if not (Cmp.le b1 e2)
    then `Bottom
    else if Cmp.le e1 e2
    then `Value (FRange.inject b1 e1)
    else FRange.inject_after_tighten prec ~nan:false b1 e2

  (* This is the "real" backward transformer for [le], which does not distinguish
     [0.] and [-0.]. Thus we enlarge the bound in the "worst" direction. *)
  let backward_le prec (b1, e1) e2 =
    let e2 = if Cmp.is_a_zero e2 then Cst.pos_zero prec else e2 in
    backward_le_aux prec (b1, e1) e2

  let backward_lt prec ((b1, e1) as f1) e2 =
    if Cmp.le_ieee e2 b1
    then `Bottom
    else if F.is_exact prec || Cmp.equal b1 e1
    then backward_le prec f1 (Cmp.prev_float_ieee prec e2)
    else
      (* On real we cannot be more precise than [le], except on zeros: at
         least get rid of the "bad" zero *)
      let e2 = if Cmp.is_a_zero e2 then Cst.neg_zero prec else e2 in
      backward_le_aux prec f1 e2

  (* see comments in {!backward_le_aux} *)
  let backward_ge_aux prec (b1, e1) b2 =
    if not (Cmp.le b2 e1)
    then `Bottom
    else if Cmp.le b2 b1
    then `Value (FRange.inject b1 e1)
    else FRange.inject_after_tighten prec ~nan:false b2 e1

  (* see comments in {!backward_le} *)
  let backward_ge prec (b1, e1) b2 =
    let b2 = if Cmp.is_a_zero b2 then Cst.neg_zero prec else b2 in
    backward_ge_aux prec (b1, e1) b2

  (* see comments in {!backward_gt} *)
  let backward_gt prec ((b1, e1) as f1) b2 =
    if Cmp.le_ieee e1 b2
    then `Bottom
    else if F.is_exact prec || Cmp.equal b1 e1
    then backward_ge prec f1 (Cmp.next_float_ieee prec b2)
    else
      let b2 = if Cmp.is_a_zero b2 then Cst.pos_zero prec else b2 in
      backward_ge_aux prec f1 b2

  (** The operands cannot be {!Nan} *)
  let backward_comp_left_true_finite op prec f1' f2' =
    let f1 = if_not_nan f1' in
    let (b2, e2) = if_not_nan f2' in
    match op with
    | Comp.Le -> backward_le prec f1 e2
    | Comp.Ge -> backward_ge prec f1 b2
    | Comp.Lt -> backward_lt prec f1 e2
    | Comp.Gt -> backward_gt prec f1 b2
    | Comp.Eq ->
      (* -0 and +0 must not be distinguished here *)
      let f2 =
        if contains_a_zero f2'
        then join f2' (FRange.inject (Cst.neg_zero prec) (Cst.pos_zero prec))
        else f2'
      in
      narrow f1' f2
    | Comp.Ne ->
      (* compute (f1 ∩ [-infty,min[ ) ∪ (f1 ∩ ]max,infty]) *)
      let before_or_after min max =
        Bottom.join join
          (backward_lt prec f1 min) (backward_gt prec f1 max)
      in
      (* As usual, we cannot reduce if [f2] is not a singleton, except that
         the two zeros are a kind of singleton. Checking whether [f2] is on
         a frontier of [f1] is not obvious because of the multiple cases
         (and [allmodes]) so we use the transformers for [lt] instead. *)
      if is_a_zero f2' then before_or_after (Cst.neg_zero prec) (Cst.pos_zero prec)
      else if is_singleton f2' then before_or_after b2 b2
      else `Value f1'

  (* Applies [backward f1 f2] and removes NaN from [f1] and [f2]. *)
  let backward_comp_no_nan backward_finite f1 f2 =
    match f1, f2 with
    | FRange.NaN, _ | _, FRange.NaN -> `Bottom
    | FRange.Itv (b, e, nan), FRange.Itv _ ->
      let f1 = if nan then FRange.inject ~nan:false b e else f1 in
      backward_finite f1 f2

  (* Applies [backward f1 f2] but preserves NaN from [f1] and [f2]. *)
  let backward_comp_with_nan backward_finite f1 f2 =
    if contains_nan f2
    then `Value f1
    else
      match f1 with
      | FRange.NaN -> `Value f1
      | FRange.Itv (_, _, nan) ->
        let nan = if nan then `Value FRange.nan else `Bottom in
        Bottom.join join (backward_finite f1 f2) nan

  let backward_comp_left_true op prec =
    let backward_finite = backward_comp_left_true_finite op prec in
    if op = Comp.Ne
    then backward_comp_with_nan backward_finite
    else backward_comp_no_nan backward_finite

  let backward_comp_left_false op prec =
    let backward_finite = backward_comp_left_true_finite (Comp.inv op) prec in
    if op = Comp.Ne
    then backward_comp_no_nan backward_finite
    else backward_comp_with_nan backward_finite


  (* ------------------------------------------------------------------------
                        Simple arithmetic operations
     ------------------------------------------------------------------------ *)

  (* The functions defined using [exact_aux] below are, among other
     properties, (1) exact (the result as a real can always be
     represented exactly, in the good type), and (2) total. In
     particular, given a float 'x', 'ff x == (float)(f (double)x)'.
     Thus, in this module, the 'f' functions are also the non-f
     (since float32 are represented using double) *)
  let (>>) t f = match t with
    | FRange.NaN -> t
    | FRange.Itv (b, e, nan) -> FRange.inject ~nan (f b) (f e)

  let floor t = t >> F.floor
  let ceil t = t >> F.ceil
  let trunc t = t >> F.trunc
  let fround t = t >> F.fround

  let neg = function
    | FRange.Itv (b, e, nan) ->
      (* do not round because exact operation *)
      FRange.inject ~nan (F.neg e) (F.neg b)
    | FRange.NaN -> FRange.nan

  let abs prec = function
    | FRange.Itv (b, e, nan) as f ->
      if contains Cmp.pos_zero f then
        let zero = Cst.pos_zero prec in
        FRange.inject ~nan zero (Cmp.max (F.abs b) (F.abs e))
      else (* f is either strictly positive or strictly negative *)
      if F.compare e Cmp.pos_zero < 0 then neg f else f
    | FRange.NaN as f -> f

  (* This monad returns a NaN if one operand can only be NaN, and lets the
     second function perform the computation if both operands contain a
     non-empty floating-point interval. *)
  let ( >>% ) = fun (x,y) f -> match x, y with
    | FRange.NaN, _ | _, FRange.NaN -> FRange.nan
    | FRange.Itv (b1, e1, nan1), FRange.Itv (b2, e2, nan2) ->
      let nan = nan1 || nan2 in
      f ~nan (b1, e1) (b2, e2)

  (* Auxiliary function used for the forward semantics of add, mul and div.
     For a monotonic function [op], the bounds of [[b1..e1] op [b2..e2]] are the
     minimum and maximum of [b1 op b2], [b1 op e2], [e1 op b2] and [e1 op e2].
     NaN can be created from \infty - \infty, 0 * \infty, 0/0 and \infty /
     \infty, in which case the result contains NaN, and new operations are
     performed to take into account the results of values near \infty and 0.
     Beware that NaN and discontinuities occuring between the bounds of the
     arguments (i.e. on zeros, as an infinity is always a bound) should be
     checked and processed by the caller. *)
  let monotonic op prec x y =
    (x, y) >>% fun ~nan (b1, e1) (b2, e2) ->
    let nan = ref nan in
    (* Results of [op] applied to the bounds of the intervals, excluding NaN. *)
    let results = ref [] in
    (* When [a op b = NaN], performs new operations to take into account values
       near [a] (and near [b] with the same reasoning). For such a NaN from add,
       mul or div, [c op b] is constant for all values c <> a of the same sign.
       Thus, we can replace [a] by any of these values.
       - if [x] is a singleton or [-0 .. +0], there are no values other than the
         bounds to take into account;
       - otherwise, if [a] is infty, replace it by 1 with the sign of [a]
         (no risk of NaN with 1 on add, mul and div);
       - otherwise, if [a] is zero, the other bound [c] of the interval has the
         same sign as the values near [a] in the interval; as [c op b] is also
         computed, no need to perform a new operation. *)
    let treat_nan_result rnd a b =
      nan := true;
      if F.is_infinite a && not (Cmp.equal b1 e1)
      then results := op rnd prec (Cst.sign prec a) b :: !results;
      if F.is_infinite b && not (Cmp.equal b2 e2)
      then results := op rnd prec a (Cst.sign prec b) :: !results;
    in
    let op rnd x y =
      let r = op rnd prec x y in
      if F.is_nan r
      then treat_nan_result rnd x y
      else results := r :: !results
    in
    let compute rnd =
      results := [];
      op rnd e1 e2; op rnd e1 b2; op rnd b1 e2; op rnd b1 b2;
    in
    let rounding_mode = if F.is_exact prec then Near else Down in
    compute rounding_mode;
    let pos_inf = Cst.pos_infinity prec in
    let min = List.fold_left Cmp.min pos_inf !results in
    if not (F.is_exact prec)
    then compute Up;
    let neg_inf = Cst.neg_infinity prec in
    let max = List.fold_left Cmp.max neg_inf !results in
    if min > max
    then (assert !nan; FRange.nan)
    else FRange.inject ~nan:!nan min max

  let add prec = monotonic F.add prec
  let sub prec = monotonic F.sub prec

  let mul prec x y =
    let r = monotonic F.mul prec x y in
    (* A NaN may occur between the bounds of the intervals, on 0 * \infty. *)
    if (contains_infinity x && contains_a_zero y) ||
       (contains_infinity y && contains_a_zero x)
    then FRange.add_nan r
    else r

  let div prec x y =
    let r = monotonic F.div prec x y in
    (* A NaN may occur between the bounds of the intervals, on 0/0. *)
    let nan = (contains_a_zero x && contains_a_zero y) in
    (* Treat the discontinuity around 0: divisions by 0 produce infinites. *)
    let pos_inf =
      contains Cmp.pos_zero y && contains_strictly_pos x ||
      contains Cmp.neg_zero y && contains_strictly_neg x
    and neg_inf =
      contains Cmp.pos_zero y && contains_strictly_neg x ||
      contains Cmp.neg_zero y && contains_strictly_pos x
    in
    let r = if pos_inf then FRange.add_pos_infinity prec r else r in
    let r = if neg_inf then FRange.add_neg_infinity prec r else r in
    if nan then FRange.add_nan r else r

  (* Could be improved a lot, cf [Marre10]. *)
  let backward_add_one prec ~other ~result =
    (* No reduction when the result contains an infinity, and when the result and
       the other operand contain NaN (as x + NaN = NaN for any x). *)
    if contains_infinity result || (contains_nan other && contains_nan result)
    then `Value (top prec)
    else
      (* Values that can lead to NaN in the result. *)
      let reduce_for_nan t =
        let t =
          if contains_pos_infinity other
          then FRange.add_neg_infinity prec t
          else t
        in
        if contains_neg_infinity other
        then FRange.add_pos_infinity prec t
        else t
      in
      let reduced_for_nan =
        if contains_nan result then `Value (reduce_for_nan FRange.nan) else `Bottom
      in
      (* Values that can lead to finite values in the result. *)
      let reduced_for_finite_values =
        match result, other with
        | FRange.NaN, _ | _, FRange.NaN  -> `Bottom
        | FRange.Itv (bres, eres, _), FRange.Itv (bother, eother, _) ->
          let bres = Cmp.prev_float_ieee prec bres in
          let eres = Cmp.next_float_ieee prec eres in
          let round = if F.is_exact prec then Up else Near in
          let b = F.sub round prec bres eother in
          let round = if F.is_exact prec then Down else Near in
          let e = F.sub round prec eres bother in
          if Cmp.le b e then `Value (FRange.inject ~nan:false b e) else `Bottom
      in
      Bottom.join join reduced_for_finite_values reduced_for_nan

  let backward_add fkind ~left ~right ~result =
    backward_add_one fkind ~other:right ~result >>- fun left' ->
    backward_add_one fkind ~other:left ~result >>- fun right' ->
    `Value (left', right')

  let backward_sub fk ~left ~right ~result =
    let right = neg right in
    backward_add fk ~left ~right ~result
    >>-: fun (left, right) -> (left, neg right)

  (* ------------------------------------------------------------------------
                             Exp Log Sqrt Pow Fmod
     ------------------------------------------------------------------------ *)

  let (>>:) t f = match t with
    | FRange.NaN -> t
    | FRange.Itv (b, e, nan) -> f ~nan b e

  let approx f prec ~nan b e =
    let round = if F.is_exact prec then Near else Down in
    let min = f round prec b in
    let round = if F.is_exact prec then Near else Up in
    let max = f round prec e in
    FRange.inject ~nan min max

  let exp prec t = t >>: approx F.exp prec

  let log_aux log prec t =
    t >>: fun ~nan b e ->
    if Cmp.(lt e neg_zero)
    then FRange.nan
    else
      let nan = nan || Cmp.(lt b neg_zero) in
      let b = Cmp.max (Cst.neg_zero prec) b in
      approx log prec ~nan b e

  let log prec = log_aux F.log prec
  let log10 prec = log_aux F.log10 prec

  (* [sqrt_f] is the actual function computing the (exact) square root,
     in single precision (sqrtf) or double precision (sqrt). *)
  let sqrt prec t =
    t >>: fun ~nan b e ->
    if Cmp.(lt_ieee e neg_zero)
    then FRange.nan
    else
      let nan, b =
        if Cmp.(ge_ieee b neg_zero)
        then nan, b
        else true, Cst.neg_zero prec
      in
      approx F.sqrt prec ~nan b e


  let value_if prec condition f = if condition then `Value (f prec) else `Bottom

  (* Returns the minimal or maximal (according to [min_or_max]) results of the
     binary operation [op] applied to the bounds of the intervals [b1..e1] and
     [b2..e2]. *)
  let extremum min_or_max op (b1, e1) (b2, e2) =
    let extremum4 a b c d = min_or_max a (min_or_max b (min_or_max c d)) in
    extremum4 (op b1 b2) (op b1 e2) (op e1 b2) (op e1 e2)

  (* Returns the minimum and maximum results of the binary operation [op] applied
     to the bounds of the intervals [b1..e1] and [b2..e2]. *)
  let extrema op (b1, e1) (b2, e2) =
    let a = op b1 b2 and b = op b1 e2 and c = op e1 b2 and d = op e1 e2 in
    Cmp.min a (Cmp.min b (Cmp.min c d)), Cmp.max a (Cmp.max b (Cmp.max c d))

  (* Computes [pow_f] on a negative [bx; ex] interval (including infinites).
     Processes by disjunction over even and odd integers enclosed within [by; ey].
     [pow] is then monotonic on even integers (including zeros and infinities),
     and on odd integers (except on infinities). *)
  let pow_negative_x pow_f prec (bx, ex as x) (by, ey as y) =
    let even, odd = split_by_parity prec y in
    (* Even integers [y] lead to positive results, while odd ones lead to negative
       results.  When [y] contains both even and odd integers, the minimum result
       is in odd integers, and the maximum in even integers. *)
    let min, max = match even, odd with
      | None, None -> Cst.pos_infinity prec, Cst.neg_infinity prec
      | Some even, None -> extrema pow_f x even
      | None, Some odd  -> extrema pow_f x odd
      | Some even, Some odd ->
        extremum Cmp.min pow_f x odd, extremum Cmp.max pow_f x even
    in
    let nonint_y = contains_finite_noninteger prec y in
    (* pow creates NaN when [x] is a negative non-zero finite value, and [y] a
       non integer value. *)
    let nan = contains_strict_neg_finite x && nonint_y in
    (* Special cases of neg_infinity and neg_zero for [x], that do not produce
       a NaN on non integer [y], unlike strictly negative finite values [x]. *)
    let neg_nonint_y = nonint_y && Cmp.(lt by neg_zero) in
    let pos_nonint_y = nonint_y && Cmp.(gt ey pos_zero) in
    let neg_infinity_x = Cmp.is_neg_infinity bx in
    let zero_x = Cmp.(equal ex neg_zero) in
    Bottom.join_list join
      [ if Cmp.le min max then `Value (FRange.inject min max) else `Bottom;
        value_if prec nan (fun _ -> FRange.nan);
        value_if prec (neg_infinity_x && neg_nonint_y) pos_zero;
        value_if prec (neg_infinity_x && pos_nonint_y) pos_infinity;
        value_if prec (zero_x && neg_nonint_y) pos_infinity;
        value_if prec (zero_x && pos_nonint_y) pos_zero ]

  (* Computes pow on a positive [bx; ex] interval (including infinites): the
     function is continuous and monotonic. *)
  let pow_positive_x pow_f _prec x y =
    let min, max = extrema pow_f x y in
    FRange.inject min max

  let pow prec x y =
    if is_one x || is_a_zero y then one prec
    else
      (x, y) >>% fun ~nan itv_x itv_y ->
      let pow_f = F.pow Near prec in
      (* Split the x interval around zeros, as pow is discontinuous on zeros. *)
      let neg_x, pos_x = split_by_sign prec itv_x in
      let pos_x_res = pos_x >>-: fun x -> pow_positive_x pow_f prec x itv_y in
      let neg_x_res = neg_x >>- fun x -> pow_negative_x pow_f prec x itv_y in
      let nan_res = if nan then `Value FRange.nan else `Bottom in
      Bottom.non_bottom
        (Bottom.join_list join [ pos_x_res; neg_x_res; nan_res ])

  (* Is [fmod] continuous on positive intervals [b1..e1] and [b2..e2]?
     This is the case if for all x, y in these intervals, the rounded quotient
     [floor(x/y)] is constant, as [fmod x y = x - floor(x/y) * y] when [x] and [y]
     are positive.
     Also checks that [x/y < 2^53], otherwise truncation to an integer may return
     an incorrect result. Note: to avoid issues with rounding, we conservatively
     set the limit to 2^51 instead of 2^53 (and to 2^21 instead of 2^23 in single
     precision). *)
  let is_fmod_continuous prec (b1, e1) (b2, e2) =
    (* Discontinuity of [fmod x y] on infinite [x] and on zero [y]. *)
    F.is_finite e1 && Cmp.(gt_ieee b2 pos_zero) &&
    let four = F.of_float Near prec 4. in
    let max_i = F.div Near prec (Cst.max_precise_integer prec) four in
    let f1 = F.floor (F.div Zero prec b1 e2) in
    let f2 = F.floor (F.div Zero prec e1 b2) in
    Cmp.equal f1 f2 && Cmp.le f1 max_i && F.is_exact prec

  (* Forward semantics of fmod on positive intervals. [x] must contain finite
     values, and [y] must contain non-zero values, in which case finite values
     are produced. This function does not check the creation of NaN. *)
  let positive_fmod prec (b1, e1 as x) (b2, e2 as y) =
    let fmod = F.fmod Near prec in
    (* Singleton case. [x] cannot be infinite, and [y] cannot be zero. *)
    if Cmp.equal b1 e1 && Cmp.equal b2 e2
    then let c = fmod b1 b2 in c, c
    (* If all values of x are smaller than all values of [y], [x] is unchanged. *)
    else if Cmp.lt_ieee e1 b2
    then x
    (* If fmod is continuous on the intervals [x] and [y], it is also monotonic,
       and the bounds of the result are the remainders of the bounds. *)
    else if is_fmod_continuous prec x y
    then fmod b1 e2, fmod e1 b2
    (* Otherwise, fmod always satisfies 0 <= [fmod x y] <= x and [fmod x y] < y. *)
    else let max = Cmp.min e1 (F.prev_float prec e2) in (Cst.pos_zero prec, max)

  let fmod prec x y =
    (* [fmod x (-y)] = [fmod x y], so use only positive [y]. *)
    let y = abs prec y in
    (x, y) >>% fun ~nan (b1, e1) (b2, e2) ->
    (* If [x] is an infinite singleton, or if [y] contains only zero, only NaN
       can be created. *)
    if (Cmp.equal b1 e1 && F.is_infinite b1)
    || (Cmp.equal b2 e2 && Cmp.is_a_zero b2)
    then FRange.nan
    else
      let nan = nan || contains_infinity x || contains_a_zero y in
      let positive_fmod x =
        let b, e = positive_fmod prec x (b2, e2) in
        FRange.inject ~nan b e
      in
      (* Process by disjunction on the sign of [x], and join the results for
         negative x and for positive x. *)
      let neg_x, pos_x = split_by_sign prec (b1, e1) in
      let neg_itv (b, e) = F.neg e, F.neg b in
      let neg_r = neg_x >>-: neg_itv >>-: positive_fmod >>-: neg in
      let pos_r = pos_x >>-: positive_fmod in
      Bottom.non_bottom (Bottom.join join neg_r pos_r)

  (* ------------------------------------------------------------------------
                                Trigonometry
     ------------------------------------------------------------------------ *)

  (* It is wrong to use m_pi as the local minimum as was previously done,
     because:
     A = 3.14159265358979312 and A < m_pi, cos A = 1. and cos m_pi = 0.999
     Moreover it is not due to the imprecision in the value of pi:
     A < pred m_pi in simple.

     So we use a quarter of the interval [-pi:pi] to be safe. (But still,
     nothing proves that cos or sin are monotonic on those ranges.) *)
  let m_pi = 3.1415929794311523 (* single-precision *)
  let cos_sin_security = F.of_float Near Float_sig.Single (m_pi /. 4.)
  let cos_sin op prec t =
    t >>: fun ~nan b e ->
    (* Special case when at least a bound is infinite. *)
    if F.is_infinite b || F.is_infinite e
    then if Cmp.equal b e then FRange.nan else minus_one_one prec ~nan:true
    (* [b] and [e] are finite. Precise case for a singleton. *)
    else if Cmp.equal b e
    then let c = op Near prec b in FRange.inject ~nan c c
    else if Cmp.le_ieee b (F.neg cos_sin_security)
         || Cmp.ge_ieee e cos_sin_security
    then minus_one_one prec ~nan
    else
      let xl =
        if Cmp.(lt b pos_zero) && Cmp.(lt pos_zero e)
        then [b; e; Cst.pos_zero prec]
        else [b; e]
      in
      let l = List.map (op Near prec) xl in
      let min_f = List.fold_left Cmp.min (Cst.pos_infinity prec) l in
      let max_f = List.fold_left Cmp.max (Cst.neg_infinity prec) l in
      FRange.inject ~nan min_f max_f

  let cos prec = cos_sin F.cos prec
  let sin prec = cos_sin F.sin prec

  let atan2 prec x y =
    (x, y) >>% fun ~nan (b1, e1) (b2, e2) ->
    let op = F.atan2 Near prec in
    if Cmp.equal b1 e1 && Cmp.equal b2 e2
    then let c = op b1 b2 in FRange.inject ~nan c c
    else
      (* Unless y ([b1,e1]) crosses the x-axis, atan2 is continuous,
         and its minimum/maximum are at the ends of the intervals of x and y.
         Otherwise, the result is [-pi,+pi]. *)
    if Cmp.(lt b1 pos_zero) && Cmp.(gt e1 neg_zero)
    then minus_pi_pi prec ~nan
    else
      let a1, a2, a3, a4 = op b1 b2, op b1 e2, op e1 b2, op e1 e2 in
      let b = Cmp.min a1 (Cmp.min a2 (Cmp.min a3 a4)) in
      let e = Cmp.max a1 (Cmp.max a2 (Cmp.max a3 a4)) in
      match prec with
      | Float_sig.Single ->
        (* Rounding of atan2f in single-precision may go against monotony and
           reach the next (previous) float after (before) the bounds. *)
        FRange.inject ~nan
          (Cmp.prev_float_ieee prec b)
          (Cmp.next_float_ieee prec e)
      | _ -> FRange.inject ~nan b e

  (* ------------------------------------------------------------------------
                                   Casts
     ------------------------------------------------------------------------ *)

  let forward_cast ~dst = function
    | FRange.NaN -> nan
    | FRange.Itv (b, e, nan) ->
      let round = F.round_to_precision Near dst in
      inject ~nan (round b) (round e)

  (* This function must make sure to return a result with float 32 bounds *)
  let backward_cast ~src = function
    | FRange.NaN -> `Value nan
    | FRange.Itv (b, e, nan) -> FRange.inject_after_tighten src ~nan b e

  let cast_int_to_float prec min max =
    let min = match min with
      | None -> Cst.neg_infinity prec
      | Some v ->
        let round = if F.is_exact prec then Near else Down in
        F.of_float round prec (Integer.to_float v)
    in
    let max = match max with
      | None -> Cst.neg_infinity prec
      | Some v ->
        let round = if F.is_exact prec then Near else Up in
        F.of_float round prec (Integer.to_float v)
    in
    FRange.inject min max

  (* Bitwise reinterpretation of a floating-point value into consecutive
     ranges of integer. This operation is exact in terms of
     concretization. 'Parametric' in the number of bits. *)
  let bits_of_float_list ~prec ~bits_of_float ~max_int =
    let neg_infinity = Cst.neg_infinity prec
    and pos_infinity = Cst.pos_infinity prec in
    let itvs_nan =
      let smallest_neg_nan = Integer.succ (bits_of_float neg_infinity) in
      let biggest_neg_nan  = Integer.minus_one in
      let smallest_pos_nan = Integer.succ (bits_of_float pos_infinity) in
      let biggest_pos_nan  = max_int in
      [(smallest_neg_nan, biggest_neg_nan);
       (smallest_pos_nan, biggest_pos_nan)]
    in
    function
    | FRange.NaN -> itvs_nan
    | FRange.Itv (b, e, nan) ->
      let nan = if nan then itvs_nan else [] in
      let neg, pos = split_by_sign prec (b, e) in
      let neg = neg >>-: fun (b, e) -> bits_of_float e, bits_of_float b in
      let pos = pos >>-: fun (b, e) -> bits_of_float b, bits_of_float e in
      Bottom.add_to_list pos (Bottom.add_to_list neg nan)

  let bits_of_float64_list =
    let bits_of_float f =
      Integer.of_int64 (Int64.bits_of_float (F.to_float f))
    in
    let max_int = Integer.of_int64 Int64.max_int in
    bits_of_float_list ~prec:Double ~bits_of_float ~max_int

  let bits_of_float32_list =
    let bits_of_float f =
      Integer.of_int32 (Int32.bits_of_float (F.to_float f))
    in
    let max_int = Integer.of_int32 Int32.max_int in
    bits_of_float_list ~prec:Single ~bits_of_float ~max_int

  (* ------------------------------------------------------------------------
                                   Subdivision
     ------------------------------------------------------------------------ *)

  (* [avg] and [split] implement two different strategies for cutting a
     floating-point interval in half: [avg] computes the mathematical average of
     the two bounds, while [split] balances the number of representable values
     of the given precision in each resulting intervals. *)

  (* Computes the average between two ocaml doubles. *)
  let avg x y =
    let fx = F.to_float x and fy = F.to_float y in
    if F.is_negative x && F.is_negative y
    then fy +. (fx -. fy) /. 2.
    else (fx +. fy) /. 2.

  (* assumption: [0. <= x <= y]. returns the median of the range [x..y]
     in number of double values. *)
  let split_positive prec x y =
    let ix = Int64.bits_of_float (F.to_float x) in
    let iy = Int64.bits_of_float (F.to_float y) in
    let f = Int64.(float_of_bits (add ix (div (sub iy ix) 2L))) in
    F.of_float Near prec f

  (* assumption: [x <= y] *)
  let _split prec x y =
    match F.is_negative x, F.is_negative y with
    | false, false -> split_positive prec x y
    | true, true -> F.neg (split_positive prec (F.neg x) (F.neg y))
    | true, false -> Cst.neg_zero prec
    | false, true -> assert false

  exception Can_not_subdiv = Abstract_interp.Can_not_subdiv

  let subdivide prec t =
    assert (prec = Single || prec = Double); (* See Value/Value#105 *)
    match t with
    | FRange.NaN -> raise Can_not_subdiv
    | FRange.Itv (b, e, nan) ->
      if nan
      then FRange.inject ~nan:false b e, FRange.nan
      else if Cmp.equal b e
      then raise Can_not_subdiv
      else
        let midpoint, smidpoint =
          (* Infinities are interesting points to consider specially. *)
          if Cmp.is_neg_infinity b
          then Cst.neg_infinity prec, F.next_float prec (Cst.neg_infinity prec)
          else if Cmp.is_pos_infinity e
          then F.prev_float prec (Cst.pos_infinity prec), Cst.pos_infinity prec
          else if Cmp.equal (F.next_float prec b) e
          then b, e
          else
            let midpoint = avg b e in
            let midpoint = F.of_float Down prec midpoint in
            let smidpoint =
              if F.is_exact prec then F.next_float prec midpoint else midpoint
            in
            if Cmp.le smidpoint e
            then midpoint, smidpoint
            else midpoint, e
        in
        if Cmp.le smidpoint b || Cmp.le e midpoint then raise Can_not_subdiv;
        let i1 = FRange.inject ~nan b midpoint in
        assert (is_included i1 t);
        let i2 = FRange.inject ~nan smidpoint e in
        assert (is_included i2 t);
        i1, i2

end