(* 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 (</) = lt_nat in
    fun (n1,e1) (n2,e2) ->
      if e1 <= e2 then
	n1 </ denormalize_nat (n2, e2 - e1)
      else
	denormalize_nat (n1, e1 - e2) </ n2;;



(* neg *)

let neg_float (s,n,e) = (not s, n, e);;


(* abs *)

let abs_float (_,n,e) = (false, n, e);;
  


(* lt0, gt0 *)

let lt0_float (s,n,e) =
  if not s then false else gt0_nat n;;

let gt0_float (s,n,e) =
  if s then false else gt0_nat n;;


(* le0, ge0 *)

let le0_float (s,n,e) =
  if s then true else eq0_nat n;;

let ge0_float (s,n,e) =
  if s then eq0_nat n else true;;


(* lt *)

let lt_float (s1,n1,e1) (s2,n2,e2) =
  if not s1 then
    if s2 then false else num_exp_lt (n1,e1) (n2,e2)
  else
    if s2 then num_exp_lt (n2,e2) (n1,e1) 
    else
      (* TF *)
      if eq0_nat n1 then gt0_nat n2 else true;;


let le_float (s1,n1,e1) (s2,n2,e2) =
  if s1 then
    if s2 then num_exp_le (n2,e2) (n1,e1) else true
  else
    if not s2 then num_exp_le (n1,e1) (n2,e2)
    else
      (* FT *)
      if eq0_nat n2 then eq0_nat n1 else false;;
	
      
    
(* min, max *)

let min_float f1 f2 =
  if le_float f1 f2 then f1 else f2;;

let max_float f1 f2 =
  if le_float f1 f2 then f2 else f1;;


(* mul *)

let badd b1 b2 =
  if b1 then not b2 else b2;;


let mul_float_eq (s1,n1,e1) (s2,n2,e2) =
  let s = badd s1 s2 in
  let n = mul_nat n1 n2 in
  let e = e1 + e2 - min_exp in
    if e < 0 then
      failwith "mul_float_eq: underflow"
    else
      (s, n, e);;


let mul_float_lo pp (s1,n1,e1) (s2,n2,e2) =
  let s = badd s1 s2 in
  let n' = mul_nat n1 n2 in
  let n, e' = if s1 = s2 then lo_nat pp n' else hi_nat pp n' in
  let e = e1 + e2 + e' - min_exp in
    if e < 0 then
      failwith "mul_float_lo: underflow"
    else
      (s, n, e);;


let mul_float_hi pp (s1,n1,e1) (s2,n2,e2) =
  let s = badd s1 s2 in
  let n' = mul_nat n1 n2 in
  let n, e' = if s1 = s2 then hi_nat pp n' else lo_nat pp n' in
  let e = e1 + e2 + e' - min_exp in
    if e < 0 then
      failwith "mul_float_hi: underflow"
    else
      (s, n, e);;


(* div *)

let div_float_lo pp (s1,n1,e1) (s2,n2,e2) =
  let s = badd s1 s2 in
  let k = 2 * pp in
  let nn1 = denormalize_nat (n1, k) in
  let n' = div_nat nn1 n2 in
  let n, e' = if s1 = s2 then lo_nat pp n' else hi_lt_nat pp n' in
  let e = min_exp + e' + e1 - e2 - k in
    if e < 0 then
      failwith "div_float_lo: underflow"
    else
      (s, n, e);;


let div_float_hi pp (s1,n1,e1) (s2,n2,e2) =
  let s = badd s1 s2 in
  let k = 2 * pp in
  let nn1 = denormalize_nat (n1, k) in
  let n' = div_nat nn1 n2 in
  let n, e' = if s1 = s2 then hi_lt_nat pp n' else lo_nat pp n' in
  let e = min_exp + e' + e1 - e2 - k in
    if e < 0 then
      failwith "div_float_hi: underflow"
    else
      (s, n, e);;


(* add *)

let add_float_lo pp (s1,n1,e1) (s2,n2,e2) =
  if s1 = s2 then
    let n', e' = num_exp_add (n1,e1) (n2,e2) in
    let n, e'' = if s1 then hi_nat pp n' else lo_nat pp n' in
      (s1, n, e' + e'')
  else
    if s1 then
      let (n', e'), flag = num_exp_sub (n2,e2) (n1,e1) in
	if flag then
	  let n, e'' = lo_nat pp n' in
	    (false, n, e' + e'')
	else
	  let n, e'' = hi_nat pp n' in
	    (true, n, e' + e'')
    else
      let (n', e'), flag = num_exp_sub (n1,e1) (n2,e2) in
	if flag then
	  let n, e'' = lo_nat pp n' in
	    (false, n, e' + e'')
	else
	  let n, e'' = hi_nat pp n' in
	    (true, n, e' + e'');;


let add_float_hi pp (s1,n1,e1) (s2,n2,e2) =
  if s1 = s2 then
    let n', e' = num_exp_add (n1,e1) (n2,e2) in
    let n, e'' = if s1 then lo_nat pp n' else hi_nat pp n' in
      (s1, n, e' + e'')
  else
    if s1 then
      let (n', e'), flag = num_exp_sub (n2,e2) (n1,e1) in
	if flag then
	  let n, e'' = hi_nat pp n' in
	    (false, n, e' + e'')
	else
	  let n, e'' = lo_nat pp n' in
	    (true, n, e' + e'')
    else
      let (n', e'), flag = num_exp_sub (n1,e1) (n2,e2) in
	if flag then
	  let n, e'' = hi_nat pp n' in
	    (false, n, e' + e'')
	else
	  let n, e'' = lo_nat pp n' in
	    (true, n, e' + e'');;
      

(* sub *)

let sub_float_lo pp f1 f2 = add_float_lo pp f1 (neg_float f2);;
let sub_float_hi pp f1 f2 = add_float_hi pp f1 (neg_float f2);;


(* sqrt *)

let rec sqrt_float_lo pp (s,n1,e1) =
  if s then
    failwith "sqrt_float_lo: negative argument"
  else
    if e1 land 1 = 1 then
      sqrt_float_lo pp (s, denormalize_nat (n1, 1), e1 - 1)
    else
      let p2 = pp * 2 in
      let f1' = denormalize_nat (n1, p2) in
      let f1 = Big_int.sqrt_big_int (big_int_of_num (dest_nat f1')) in
      let n, e' = lo_nat pp (mk_nat (num_of_big_int f1)) in
      let e = ((e1 + min_exp) lsr 1) + e' - pp in
	if e < 0 then
	  failwith "sqrt_float_lo: underflow"
	else
	  (s, n, e);;


let rec sqrt_float_hi pp (s,n1,e1) =
  if s then
    failwith "sqrt_float_hi: negative argument"
  else
    if e1 land 1 = 1 then
      sqrt_float_hi pp (s, denormalize_nat (n1, 1), e1 - 1)
    else
      let p2 = pp * 2 in
      let x = (big_int_of_num o dest_nat o denormalize_nat) (n1, p2) in
      let f1' = Big_int.sqrt_big_int x in
      let f1 = (mk_nat o num_of_big_int) f1' in
      let n, e' = 
	let ( * ) = Big_int.mult_big_int and
	    (==) = Big_int.eq_big_int in
	  hi_nat pp (if f1' * f1' == x then f1 else suc_nat f1) in
      let e = ((e1 + min_exp) lsr 1) + e' - pp in
	if e < 0 then
	  failwith "sqrt_float_hi: underflow"
	else
	  (s, n, e);;

end;;


(* Interval arithmetic *)


module type Informal_interval_sig =
  sig
    type interval
    val one_interval : interval
    val two_interval : interval
    val mk_interval : Informal_float.ifloat * Informal_float.ifloat -> 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;;



*)