package codex

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

Source file floating_point.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
(**************************************************************************)
(*  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).                      *)
(*                                                                        *)
(**************************************************************************)

type c_rounding_mode =
    FE_ToNearest | FE_Upward | FE_Downward | FE_TowardZero

let string_of_c_rounding_mode = function
  | FE_ToNearest -> "FE_NEAREST"
  | FE_Upward -> "FE_UPWARD"
  | FE_Downward -> "FE_DOWNWARD"
  | FE_TowardZero -> "FE_TOWARDZERO"

(* replace "noalloc" with [@@noalloc] for OCaml version >= 4.03.0 *)
[@@@ warning "-3"]

external set_round_downward: unit -> unit = "fc_ival_set_round_downward" "noalloc"
external set_round_upward: unit -> unit = "fc_ival_set_round_upward" "noalloc"
external set_round_nearest_even: unit -> unit = "fc_ival_set_round_nearest_even" "noalloc"
external set_round_toward_zero : unit -> unit = "fc_ival_set_round_toward_zero" "noalloc"
external get_rounding_mode: unit -> c_rounding_mode = "fc_ival_get_rounding_mode" "noalloc"
external set_rounding_mode: c_rounding_mode -> unit = "fc_ival_set_rounding_mode" "noalloc"

[@@@ warning "+3"]

external round_to_single_precision_float: float -> float = "fc_ival_round_to_float" 
external sys_single_precision_of_string: string -> float = 
    "fc_ival_single_precision_of_string"
(* TODO two functions above: declare "float", 
   must have separate version for bytecode, see OCaml manual *)

let max_single_precision_float = Int32.float_of_bits 0x7f7fffffl
let most_negative_single_precision_float = -. max_single_precision_float

type parsed_float = {
  f_nearest : float ;
  f_lower : float ;
  f_upper : float ;
}

let inf ~man_size ~max_exp = 
  let biggest_not_inf = ldexp (2.0 -. ldexp 1.0 (~- man_size)) max_exp in
  { 
    f_lower = biggest_not_inf ;
    f_nearest = infinity ;
    f_upper = infinity ;
  }

(* [s = num * 2^exp / den] hold *)
let make_float ~num ~den ~exp ~man_size ~min_exp ~max_exp = 
  assert (Integer.gt num Integer.zero);
  assert (Integer.gt den Integer.zero);
(*
  Format.printf "make_float: num den exp:@\n%a@\n@\n%a@\n@\n%d@.min_exp:%d max_exp:%d@."
    Datatype.Integer.pretty num Datatype.Integer.pretty den exp min_exp max_exp;
*)
  let size_bi = Integer.of_int man_size in
  let ssize_bi = Integer.of_int (succ man_size) in
  let min_exp = min_exp - man_size in

  let den = ref den in
  let exp = ref exp in
  while 
    Integer.ge num (Integer.shift_left !den ssize_bi)
    || !exp < min_exp
  do
    den := Integer.shift_left !den Integer.one;
    incr exp
  done;
  let den = !den in
  let shifted_den = Integer.shift_left den size_bi in
  let num = ref num in
  while 
    Integer.lt !num shifted_den && !exp > min_exp
  do
    num := Integer.shift_left !num Integer.one;
    decr exp
  done;
  let num = !num in
  let exp = !exp in
(* 
  Format.printf "make_float2: num den exp:@\n%a@\n@\n%a@\n@\n%d@."
    Datatype.Integer.pretty num Datatype.Integer.pretty den exp;
*)
  if exp > max_exp - man_size then inf ~man_size ~max_exp
  else
    let man,rem = Integer.e_div_rem num den in
    let rem2 = (* twice the remainder *)
      Integer.shift_left rem Integer.one
    in
    let man = Integer.to_int64 man in
(* Format.printf "pre-round: num den man rem:@\n%a@\n@\n%a@\n@\n%Ld@\n@\n%a@."
        Datatype.Integer.pretty num Datatype.Integer.pretty den
        man Datatype.Integer.pretty rem; *)
    let lowb = ldexp (Int64.to_float man) exp in
    if Integer.is_zero rem2 then {
      f_lower = lowb ;
      f_nearest = lowb ;
      f_upper = lowb ;
    } else
      let upb = ldexp (Int64.to_float (Int64.succ man)) exp in
      if Integer.lt rem2 den ||
        (Integer.equal rem2 den && (Int64.logand man Int64.one) = 0L)
      then {
	f_lower = lowb ;
	f_nearest = lowb ;
	f_upper = upb ;
      }
      else {
	f_lower = lowb ;
	f_nearest = upb ;
	f_upper = upb ;
      }

let reg_exp = "[eE][+]?\\(-?[0-9]+\\)"
let reg_dot = "[.]"
let reg_numopt = "\\([0-9]*\\)"
let reg_num = "\\([0-9]+\\)"

let numdotfrac = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt)
let numdotfracexp = Str.regexp (reg_numopt ^ reg_dot ^ reg_numopt ^ reg_exp)
let numexp = Str.regexp (reg_num ^ reg_exp)

exception Shortcut of parsed_float

let zero = { f_lower = 0.0 ; f_nearest = 0.0 ; f_upper = 0.0 }

(* [man_size] is the size of the mantissa, [min_exp] the frontier exponent
   between normalized and denormalized numbers *)
let parse_float ~man_size ~min_exp ~max_exp s =
  (*  Format.printf "parse: %s@." s; *)
  let match_exp group =
    let s = Str.matched_group group s in
    try
      int_of_string s
    with Failure _ ->
      (* Format.printf "Error in exponent: %s@." s; *)
      if s.[0] = '-'
      then raise (Shortcut { 
		    f_lower = 0.0 ; 
		    f_nearest = 0.0 ;
		    f_upper = ldexp 1.0 (min_exp - man_size) ;
		  })
      else raise (Shortcut (inf ~man_size ~max_exp))
  in
    try
      (* At the end of the function, [s = num * 2^exp / den] *)
      let num, den, exp =
	if Str.string_match numdotfracexp s 0
	then
	  let n = Str.matched_group 1 s in
	  let frac = Str.matched_group 2 s in
	  let len_frac = String.length frac in
	  let num = Integer.of_string (n ^ frac) in
	  let den = Integer.power_int_positive_int 5 len_frac in
          if Integer.is_zero num then raise (Shortcut zero);
	  let exp10 = match_exp 3
	  in
	  if exp10 >= 0
	  then
	    Integer.mul num (Integer.power_int_positive_int 5 exp10),
            den,
            exp10 - len_frac
	  else
	    num,
            Integer.mul den (Integer.power_int_positive_int 5 (~- exp10)),
            exp10 - len_frac
	else if Str.string_match numdotfrac s 0
	then 
	  let n = Str.matched_group 1 s in
	  let frac = Str.matched_group 2 s in
	  let len_frac = String.length frac in
	  Integer.of_string (n ^ frac),
          Integer.power_int_positive_int 5 len_frac,
          ~- len_frac
	else if Str.string_match numexp s 0
	then
	  let n = Str.matched_group 1 s in
	  let num = Integer.of_string n in
          if Integer.is_zero num then raise (Shortcut zero);
	  let exp10 = match_exp 2 in
	  if exp10 >= 0
	  then
	    Integer.mul num (Integer.power_int_positive_int 5 exp10),
            Integer.one,
            exp10
   	  else
	    num,
           (Integer.power_int_positive_int 5 (~- exp10)),
           exp10
	else (Format.printf "Could not parse floating point number %S@." s;
	      assert false)
      in
      if Integer.is_zero num 
      then zero
      else
	make_float ~num ~den ~exp ~man_size ~min_exp ~max_exp
    with Shortcut r -> r

let is_hex s =
  let l = String.length s in
  l >= 2 && s.[0] = '0' && (s.[1] = 'x' || s.[1] = 'X')

let opp_parse_float f =
  { f_lower = -. f.f_upper ; f_nearest = -. f.f_nearest ; f_upper = -. f.f_lower }

let rec single_precision_of_string s =
  if s.[0] = '-' then
     opp_parse_float (single_precision_of_string (String.sub s 1 (String.length s - 1)))
  else if is_hex s
  then
    try
      let f = sys_single_precision_of_string s in
      { f_lower = f ; f_nearest = f ; f_upper = f }
    with Failure _ ->
      Codex_log.fatal "could not parse single-precision float string: %s" s
  else (* decimal *)
    parse_float ~man_size:23 ~min_exp:(-126) ~max_exp:127 s

(* May raise Failure("float_of_string"). *)
let rec double_precision_of_string s = 
  if s.[0] = '-' then
    opp_parse_float (double_precision_of_string (String.sub s 1 (String.length s - 1)))
  else if is_hex s
  then 
    let f = float_of_string s in
    { f_lower = f ; f_nearest = f ; f_upper = f }
  else (* decimal *)
    parse_float ~man_size:52 ~min_exp:(-1022) ~max_exp:1023 s

(* let parse_kind kind string =
 *   match kind with
 *   | Cil_types.FFloat -> single_precision_of_string string
 *   | Cil_types.FDouble
 *   | Cil_types.FLongDouble -> double_precision_of_string string *)

(* let parse string =
 *   let l = String.length string - 1 in
 *   if l < 0
 *   then Codex_log.fatal ~current:true
 *       "Parsing an empty string as a floating-point constant."
 *   else
 *     let last = Transitioning.Char.uppercase_ascii string.[l] in
 *     let suffix, kind =
 *       match last with
 *       | 'F' -> true, Cil_types.FFloat
 *       | 'D' -> true, Cil_types.FDouble
 *       | 'L' -> true, Cil_types.FLongDouble
 *       | _ -> false, Cil_types.FDouble
 *     in
 *     let baseint = if suffix then String.sub string 0 l else string in
 *     try
 *       let basefloat = parse_kind kind baseint in
 *       kind, basefloat
 *     with Failure _ -> (\* should never happen, suffix already stripped *\)
 *       Kernel.fatal ~current:true
 *         "Unexpected error parsing floating-point constant: %s." string *)

let pretty_normal ~use_hex fmt f =
  let double_norm = Int64.shift_left 1L 52 in
  let double_mask = Int64.pred double_norm in
  let i = Int64.bits_of_float f in
  let s = 0L <> (Int64.logand Int64.min_int i) in
  let i = Int64.logand Int64.max_int i in
  let exp = Int64.to_int (Int64.shift_right_logical i 52) in
  let man = Int64.logand i double_mask in
  let s = if s then "-" else "" in
  if exp = 2047
  then begin
      if man = 0L
      then 
	Format.fprintf fmt "%sinf" s
      else
	Format.fprintf fmt "NaN"
    end
  else
  let firstdigit, exp =
    if exp <> 0
    then 1, (exp - 1023)
    else 0, -1022
  in
  if not use_hex
  then begin
      let firstdigit, man, exp =
	if 0 < exp && exp <= 12
	then begin
	    Int64.to_int
	      (Int64.shift_right_logical
		  (Int64.logor man double_norm)
		  (52 - exp)),
  	    Int64.logand (Int64.shift_left man exp) double_mask,
	    0
	  end
	else firstdigit, man, exp
      in
      let d =
	Int64.float_of_bits
	  (Int64.logor 0x3ff0000000000000L man)
      in
      let d, re =
	if d >= 1.5
	then d -. 1.5, 5000000000000000L
	else d -. 1.0, 0L
      in
      let d = d *. 1e16 in
      let decdigits = Int64.add re (Int64.of_float d) in
      if exp = 0 || (firstdigit = 0 && decdigits = 0L && exp = -1022)
      then
	Format.fprintf fmt "%s%d.%016Ld"
	  s
	  firstdigit
	  decdigits
      else
	Format.fprintf fmt "%s%d.%016Ld*2^%d"
	  s
	  firstdigit
	  decdigits
	  exp
    end
  else
    Format.fprintf fmt "%s0x%d.%013Lxp%d"
      s
      firstdigit
      man
      exp

let pretty fmt f =
  let use_hex = (* Kernel.FloatHex.get() in *)true in
  pretty_normal ~use_hex fmt f
  (* (\* should always arrive here with nearest_even *\)
   * if get_rounding_mode () <> FE_ToNearest then begin
   *   (assert false) "pretty: rounding mode (%s) <> FE_TONEAREST"
   *     (string_of_c_rounding_mode (get_rounding_mode ()));
   *   set_round_nearest_even();
   * end;
   * if use_hex || (Kernel.FloatNormal.get ())
   * then
   *   pretty_normal ~use_hex fmt f
   * else begin
   *     let r = Format.sprintf "%.*g" 12 f in
   *     if (String.contains r '.' || String.contains r 'e' ||
   *            String.contains r 'E')
   *       || (match classify_float f with
   *       | FP_normal | FP_subnormal | FP_zero -> false
   *       | FP_infinite | FP_nan -> true)
   *     then Format.pp_print_string fmt r
   *     else Format.fprintf fmt "%s." r
   *   end *)


type sign = Neg | Pos

exception Float_Non_representable_as_Int64 of sign

(* If the argument [x] is not in the range [min_64_float, 2*max_64_float],
   raise Float_Non_representable_as_Int64. This is the most reasonable as
   a floating-point number may represent an exponentially large integer. *)
let truncate_to_integer =
  let min_64_float = -9.22337203685477581e+18 
           (* Int64.to_float (-0x8000000000000000L) *) 
  in
  let max_64_float = 9.22337203685477478e+18 
(*    let open Int64 in
    float_of_bits (pred (bits_of_float (to_float max_int))) *)
  in
  fun x ->
    let max_64_float = (fun x -> x)(* Extlib.id *) max_64_float in
    if x < min_64_float
    then raise (Float_Non_representable_as_Int64 Neg);
    if x > (max_64_float +. max_64_float)
    then raise (Float_Non_representable_as_Int64 Pos);
    if x <= max_64_float then
      Integer.of_int64 (Int64.of_float x)
    else 
      Integer.add 
	(Integer.of_int64 (Int64.of_float (x +. min_64_float)))
	(Integer.two_power_of_int 63)

let bits_of_max_double =
  Integer.of_int64 (Int64.bits_of_float max_float)
let bits_of_most_negative_double =
  Integer.of_int64 (Int64.bits_of_float (-. max_float))

(** See e.g. http://www.h-schmidt.net/FloatConverter/IEEE754.html *)
let bits_of_max_float = Integer.of_int64 0x7F7FFFFFL
let bits_of_most_negative_float =
  let v = Int64.of_int32 0xFF7FFFFFl in(* cast to int32 to get negative value *)
  Integer.of_int64 v

external fround: float -> float = "fc_ival_c_round"
external trunc: float -> float = "fc_ival_c_trunc"

(** Single-precision (32-bit) functions. We round the result computed
    as a double, since float32 functions are rarely precise. *)

external expf: float -> float = "fc_ival_c_expf"
external logf: float -> float = "fc_ival_c_logf"
external log10f: float -> float = "fc_ival_c_log10f"
external powf: float -> float -> float = "fc_ival_c_powf"
external sqrtf: float -> float = "fc_ival_c_sqrtf"
external fmodf: float -> float -> float = "fc_ival_c_fmodf"
external cosf: float -> float = "fc_ival_c_cosf"
external sinf: float -> float = "fc_ival_c_sinf"
external atan2f: float -> float -> float = "fc_ival_c_atan2f"


(** C math-like functions *)

let isnan f =
  match classify_float f with
  | FP_nan -> true
  | _ -> false

let isfinite f =
  match classify_float f with
  | FP_nan | FP_infinite -> false
  | _ -> true

let min_denormal = Int64.float_of_bits 1L
let neg_min_denormal = -. min_denormal
let min_single_precision_denormal = Int32.float_of_bits 1l
let neg_min_single_precision_denormal = -. min_single_precision_denormal

(* auxiliary functions for nextafter/nextafterf *)
let min_denormal_float ~is_f32 =
  if is_f32 then min_single_precision_denormal else min_denormal
let nextafter_aux ~is_f32 fincr fdecr x y =
  if x = y (* includes cases "(0.0, -0.0) => -0.0" and its symmetric *)
  then y
  else if isnan x || isnan y then nan
  else if x = 0.0 (* or -0.0 *) then
    if x < y then min_denormal_float is_f32 else -. (min_denormal_float is_f32)
    (* the following conditions might be simpler if we had unsigned ints
       (uint32/uint64) *)
  else if x = neg_infinity (* && y = neg_infinity *) then fdecr x
  else if (x < y && x > 0.0) || (x > y && x < 0.0) then fincr x else fdecr x
let incr_f64 f =
  Int64.float_of_bits (Int64.succ (Int64.bits_of_float f))
let decr_f64 f =
  if f = infinity then max_float
  else Int64.float_of_bits (Int64.pred (Int64.bits_of_float f))
let incr_f32 f =
  if f = neg_infinity then most_negative_single_precision_float
  else Int32.float_of_bits (Int32.succ (Int32.bits_of_float f))
let decr_f32 f =
  if f = infinity then max_single_precision_float
  else Int32.float_of_bits (Int32.pred (Int32.bits_of_float f))

let nextafter x y =
  nextafter_aux ~is_f32:false incr_f64 decr_f64 x y

let nextafterf x y =
  nextafter_aux ~is_f32:true incr_f32 decr_f32 x y


(*
Local Variables:
compile-command: "make -C ../../.. byte"
End:
*)