(* Dependencies *) needs "../formal_lp/arith/misc.hl";; needs "../formal_lp/arith/arith_options.hl";; (* Natural numbers *) module type Informal_nat_sig = sig type nat val mk_nat : num -> nat val mk_small_nat : int -> nat val dest_nat : nat -> num val eq_nat : nat -> nat -> bool val suc_nat : nat -> nat val pre_nat : nat -> nat val eq0_nat : nat -> bool val gt0_nat : nat -> bool val lt_nat : nat -> nat -> bool val le_nat : nat -> nat -> bool val add_nat : nat -> nat -> nat val sub_nat : nat -> nat -> nat (* If sub_and_le_nat m n = (m - n, true) if n <= m; (n - m, false) if m < n *) val sub_and_le_nat : nat -> nat -> nat * bool val mul_nat : nat -> nat -> nat val div_nat : nat -> nat -> nat val even_nat : nat -> bool val odd_nat : nat -> bool (* normalize_nat m = (n, e) s.t. m = n * base^e, e >= 0 *) val normalize_nat : nat -> nat * int val denormalize_nat : nat * int -> nat (* hi_nat p m = (n, e) s.t. m <= n * base^e and n contains at most p "digits" *) val hi_nat : int -> nat -> nat * int val hi_lt_nat : int -> nat -> nat * int (* lo_nat p m = (n, e) s.t. n * base^e <= m and n contains at most p "digits" *) val lo_nat : int -> nat -> nat * int end;; module Informal_nat : Informal_nat_sig = struct open Arith_misc;; open Arith_options;; open Big_int;; type nat = big_int;; let mk_nat n = let result = big_int_of_num n in if sign_big_int result < 0 then zero_big_int else result;; let mk_small_nat n = if n < 0 then zero_big_int else big_int_of_int n;; let dest_nat = num_of_big_int;; let eq_nat = eq_big_int;; let suc_nat = succ_big_int;; let pre_nat n = let result = pred_big_int n in if sign_big_int result < 0 then zero_big_int else result;; let eq0_nat n = sign_big_int n = 0;; let gt0_nat n = sign_big_int n > 0;; let lt_nat = lt_big_int;; let le_nat = le_big_int;; let add_nat = add_big_int;; let sub_nat m n = let result = sub_big_int m n in if sign_big_int result < 0 then zero_big_int else result;; let sub_and_le_nat m n = let result = sub_big_int m n in if sign_big_int result >= 0 then (result, true) else (abs_big_int result, false);; let mul_nat = mult_big_int;; let div_nat = div_big_int;; let two_big_int = big_int_of_int 2;; let even_nat n = sign_big_int (mod_big_int n two_big_int) = 0;; let odd_nat n = sign_big_int (mod_big_int n two_big_int) > 0;; (*******************************) (* num_exp *) let base_nat = mk_small_nat base;; (* normalize_nat m = (n, e) s.t. m = n * base^e, e >= 0 *) let normalize_nat = let rec normalize n e = let q, r = quomod_big_int n base_nat in if sign_big_int r > 0 then (n, e) else normalize q (succ e) in fun n -> if sign_big_int n = 0 then (n, 0) else normalize n 0;; let denormalize_nat (n, e) = mult_big_int n (power_int_positive_int base e);; let lo_nat pp = let max = power_int_positive_int base pp in let rec lo m e = if lt_big_int m max then (m, e) else let q = div_big_int m base_nat in lo q (succ e) in fun m -> if sign_big_int m = 0 then (m, 0) else let n1, e1 = lo m 0 in let n, e2 = normalize_nat n1 in n, e1 + e2;; let hi_nat pp = if pp <= 0 then failwith "hi_nat: pp <= 0" else let max = power_int_positive_int base pp in let rec hi m e = if lt_big_int m max then (m, e) else let q, r = quomod_big_int m base_nat in if sign_big_int r = 0 then hi q (succ e) else hi (succ_big_int q) (succ e) in fun m -> if sign_big_int m = 0 then (m, 0) else let n1, e1 = hi m 0 in let n, e2 = normalize_nat n1 in n, e1 + e2;; let hi_lt_nat pp m = hi_nat pp (succ_big_int m);; end;; (* Floating point numbers *) module type Informal_float_sig = sig type ifloat val mk_float : num -> int -> ifloat val mk_num_float : num -> ifloat val mk_small_num_float : int -> ifloat val dest_float : ifloat -> bool * num * int val sign_float : ifloat -> bool (* Compares representations, not numbers themselves *) val eq_float : ifloat -> ifloat -> bool val lo_float : int -> ifloat -> ifloat val hi_float : int -> ifloat -> ifloat val neg_float : ifloat -> ifloat val abs_float : ifloat -> ifloat val lt0_float : ifloat -> bool val gt0_float : ifloat -> bool val le0_float : ifloat -> bool val ge0_float : ifloat -> bool val lt_float : ifloat -> ifloat -> bool val le_float : ifloat -> ifloat -> bool val min_float : ifloat -> ifloat -> ifloat val max_float : ifloat -> ifloat -> ifloat val mul_float_eq : ifloat -> ifloat -> ifloat val mul_float_lo : int -> ifloat -> ifloat -> ifloat val mul_float_hi : int -> ifloat -> ifloat -> ifloat val div_float_lo : int -> ifloat -> ifloat -> ifloat val div_float_hi : int -> ifloat -> ifloat -> ifloat val add_float_lo : int -> ifloat -> ifloat -> ifloat val add_float_hi : int -> ifloat -> ifloat -> ifloat val sub_float_lo : int -> ifloat -> ifloat -> ifloat val sub_float_hi : int -> ifloat -> ifloat -> ifloat val sqrt_float_lo : int -> ifloat -> ifloat val sqrt_float_hi : int -> ifloat -> ifloat end;; module Informal_float : Informal_float_sig = struct open Arith_options;; open Informal_nat;; type ifloat = bool * nat * int;; (* Creates a non-negative float *) let mk_float n e : ifloat = false, mk_nat n, e;; let mk_num_float n = false, mk_nat n, min_exp;; let mk_small_num_float n = false, mk_small_nat n, min_exp;; let dest_float ((s, n, e) : ifloat) = s, dest_nat n, e;; let sign_float ((s,_,_) : ifloat) = s;; let eq_float (s1,n1,e1) (s2,n2,e2) = s1 = s2 && eq_nat n1 n2 && e1 = e2;; let lo_float pp (s,n,e) = let n1, e1 = if s then hi_nat pp n else lo_nat pp n in (s, n1, e + e1);; let hi_float pp (s,n,e) = let n1, e1 = if s then lo_nat pp n else hi_nat pp n in (s, n1, e + e1);; (* Auxiliary num_exp functions *) let num_exp_add = let (+) = add_nat in fun (n1,e1) (n2,e2) -> if e1 <= e2 then n1 + denormalize_nat (n2, e2 - e1), e1 else n2 + denormalize_nat (n1, e1 - e2), e2;; (* Returns (n,e),true if (n1,e1) >= (n2,e2) and (n,e) = (n1,e1) - (n2,e2) Returns (n,e),false if (n1,e1) <= (n2,e2) and (n,e) = (n2,e2) - (n1,e1) *) let num_exp_sub = let (--) = sub_and_le_nat in fun (n1,e1) (n2,e2) -> if e2 <= e1 then let a = denormalize_nat (n1, e1 - e2) and b = n2 in let sub, flag = a -- b in (sub, e2), flag else let a = n1 and b = denormalize_nat (n2, e2 - e1) in let sub, flag = a -- b in (sub, e1), flag;; let num_exp_le = let (<=/) = le_nat in fun (n1,e1) (n2,e2) -> if e1 <= e2 then n1 <=/ denormalize_nat (n2, e2 - e1) else denormalize_nat (n1, e1 - e2) <=/ n2;; let num_exp_lt = let ( if e1 <= e2 then n1 interval val mk_num_interval : num -> interval val mk_small_num_interval : int -> interval val dest_interval : interval -> Informal_float.ifloat * Informal_float.ifloat val round_interval : int -> interval -> interval val neg_interval : interval -> interval val mul_interval : int -> interval -> interval -> interval val div_interval : int -> interval -> interval -> interval val add_interval : int -> interval -> interval -> interval val sub_interval : int -> interval -> interval -> interval val sqrt_interval : int -> interval -> interval val inv_interval : int -> interval -> interval val pow_interval : int -> int -> interval -> interval (* Computes max(-lo, hi) *) val abs_interval : interval -> Informal_float.ifloat end;; module Informal_interval : Informal_interval_sig = struct open Informal_float;; type interval = ifloat * ifloat;; let mk_interval (lo,hi) = if lt_float hi lo then failwith "mk_interval: hi < lo" else (lo,hi);; let mk_num_interval n = let f = mk_num_float n in (f, f);; let mk_small_num_interval n = let f = mk_small_num_float n in (f, f);; let one_interval = mk_small_num_interval 1;; let two_interval = mk_small_num_interval 2;; let dest_interval ((lo,hi) : interval) = (lo,hi);; let round_interval pp (lo,hi) = (lo_float pp lo, hi_float pp hi);; let neg_interval (lo,hi) = (neg_float hi, neg_float lo);; let abs_interval (lo,hi) = max_float hi (neg_float lo);; let add_interval pp (lo1,hi1) (lo2,hi2) = (add_float_lo pp lo1 lo2, add_float_hi pp hi1 hi2);; let sub_interval pp (lo1,hi1) (lo2,hi2) = (sub_float_lo pp lo1 hi2, sub_float_hi pp hi1 lo2);; let sqrt_interval pp (lo,hi) = if sign_float lo then failwith "sqrt_interval: negative lower bound" else (sqrt_float_lo pp lo, sqrt_float_hi pp hi);; (* mul *) let mul_interval pp (l_lo,l_hi) (r_lo,r_hi) = let s1 = sign_float l_lo and s2 = sign_float l_hi and s3 = sign_float r_lo and s4 = sign_float r_hi in if s1 <> s2 && s3 <> s4 then if not s1 or not s3 then failwith "mul_interval: FT interval" else let lo1, lo2 = mul_float_lo pp l_lo r_hi, mul_float_lo pp l_hi r_lo and hi1, hi2 = mul_float_hi pp l_lo r_lo, mul_float_hi pp l_hi r_hi in (min_float lo1 lo2, max_float hi1 hi2) else let lo1, lo2, hi1, hi2 = if s1 <> s2 then if not s1 then failwith "mul_interval: FT interval" else if not s3 then l_lo, r_hi, l_hi, r_hi else l_hi, r_lo, l_lo, r_lo else if s3 <> s4 then if not s3 then failwith "mul_interval: FT interval" else if not s1 then l_hi, r_lo, l_hi, r_hi else l_lo, r_hi, l_lo, r_lo else if not s1 then if not s3 then l_lo, r_lo, l_hi, r_hi else l_hi, r_lo, l_lo, r_hi else if not s3 then l_lo, r_hi, l_hi, r_lo else l_hi, r_hi, l_lo, r_lo in (mul_float_lo pp lo1 lo2, mul_float_hi pp hi1 hi2);; (* div *) let div_interval pp (l_lo,l_hi) (r_lo,r_hi) = let s1 = sign_float l_lo and s2 = sign_float l_hi and s3 = sign_float r_lo and s4 = sign_float r_hi in if s3 <> s4 then failwith "div_interval: division by an interval containing 0" else let lo1, lo2, hi1, hi2 = if s1 = s2 then if not s1 then if not s3 then l_lo, r_hi, l_hi, r_lo else l_hi, r_hi, l_lo, r_lo else if not s3 then l_lo, r_lo, l_hi, r_hi else l_hi, r_lo, l_lo, r_hi else if not s1 then failwith "div_interval: FT interval" else if not s3 then l_lo, r_lo, l_hi, r_lo else l_hi, r_hi, l_lo, r_hi in (div_float_lo pp lo1 lo2, div_float_hi pp hi1 hi2);; (* inv *) let inv_interval pp int = div_interval pp one_interval int;; (* pow *) let pow_interval pp n int = let rec pow n = if n <= 0 then one_interval else if n = 1 then int else let i2 = pow (n - 1) in mul_interval pp int i2 in pow n;; (* Arith_misc.gen_pow (mul_interval pp) one_interval n;; *) end;; (* atn *) module type Informal_atn_sig = sig val atn_interval : int -> Informal_interval.interval -> Informal_interval.interval val acs_interval : int -> Informal_interval.interval -> Informal_interval.interval val pi_approx_array : Informal_interval.interval array val pi2_approx_array : Informal_interval.interval array end;; module Informal_atn : Informal_atn_sig = struct open Informal_float;; open Informal_interval;; let rec poly_f_interval pp l x = if length l = 0 then failwith "poly_f_interval: an empty coefficient list" else let first = hd l in if length l = 1 then first else let r = poly_f_interval pp (tl l) x in add_interval pp first (mul_interval pp x r);; let poly_f_even_interval pp l x = let xx = mul_interval pp x x in poly_f_interval pp l xx;; let poly_f_odd_interval pp l x = let even = poly_f_even_interval pp l x in mul_interval pp x even;; let halfatn_interval pp x = let xx = mul_interval pp x x in let one_xx = add_interval pp one_interval xx in let sqrt = sqrt_interval pp one_xx in let one_sqrt = add_interval pp sqrt one_interval in div_interval pp x one_sqrt;; let halfatn4_interval pp x = (halfatn_interval pp o halfatn_interval pp o halfatn_interval pp o halfatn_interval pp) x;; (* Computes an interval for 16 * sum(0..n) (halfatn4_co x) *) let atn_sum_interval = let interval_16 = mk_small_num_interval 16 in fun pp l x -> let halfatn4 = halfatn4_interval pp x in let poly = poly_f_odd_interval pp l halfatn4 in mul_interval pp interval_16 poly;; let atn0_interval pp l eps x = let sum = atn_sum_interval pp l x in let a, b = dest_interval sum in let _, d = dest_interval eps in let hi = add_float_hi pp b d in let lo = sub_float_lo pp a d in mk_interval (lo, hi);; (* Computes an interval for 2 ^ -(6n + 5) *) let compute_eps1 pp n = let pow = pow_interval pp (6 * n + 5) two_interval in inv_interval pp pow;; let mk_atn_co_table pp n = let get_val k = let l = if (k land 1) = 0 then one_interval else neg_interval (one_interval) in let r = mk_small_num_interval (2 * k + 1) in div_interval pp l r in map get_val (0--n);; (* Lookup tables *) let n_of_p pp = let x = (float_of_int (pp + 1) *. log (float_of_int Arith_options.base) /. log (2.0) -. 5.0) /. 6.0 in let n = (int_of_float o ceil) x in if n < 1 then 1 else n;; let atn_co_array = Array.init 21 (fun i -> mk_atn_co_table (i + 1) (n_of_p i));; let eps1_array = Array.init 21 (fun i -> compute_eps1 (i + 1) (n_of_p i));; let atn_interval pp x = atn0_interval pp atn_co_array.(pp) eps1_array.(pp) x;; (* pi approximation *) let pi_approx_array, pi2_approx_array = let pp = 20 in let x = one_interval in let r1 = atn_interval pp x in let r2 = mul_interval pp (mk_small_num_interval 4) r1 in let float_pi = r2 in let float_pi2 = div_interval pp float_pi two_interval in let pi_int0 = mk_small_num_interval 0 in let pi2_int0 = pi_int0 in Array.init 19 (fun i -> if i = 0 then pi_int0 else round_interval i float_pi), Array.init 19 (fun i -> if i = 0 then pi2_int0 else round_interval i float_pi2);; (* acs *) let acs0_interval pp l eps1 x = let int1 = sub_interval pp one_interval (mul_interval pp x x) in let int2 = div_interval pp x (sqrt_interval pp int1) in let atn_int = atn0_interval pp l eps1 int2 in sub_interval pp pi2_approx_array.(pp + 1) atn_int;; let acs_interval pp x = acs0_interval pp atn_co_array.(pp) eps1_array.(pp) x;; end;; (* needs "../formal_lp/arith/float_atn.hl";; open Float_atn;; open Informal_atn;; open Arith_float;; open Informal_interval;; let n1 = 111 and n2 = 33;; let a_th = mk_float_interval_small_num n1 and b_th = mk_float_interval_small_num n2;; let a = mk_small_num_interval n1 and b = mk_small_num_interval n2;; let dest_int i = let lo, hi = dest_interval i in Informal_float.dest_float lo, Informal_float.dest_float hi;; let pp = 10;; float_interval_mul pp a_th b_th;; dest_int (mul_interval pp a b);; float_interval_div pp a_th b_th;; dest_int (div_interval pp a b);; let mk_f = (fst o dest_pair o rand o concl o mk_float_interval_small_num);; let mk_neg_f = (fst o dest_pair o rand o concl o float_interval_neg o mk_float_interval_small_num);; let mk_if n = mk_float (Num.num_of_int n) Arith_options.min_exp;; let n1 = 3 and n2 = 121;; let f1_tm = mk_f n1 and f2_tm = mk_f n2;; let f1 = (mk_if n1) and f2 = (mk_if n2);; let pp = 20;; float_sqrt_lo pp f1_tm;; dest_float (sqrt_float_lo pp f1);; float_sqrt_lo pp f2_tm;; dest_float (sqrt_float_lo pp f2);; float_sqrt_hi pp f1_tm;; dest_float (sqrt_float_hi pp f1);; float_sqrt_hi pp f2_tm;; dest_float (sqrt_float_hi pp f2);; open Arith_misc;; test 100 (float_sqrt_lo pp) f1_tm;; test 100 (sqrt_float_lo pp) f1;; dest_float (sqrt_float_lo pp f1);; float_le0 f1_tm;; le0_float f1;; float_le0 f2_tm;; le0_float f2;; float_lt f2_tm f1_tm;; lt_float f2 f1;; float_le f2_tm f1_tm;; le_float f2 f1;; float_mul_eq f1_tm f2_tm;; dest_float (mul_float_eq f1 f2);; let pp = 2;; float_add_lo pp f1_tm f2_tm;; dest_float (add_float_lo pp f1 f2);; float_add_hi pp f1_tm f2_tm;; dest_float (add_float_hi pp f1 f2);; float_sub_lo pp f1_tm f2_tm;; dest_float (sub_float_lo pp f1 f2);; float_sub_hi pp f1_tm f2_tm;; dest_float (sub_float_hi pp f1 f2);; float_mul_lo pp f1_tm f2_tm;; dest_float (mul_float_lo pp f1 f2);; float_mul_hi pp f1_tm f2_tm;; dest_float (mul_float_hi pp f1 f2);; float_div_lo pp f1_tm f2_tm;; dest_float (div_float_lo pp f1 f2);; float_div_hi pp f1_tm f2_tm;; dest_float (div_float_hi pp f1 f2);; open Arith_misc;; test 1000 (float_div_lo pp f1_tm) f2_tm;; test 1000 (div_float_lo pp f1) f2;; *)