bin_float.hl

(* ========================================================================== *)
(*      Formal verification of FPTaylor certificates                          *)
(*                                                                            *)
(*      Author: Alexey Solovyev, University of Utah                           *)
(*                                                                            *)
(*      This file is distributed under the terms of the MIT licence           *)
(* ========================================================================== *)

(* -------------------------------------------------------------------------- *)
(* Binary floating-point numbers                                              *)
(* Note: requires the nonlinear inequality verification tool                  *)
(*       https://github.com/monadius/formal_ineqs                             *)
(* -------------------------------------------------------------------------- *)

needs "ipow.hl";;
needs "arith/arith_float.hl";;
needs "arith/float_pow.hl";;
needs "arith/more_float.hl";;

module Bin_float = struct

open Num;;
open Arith_nat;;
open Arith_float;;
open More_float;;
open Float_pow;;
open Ipow;;
open Misc_functions;;
open Misc_vars;;

prioritize_real();;

(* --------------------------------------------- *)
(* Definitions                                   *)
(* --------------------------------------------- *)

let ipow2 = new_definition `ipow2 (e:int) = &2 ipow e`;;

let bin_float = new_definition
  `bin_float s m e = (if s then (-- &1) else &1) * &m * ipow2 e`;;

(* --------------------------------------------- *)
(* Basic properties                              *)
(* --------------------------------------------- *)

let ipow2_0 = prove
  (`ipow2 (&0) = &1`,
   REWRITE_TAC[ipow2; IPOW_0]);;

let ipow2_add = prove
  (`!e1 e2. ipow2 (e1 + e2) = ipow2 e1 * ipow2 e2`,
   ASM_SIMP_TAC[ipow2; IPOW_ADD; REAL_ARITH `~(&2 = &0)`]);;

let ipow2_sub = prove
  (`!e1 e2. ipow2 (e1 - e2) = ipow2 e1 * ipow2 (--e2)`,
   ASM_SIMP_TAC[ipow2; IPOW_SUB; REAL_ARITH `~(&2 = &0)`]);;

let ipow2_pos = prove
  (`!e. &0 < ipow2 e`,
   ASM_SIMP_TAC[ipow2; IPOW_LT_0; REAL_ARITH `&0 < &2`]);;

let ipow2_pos_le = prove
  (`!e. &0 <= ipow2 e`,
   SIMP_TAC[ipow2_pos; REAL_LT_IMP_LE]);;

let ipow2_pos_num = prove
  (`!n. ipow2 (&n) = &2 pow n`,
   REWRITE_TAC[ipow2; IPOW_NUM]);;

let ipow2_neg_num = prove
  (`!n. ipow2 (-- &n) = inv (&2) pow n`,
   REWRITE_TAC[ipow2; IPOW_NEG; INV_IPOW; IPOW_NUM]);;

let bin_float_T_eq = prove
  (`!m e. bin_float T m e = --bin_float F m e`,
   REWRITE_TAC[bin_float; REAL_MUL_LNEG]);;

let pos_bin_float = prove
  (`!m e. &0 <= bin_float F m e`,
   REWRITE_TAC[bin_float; REAL_MUL_LID] THEN REPEAT GEN_TAC THEN
     MATCH_MP_TAC REAL_LE_MUL THEN REWRITE_TAC[REAL_POS; ipow2_pos_le]);;

let neg_bin_float = prove
  (`!m e. bin_float T m e <= &0`,
   REWRITE_TAC[bin_float_T_eq; REAL_NEG_LE0; pos_bin_float]);;

let bin_float_pos = prove
  (`interval_arith (&m * ipow2 e) bounds
     <=> interval_arith (bin_float F m e) bounds`,
   REWRITE_TAC[bin_float; REAL_MUL_LID]);;

let bin_float_neg = prove
  (`interval_arith (-- (&m * ipow2 e)) bounds
     <=> interval_arith (bin_float T m e) bounds`,
   REWRITE_TAC[bin_float; REAL_MUL_LNEG; REAL_MUL_LID]);;

(* --------------------------------------------- *)
(* Constuction/destruction                       *)
(* --------------------------------------------- *)

let mk_ipow2 =
  let ipow2_const = `ipow2` in
    fun n ->
      let n_tm = mk_intconst n in
	mk_comb (ipow2_const, n_tm);;

let dest_bin_float tm =
  match tm with
    | Comb (Comb (Comb (Const ("bin_float", _), Const (s_name, _)), m_tm), e_tm) ->
	let s = (s_name = "T") in
	  s, m_tm, e_tm
    | _ -> error "dest_bin_float" [tm] [];;

let mk_small_int =
  let neg_int_const = `(--) : int->int` and
      int_const = `(&) : num->int` in
    fun n ->
      if n >= 0 then
	mk_comb (int_const, mk_small_numeral n)
      else
	mk_comb (neg_int_const, mk_comb (int_const, mk_small_numeral (-n)));;

let dest_fp =
  let rec dest f (m, e) =
    if f <= 0.0 then (m, e)
    else
      let d = int_of_float f in
      let f' = ldexp (f -. float_of_int d) 1 in
	dest f' (Int 2 */ m +/ Int d, e - 1) in
    fun f ->
      let m, e = frexp f in
      let m, sign = if m < 0.0 then -.m, true else m, false in
	sign, dest m (Int 0, e + 1);;

let mk_bin_float =
  let bin_float_const = `bin_float` in
    fun f ->
      let s, (m, e) = dest_fp f in
      let s_tm = if s then t_const else f_const in
      let m_tm = mk_numeral m in
      let e_tm = mk_small_int e in
	mk_comb (mk_comb (mk_comb (bin_float_const, s_tm), m_tm), e_tm);;

(* --------------------------------------------- *)
(* Interval arithmetic                           *)
(* --------------------------------------------- *)

let ipow2_interval =
  let ipow2_pos = (SYM o SPEC_ALL) ipow2_pos_num and
      ipow2_neg = (SYM o SPEC_ALL) ipow2_neg_num in
  let inv2_interval = 
    GEN_REWRITE_RULE LAND_CONV [GSYM float_inv2_th] More_float.inv2_interval in
    fun pp e_tm ->
      match e_tm with
	| Comb (Const ("int_of_num", _), n_tm) ->
	    let n = dest_small_numeral n_tm in
	      ONCE_REWRITE_RULE[ipow2_pos] (float_interval_pow pp n two_interval)
	| Comb (Const ("int_neg", _), Comb (Const ("int_of_num", _), n_tm)) ->
	    let n = dest_small_numeral n_tm in
	      ONCE_REWRITE_RULE[ipow2_neg] (float_interval_pow pp n inv2_interval)
	| _ -> failwith "ipow2_interval";;

let bin_float_interval =
  let e_var_int = `e : int` in
    fun pp tm ->
      let s, m_tm, e_tm = dest_bin_float tm in
      let m_th = mk_float_interval_num (dest_numeral m_tm) in
      let ipow2_th = ipow2_interval pp e_tm in
      let th1 = float_interval_mul pp m_th ipow2_th in
	if s then
	  let th2 = float_interval_neg th1 in
	  let bounds = rand (concl th2) in
	  let th0 = INST[m_tm, m_var_num; e_tm, e_var_int;
			 bounds, bounds_var] bin_float_neg in
	    EQ_MP th0 th2
	else
	  let bounds = rand (concl th1) in
	  let th0 = INST[m_tm, m_var_num; e_tm, e_var_int;
			 bounds, bounds_var] bin_float_pos in
	    EQ_MP th0 th1;;

(* --------------------------------------------- *)
(* Rational conversion                           *)
(* --------------------------------------------- *)

let bin_float_rat_conv =
    REWRITE_CONV[GSYM ipow2; bin_float; ipow2_neg_num; ipow2_pos_num] THENC 
      REAL_RAT_REDUCE_CONV;;

end;;