(* This module is a translation in ocaml of Hans's Boehm Java library CR *)

(*i
// Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED 
// 
// Permission is granted free of charge to copy, modify, use and distribute
// this software  provided you include the entirety of this notice in all
// copies made.
// 
// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY
// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT
// FOR A PARTICULAR PURPOSE OR NON-INFRINGING.   SGI ASSUMES NO RISK AS TO THE
// QUALITY AND PERFORMANCE OF THE SOFTWARE.   SHOULD THE SOFTWARE PROVE
// DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY
// SERVICING, REPAIR OR CORRECTION.  THIS DISCLAIMER OF WARRANTY CONSTITUTES
// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS
// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
// 
// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING,
// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR
// OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL,
// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE
// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK
// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL
// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF
// THE POSSIBILITY OF SUCH DAMAGES.  THIS LIMITATION OF LIABILITY SHALL NOT
// APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE
// LAW PROHIBITS SUCH LIMITATION.  SOME JURISDICTIONS DO NOT ALLOW THE
// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT
// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU.
// 
// These license terms shall be governed by and construed in accordance with
// the laws of the United States and the State of California as applied to
// agreements entered into and to be performed entirely within California
// between California residents.  Any litigation relating to these license
// terms shall be subject to the exclusive jurisdiction of the Federal Courts
// of the Northern District of California (or, absent subject matter
// jurisdiction in such courts, the courts of the State of California), with
// venue lying exclusively in Santa Clara County, California. 

// Copyright (c) 2001-2004, Hewlett-Packard Development Company, L.P. 
// 
// Permission is granted free of charge to copy, modify, use and distribute
// this software  provided you include the entirety of this notice in all
// copies made.
// 
// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY
// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT
// FOR A PARTICULAR PURPOSE OR NON-INFRINGING.   HEWLETT-PACKARD ASSUMES
// NO RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE.
// SHOULD THE SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, 
// HEWLETT-PACKARD ASSUMES NO COST OR LIABILITY FOR ANY
// SERVICING, REPAIR OR CORRECTION.  THIS DISCLAIMER OF WARRANTY CONSTITUTES
// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS
// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER.
// 
// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING,
// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR
// OTHERWISE, SHALL HEWLETT-PACKARD BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL,
// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE
// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK
// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL
// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF HEWLETT-PACKARD SHALL
// HAVE BEEN INFORMED OF THE POSSIBILITY OF SUCH DAMAGES.
// THIS LIMITATION OF LIABILITY SHALL NOT APPLY TO LIABILITY RESULTING
// FROM HEWLETT-PACKARD's NEGLIGENCE TO THE EXTENT APPLICABLE
// LAW PROHIBITS SUCH LIMITATION.  SOME JURISDICTIONS DO NOT ALLOW THE
// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT
// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU.
// 
i*)

open Gmp

type t = {
  mutable cache : (int * Z.t) option;
  approximate : int -> Z.t 
}

let create f = { cache = None; approximate = f }

let bound_log2 n =
  let np1 = float (abs n + 1) in
  truncate (ceil (log np1 /. log 2.0))

let z_zero = Z.from_int 0
let z_one = Z.from_int 1
let z_mone = Z.from_int (-1)
let z_two = Z.from_int 2
let z_three = Z.from_int 3
let z_four = Z.from_int 4
let z_six = Z.from_int 6

let z_shift_left = Z.mul2exp
let z_shift_right = Z.fdiv_q_2exp

let shift k n =
  if n == 0 then k
  else if n < 0 then z_shift_right k (-n)
  else z_shift_left k n

let scale k n =
  if n >= 0 then
    z_shift_left k n
  else
    let adj_k = Z.add_ui (shift k (n+1)) 1 in
    z_shift_right adj_k 1

exception PrecisionOverflow

let check_prec n = 
  let high = n lsr 28 in
  let high_shifted = n lsr 29 in
  (* TODO *)
  ()

let approx x p =
  check_prec p;
  match x.cache with
    | Some (min_p, ma) when p >= min_p -> 
	scale ma (min_p - p)
    | _ -> 
	let r = x.approximate p in
	x.cache <- Some (p, r);
	r

(* ceil(log2(if z < 0 then -z else z+1)) *)
let z_bit_length z =
  let rec loop k two_k =
    if Z.cmp z two_k <= 0 then k else loop (succ k) (Z.mul2exp two_k 1)
  in
  let z = if Z.sgn z < 0 then Z.neg z else Z.add_ui z 1 in
  if Z.cmp z z_one = 0 then 0 else loop 1 z_two

let known_msd x = match x.cache with
  | Some (mp, ma) ->
      let length = 
	if Z.sgn ma >= 0 then z_bit_length ma else z_bit_length (Z.neg ma)
      in
      mp + length - 1
  | None -> 
      assert false

let msd_n x n =
  let close_to_0 = match x.cache with
    | None -> true
    | Some (_, ma) -> Z.cmp ma z_one <= 0 && Z.cmp ma z_mone >= 0
  in
  if close_to_0 then begin
    ignore (approx x (n-1));
    match x.cache with
      | Some (_, ma) when Z.cmp (Z.abs ma) z_one <= 0 -> min_int
      | Some _ -> known_msd x
      | None -> assert false
  end else
    known_msd x

let iter_msd x n =
  let rec iter prec =
    if prec > n + 30 then
      let msd = msd_n x prec in
      if msd != min_int then msd
      else begin
	check_prec prec;
	iter ((prec * 3)/2 - 16)
      end
    else
      msd_n x n
  in
  iter 0

let msd x = iter_msd x min_int

let shifted x n = create (fun p -> approx x (p - n))
let shift_left x n = shifted x n
let shift_right x n = shifted x (-n)

let of_z z = create (fun p -> scale z (-p))

let of_int n = of_z (Z.from_int n)

let of_int64 n = let z = Z.from_string (Int64.to_string n) in of_z z

let zero = of_int 0
let one = of_int 1
let two = of_int 2
let four = of_int 4

(* Operations *)

let add x y =
  create (fun p -> scale (Z.add (approx x (p-2)) (approx y (p-2))) (-2))

let (+!) = add

let neg x = create (fun p -> Z.neg (approx x p))

let sub x y = add x (neg y)

let (-!) = sub

exception Zero

let mul x y =
  create 
    (fun p ->
       let half_prec = (p asr 1) - 1 in
       try
	 let x,y,msd_x =
	   let msd_x = msd_n x half_prec in
	   if msd_x == min_int then
	     let msd_y = msd_n y half_prec in
	     if msd_y == min_int then raise Zero else y,x,msd_y
	   else
	     x,y,msd_x
	 in
	 let prec2 = p - msd_x - 3 in
	 let appry = approx y prec2 in
	 if Z.sgn appry == 0 then raise Zero;
	 let msd_y = known_msd y in
	 let prec1 = p - msd_y - 3 in
	 let apprx = approx x prec1 in
	 let scale_digits = prec1 + prec2 - p in
	 scale (Z.mul apprx appry) scale_digits
       with Zero ->
	 z_zero)

let ( *! ) = mul

let inv x =
  create
    (fun p ->
       let msd = msd x in
       let inv_msd = 1 - msd in
       let digits_needed = inv_msd - p + 3 in
       let prec_needed = msd - digits_needed in
       let log_scale_factor = -p - prec_needed in
       if log_scale_factor < 0 then
	 z_zero
       else
	 let dividend = z_shift_left z_one log_scale_factor in
	 let scaled_divisor = approx x prec_needed in
	 let abs_scaled_divisor = Z.abs scaled_divisor in
	 let adj_dividend = 
	   Z.add dividend (z_shift_right abs_scaled_divisor 1) 
	 in
	 let result = Z.fdiv_q adj_dividend abs_scaled_divisor in
	 if Z.sgn scaled_divisor < 0 then Z.neg result else result)

let div x y = mul x (inv y)

let (/!) = div

let sgn_n x a = 
  let quick_try = match x.cache with
    | Some (_, ma) -> Z.sgn ma
    | None -> 0
  in
  if quick_try != 0 then 
    quick_try
  else
    let needed_prec = a - 1 in
    let appr = approx x needed_prec in
    Z.sgn appr

let sgn x =
  let rec loop a =
    check_prec a;
    let result = sgn_n x a in
    if result != 0 then result else loop (2 * a)
  in
  loop (-20)

let select s x y =
  let selector_sign = ref (Z.sgn (approx s (-20))) in
  create
    (fun p ->
       if !selector_sign < 0 then
	 approx x p
       else if !selector_sign > 0 then
	 approx y p
       else
	 let x_appr = approx x (p-1) in
	 let y_appr = approx y (p-1) in
	 let diff = Z.abs (Z.sub x_appr y_appr) in
	 if Z.cmp diff z_one <= 0 then
	   scale x_appr (-1)
	 else if sgn s < 0 then begin
	   selector_sign := -1;
	   scale x_appr (-1)
	 end else begin
	   selector_sign := 1;
	   scale y_appr (-1)
	 end)

let abs x = select x (neg x) x

let max x y = select (sub x y) y x
let min x y = select (sub x y) x y

let compare_a x y a =
  let needed_prec = a - 1 in
  let x_appr = approx x needed_prec in
  let y_appr = approx y needed_prec in
  let comp1 = Z.cmp x_appr (Z.add_ui y_appr 1) in
  if comp1 > 0 then 
    1
  else 
    let comp2 = Z.cmp x_appr (Z.add_ui y_appr (-1)) in
    if comp2 < 0 then -1 else 0

let compare x y =
  let rec loop a =
    check_prec a;
    let r = compare_a x y a in
    if r <> 0 then r else loop (2 * a)
  in
  loop (-20)

let rec pow_int x n =
  if n == 0 then
    one
  else if n < 0 then
    inv (pow_int x (-n))
  else 
    let y = pow_int (mul x x) (n / 2) in
    if n mod 2 == 0 then y else mul y x

let prescaled_exp x =
  create
    (fun p ->
       if p >= 1 then
	 z_zero
       else
	 let iterations_needed = -p/2 + 2 in
	 let calc_precision = p - bound_log2(2*iterations_needed) - 4 in
	 let op_prec = p - 3 in
	 let op_appr = approx x op_prec in
	 let scaled_1 = z_shift_left z_one (-calc_precision) in
	 let current_term = ref scaled_1 in
	 let current_sum = ref scaled_1 in
	 let n = ref 0 in
	 let max_trunc_error = z_shift_left z_one (p - 4 - calc_precision) in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   incr n;
	   current_term := scale (Z.mul !current_term op_appr) op_prec;
	   current_term := Z.fdiv_q_ui !current_term !n;
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p))

let rec exp x =
  let low_prec = -10 in
  let rough_appr = approx x low_prec in
  if Z.sgn rough_appr < 0 then
    inv (exp (neg x))
  else if Z.cmp rough_appr z_two > 0 then
    let square_root = exp (shift_right x 1) in
    mul square_root square_root
  else
    prescaled_exp x

let e = prescaled_exp one

let sqrt x =
  let fp_prec = 50 in
  let fp_op_prec = 60 in
  let rec sqrt_rec p =
    let max_prec_needed = 2*p - 1 in
    let msd = msd_n x max_prec_needed in
    if msd <= max_prec_needed then
      z_zero
    else
      let result_msd = msd / 2 in
      let result_digits = result_msd - p in
      if result_digits > fp_prec then
	let appr_digits = result_digits/2 + 6 in
	let appr_prec = result_msd - appr_digits in
	let last_appr = sqrt_rec appr_prec in
	let prod_prec = 2 * appr_prec in
	let op_appr = approx x prod_prec in
	let prod_prec_scaled_numerator = 
	  Z.add (Z.mul last_appr last_appr) op_appr
	in
	let scaled_numerator = 
	  scale prod_prec_scaled_numerator (appr_prec - p)
	in
	assert (Z.cmp last_appr z_zero != 0);
	let shifted_result = Z.fdiv_q scaled_numerator last_appr in
	z_shift_right (Z.add shifted_result z_one) 1
      else begin
	let op_prec = (msd - fp_op_prec) land (lnot 1) in
	let working_prec = op_prec - fp_op_prec in
	let scaled_bi_appr = z_shift_left (approx x op_prec) fp_op_prec in
	let scaled_appr = Z.float_from scaled_bi_appr in
	if scaled_appr < 0.0 then invalid_arg "Cr.sqrt";
	let scaled_fp_sqrt = sqrt scaled_appr in
	let scaled_sqrt = 
	  Z.from_string (Int64.to_string (Int64.of_float scaled_fp_sqrt)) 
	in
	let shift_count = working_prec/2 - p in
	shift scaled_sqrt shift_count
      end
  in
  create sqrt_rec


let prescaled_ln x = 
  create
    (fun p ->
       if p >= 0 then
	 z_zero
       else
	 let iterations_needed = -p in
	 let calc_precision = p - bound_log2(2 * iterations_needed) - 4 in
	 let op_prec = p - 3 in
	 let op_appr = approx x op_prec in
	 let scaled_1 = z_shift_left z_one (-calc_precision) in
	 let x_nth = ref (scale op_appr (op_prec - calc_precision)) in
	 let current_term = ref !x_nth in
	 let current_sum = ref !current_term in
	 let n = ref 1 in
	 let current_sign = ref 1 in
	 let max_trunc_error = z_shift_left z_one (p - 4 - calc_precision) in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   incr n;
	   current_sign := - !current_sign;
	   x_nth := scale (Z.mul !x_nth op_appr) op_prec;
	   current_term := Z.fdiv_q !x_nth (Z.of_int (!n * !current_sign));
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p))

let simple_ln x = prescaled_ln (sub x one)

(* ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80) *)
let ten_ninths = of_int 10 /! of_int 9
let ln2_1 = of_int 7 *! simple_ln ten_ninths
let twentyfive_twentyfourths = of_int 25 /! of_int 24
let ln2_2 = of_int 2 *! simple_ln twentyfive_twentyfourths
let eightyone_eightyeths = of_int 81 /! of_int 80
let ln2_3 = of_int 3 *! simple_ln eightyone_eightyeths
let ln2 = ln2_1 -! ln2_2 +! ln2_3

let low_ln_limit = Z.of_int 8
let high_ln_limit = Z.of_int (16 + 8)
let scaled_4 = Z.of_int (4 * 16)

let rec ln x =
  let low_prec = -4 in
  let rough_appr = approx x low_prec in
  if Z.cmp rough_appr z_zero < 0 then invalid_arg "Cr.ln";
  if Z.cmp rough_appr low_ln_limit <= 0 then
    neg (ln (inv x))
  else if Z.cmp rough_appr high_ln_limit >= 0 then
    if Z.cmp rough_appr scaled_4 <= 0 then
      let quarter = ln (sqrt (sqrt x)) in
      shift_left quarter 2
    else
      let extra_bits = z_bit_length rough_appr - 3 in
      let scaled_result = ln (shift_right x extra_bits) in
      add scaled_result (mul (of_int extra_bits) ln2)
  else
    simple_ln x

let log ~base:x y = ln y /! ln x

let pow x y = exp (y *! ln x)

let root n x = pow x (inv (of_int n))

let arctan_reciproqual op =
  create
    (fun p ->
       if p >= 1 then 
	 z_zero
       else
	 let iterations_needed = -p/2 + 2 in
	 let calc_precision = p - bound_log2(2 * iterations_needed) - 2 in
	 let scaled_1 = z_shift_left z_one (-calc_precision) in
	 let big_op = Z.of_int op in
	 let big_op_squared = Z.of_int (op*op) in
	 let op_inverse = Z.fdiv_q scaled_1 big_op in
	 let current_power = ref op_inverse in
	 let current_term = ref op_inverse in
	 let current_sum = ref op_inverse in
	 let current_sign = ref 1 in
	 let n = ref 1 in
	 let max_trunc_error = z_shift_left z_one (p - 2 - calc_precision) in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   n := !n + 2;
	   current_power := Z.fdiv_q !current_power big_op_squared;
	   current_sign := - !current_sign;
	   current_term := 
	     Z.fdiv_q !current_power (Z.of_int (!current_sign * !n));
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p))

let pi = 
     (of_int 48 *! arctan_reciproqual 18)
  +! (of_int 32 *! arctan_reciproqual 57)
  -! (of_int 20 *! arctan_reciproqual 239)

let half_pi = shift_right pi 1

let prescaled_cos x =
  create
    (fun p ->
       if p >= 1 then
	 z_zero
       else begin
	 let iterations_needed = -p/2 + 4 in
	 let calc_precision = p - bound_log2(2 * iterations_needed) - 4 in
	 let op_prec = p - 2 in
	 let op_appr = approx x op_prec in
	 let current_term = ref (z_shift_left z_one (-calc_precision)) in
	 let n = ref 0 in
	 let max_trunc_error = z_shift_left z_one (p - 4 - calc_precision) in
	 let current_sum = ref !current_term in
	 while Z.cmp (Z.abs !current_term) max_trunc_error >= 0 do
	   n := !n + 2;
	   current_term := scale (Z.mul !current_term op_appr) op_prec;
	   current_term := scale (Z.mul !current_term op_appr) op_prec;
	   let divisor = Z.mul (Z.of_int (- !n)) (Z.of_int (!n-1)) in
	   current_term := Z.fdiv_q !current_term divisor;
	   current_sum := Z.add !current_sum !current_term
	 done;
	 scale !current_sum (calc_precision - p)
       end)

let rec cos x =
  let rough_appr = approx x (-1) in
  let abs_rough_appr = Z.abs rough_appr in
  if Z.cmp abs_rough_appr z_six >= 0 then
    let multiplier = Z.fdiv_q_ui rough_appr 6 in
    let adjustment = mul pi (of_z multiplier) in
    if Z.sgn (Z.band multiplier z_one) != 0 then
      neg (cos (x -! adjustment))
    else
      cos (x -! adjustment)
  else if Z.cmp abs_rough_appr z_two >= 0 then
    let cos_half = cos (shift_right x 1) in
    (shift_left (cos_half *! cos_half) 1) -! one
  else
    prescaled_cos x

let sin x = cos (half_pi -! x)

let tan x = sin x /! cos x

(*s Hyperbolic functions. *)

let sinh x = let expx = exp x in (expx -! inv expx) /! two

let cosh x = let expx = exp x in (expx +! inv expx) /! two

let tanh x = 
  let expx = exp x in
  let exp_minus_x = inv expx in
  (expx -! exp_minus_x) /! (expx +! exp_minus_x)

let arcsinh x = ln (x +! sqrt (x *! x +! one))
let arccosh x = ln (x +! sqrt (x *! x -! one))
let arctanh x = ln ((one +! x) /! (one -! x)) /! two

let of_float n =
  begin match classify_float n with
    | FP_nan | FP_infinite -> invalid_arg "Cr.of_float"
    | _ -> ()
  end;
  let negative = n < 0.0 in
  let bits = Int64.bits_of_float (abs_float n) in
  let mantissa = Int64.logand bits 0xfffffffffffffL in
  let biased_exp = Int64.to_int (Int64.shift_right_logical bits 52) in
  let exp = biased_exp - 1075 in
  let mantissa = 
    if biased_exp != 0 then
      Int64.add mantissa (Int64.shift_left Int64.one 52)
    else
      Int64.shift_left mantissa 1
  in
  let result = shift_left (of_int64 mantissa) exp in
  if negative then neg result else result

let to_string ?(radix=10) x n =
  if n < 0 then invalid_arg "Cr.to_string";
  let scaled_x =
    if radix == 16 then 
      shift_left x (4 * n)
    else
      let scale_factor = Z.ui_pow_ui radix n in
      mul x (of_z scale_factor)
  in
  let scaled_int = approx scaled_x 0 in
  let scaled_string = Z.to_string_base ~base:radix (Z.abs scaled_int) in
  let result = 
    if n == 0 then
      scaled_string
    else
      let len = String.length scaled_string in
      let len, scaled_string =
	if len <= n then
	  n + 1, String.make (n + 1 - len) '0' ^ scaled_string
	else
	  len, scaled_string
      in
      let whole = String.sub scaled_string 0 (len - n) in
      let fraction = String.sub scaled_string (len - n) n in
      whole ^ "." ^ fraction
  in
  if Z.sgn scaled_int < 0 then "-" ^ result else result

let of_string ?(radix=10) s =
  try
    begin
      try
	let n = String.length s in
	let p = String.index s '.' in
	let dec = n - p - 1 in
	let s' = (String.sub s 0 p) ^ (String.sub s (p + 1) dec) in
	div (of_z (Z.from_string_base radix s')) (of_z (Z.pow_ui_ui radix dec))
      with Not_found -> 
	of_z (Z.from_string_base radix s)
    end
  with Invalid_argument _ -> invalid_arg "Cr.of_string"

let to_q x n =
  let xn = approx x (-n) in
  Q.div (Q.from_z xn) (Q.from_z (z_shift_left z_one n))

let to_float x n = Q.float_from (to_q x n)

(* Inverse of a monotone function (see file UnaryCRFunction.java) *)

let sloppy_compare x y =
  let difference = Z.sub x y in
  if Z.cmp difference z_one > 0 then 1 
  else if Z.cmp difference z_mone < 0 then -1 else 0

let trace = false
let printf = Format.printf

let inverse_monotone f ~low ~high =
  let f_low = f low in
  let f_high = f high in
  let f,negated,f_low,f_high =
    if compare f_low f_high > 0 then
      (fun x -> neg (f x)), true, neg f_low, neg f_high
    else
      f, false, f_low, f_high
  in
  let max_msd = msd (max (abs low) (abs high)) in
  let max_arg_prec = msd (high -! low) - 4 in
  let deriv_msd = msd ((f_high -! f_low) /! (high -! low)) in
  fun x ->
    let arg = if negated then neg x else x in
    let rec r = { 
      cache = None;
      approximate = fun p ->
	let digits_needed = max_msd - p in
	if digits_needed < 0 then 
	  z_zero
	else
	  let extra_arg_prec = 4 in
	  let working_arg_prec = 
	    Stdlib.min (p - extra_arg_prec) max_arg_prec
	  in
	  let working_eval_prec = ref (working_arg_prec + deriv_msd - 20) in
	  let low_appr = Z.add_ui (approx low working_arg_prec) 1 in
	  let high_appr = Z.sub_ui (approx high working_arg_prec) 1 in
	  let arg_appr = ref (approx arg !working_eval_prec) in
	  let have_good_appr = match r.cache with
	    | Some (min_prec, _) -> min_prec < max_msd
	    | None -> false
	  in
	  let small_steps = ref 0 in
	  let l,f_l,h,f_h,at_left,at_right =
	    if digits_needed < 30 && not have_good_appr then begin
	      if trace then printf "Setting interval to entire domain@.";
	      small_steps := 2;
	      low_appr, approx f_low !working_eval_prec,
	      high_appr, approx f_high !working_eval_prec,
	      true, true
	    end else
	      let rough_prec = p + digits_needed/2 in
	      let rough_prec = match r.cache with
		| Some (min_prec, _) when
                  digits_needed < 30 || min_prec<p+3*digits_needed/4 ->
		    min_prec
		| _ -> rough_prec
	      in
	      let rough_appr = approx r rough_prec in
	      if trace then begin
		printf "Setting interval based on prev. appr@.";
		printf "prev. prec = %d appr = %a@." 
		  rough_prec Z.print rough_appr
	      end;
	      let h = z_shift_left (Z.add_ui rough_appr 1) 
		(rough_prec - working_arg_prec)
	      in
	      let l = z_shift_left (Z.sub_ui rough_appr 1)
		(rough_prec - working_arg_prec)
	      in
	      let h,f_h,at_right = 
		if Z.cmp h high_appr > 0 then
		  high_appr, approx f_high !working_eval_prec, true
		else
		  let h_cr = shift_left (of_z h) working_arg_prec in
		  let f_h = approx (f h_cr) !working_eval_prec in
		  h, f_h, false
	      in
	      let l,f_l,at_left =
		if Z.cmp l low_appr < 0 then
		  low_appr, approx f_low !working_eval_prec, true
		else
		  let l_cr = shift_left (of_z l) working_arg_prec in
		  let f_l = approx (f l_cr) !working_eval_prec in
		  l, f_l, false
	      in
	      l,f_l,h,f_h,at_left,at_right
	  in
	  let l,f_l,h,f_h,at_left,at_right =
	    ref l, ref f_l, ref h, ref f_h, ref at_left, ref at_right
	  in
	  let difference = ref (Z.sub !h !l) in
	  let rec loop i =
	    if trace then begin
	      printf "***Iteration: %d@." i;
	      printf "Arg prec = %d eval prec = %d arg appr. = %a@."
		working_arg_prec !working_eval_prec Z.print !arg_appr;
	      printf "l = %a h = %a@." Z.print !l Z.print !h;
	      printf "f(l) = %a; f(h) = %a@." Z.print !f_l Z.print !f_h
	    end;
	    if Z.cmp !difference z_six < 0 then
	      scale !h (-extra_arg_prec)
	    else
	      let f_difference = Z.sub !f_h !f_l in
	      let guess = 
		if !small_steps >= 2 || Z.sgn f_difference = 0 then
		  z_shift_right (Z.add !l !h) 1 
		else
		  let arg_difference = Z.sub !arg_appr !f_l in
		  let t = Z.mul arg_difference !difference in
		  let adj = Z.fdiv_q t f_difference in
		  let adj =
		    if Z.cmp adj (z_shift_right !difference 2) < 0 then
		      z_shift_left adj 1
		    else if Z.cmp adj 
		      (z_shift_right (Z.mul_ui !difference 3) 2) > 0 then
			Z.sub !difference 
			  (z_shift_left (Z.sub !difference adj) 1)
		    else
		      adj
		  in
		  let adj = if Z.sgn adj <= 0 then z_two else adj in
		  let adj = 
		    if Z.cmp adj !difference >= 0 then
		      Z.sub_ui !difference 2 
		    else 
		      adj
		  in
		  if Z.sgn adj <= 0 then Z.add_ui !l 2 else Z.add !l adj 
	      in
	      let guess = ref guess in
	      let tweak = ref z_two in
	      let rec loop2 adj_prec =
		let guess_cr = shift_left (of_z !guess) working_arg_prec in
		if trace then begin
		  printf "Evaluating at %s with precision %d@."
		    (to_string guess_cr 10) !working_eval_prec;
		end;
		let f_guess_cr = f guess_cr in
		if trace then begin
		  printf "fn value = %s@." (to_string f_guess_cr 10)
		end;
		let f_guess = approx f_guess_cr !working_eval_prec in
		let outcome = sloppy_compare f_guess !arg_appr in
		if outcome <> 0 then begin (* break *) 
		  if outcome > 0 then begin
		    h := !guess; f_h := f_guess; at_right := false
		  end else begin
		    l := !guess; f_l := f_guess; at_left := false
		  end;
		  let new_difference = Z.sub !h !l in
		  if Z.cmp new_difference (z_shift_right !difference 1) >= 0
		  then incr small_steps
		  else small_steps := 0;
 		  difference := new_difference;
		  loop (i+1)
		end else begin
		  if adj_prec then begin
		    let adjustment = 
		      if deriv_msd > 0 then -20 else deriv_msd - 20
		    in
		    let l_cr = shift_left (of_z !l) working_arg_prec in
		    let h_cr = shift_left (of_z !h) working_arg_prec in
		    working_eval_prec := !working_eval_prec + adjustment;
		    if trace then begin
		      printf "New eval prec = %d %s%s@." !working_eval_prec
			(if !at_left then "(at left)" else "")
			(if !at_right then "(at right)" else "")
		    end;
		    if !at_left then
		      f_l := approx f_low !working_eval_prec
		    else
		      f_l := approx (f l_cr) !working_eval_prec;
		    if !at_right then
		      f_h := approx f_high !working_eval_prec
		    else
		      f_h := approx (f h_cr) !working_eval_prec;
		    arg_appr := approx arg !working_eval_prec;
		  end else begin
		    if trace then printf "tweaking guess@.";
		    let new_guess = Z.add !guess !tweak in
		    if Z.cmp new_guess !h >= 0 then
		      guess := Z.sub !guess !tweak
		    else
		      guess := new_guess;
		    tweak := Z.neg !tweak
		  end;
		  loop2 (not adj_prec)
		end
	      in
	      loop2 false
	  in
	  loop 0
    }
    in
    r

(* Application to the inverse trigonometric functions *)

let arcsin = inverse_monotone sin ~low:(neg half_pi) ~high:half_pi

let arccos x = half_pi -! arcsin x

(* uses the identity [(sin x)^2 = (tan x)^2/(1 + (tan x)^2)] *)
let arctan x =
  let x2 = x *! x in
  let abs_sin_atan = sqrt (x2 /! (one +! x2)) in
  let sin_atan = select x (neg abs_sin_atan) abs_sin_atan in
  arcsin sin_atan

(*s Format pretty-printer. *)

let print_precision = ref 10
let set_print_precision = (:=) print_precision
let print fmt x = Format.fprintf fmt "%s" (to_string x !print_precision)

(* Infix notations *)

module Infixes = struct
  let (+!) = add
  let (-!) = sub
  let ( *! ) = mul
  let (/!) = div
end