(*************************************
 ** Ulps.sml
 ** sml
 **
 ** Aleksandar Nanevski
 **
 ** Symbolic operations over the 
 ** machine epsilon
 *************************************)

structure Ulps :
sig
    include RATIONAL
    val E : t
end
=
struct
    open NormalizedRational

    val Izero = IntInf.fromInt 0
    val Ione = IntInf.fromInt 1
    val Itwo = IntInf.fromInt 2

    val e = IntInf.pow (Itwo, 53)
    val E = one / (fromIntInf e) 


    (* print the errors in Shewchuk style *)
    (* as sums of machine epsilons *)

    fun approx2 p =
	let val {num=a, denom=b} = toFrac p
	    fun approx' (a, b, exp) = 
		let open IntInf
		in
		    if a = Izero then (a, a, exp)
		    else
			if a >= b then 
			    let val (q, r) = quotrem (a, b)
				val (r, s) = quotrem (r * e, b)
			    in
				if s = Izero then (q, r, exp)
				else (q, r + Ione, exp)
			    end
			else approx' (a * e, b, Int.+(exp, 1))
		end
	in
	    approx' (a, b, 0)
	end


    fun toString p = 
	let val (q, r, exp) = approx2 p
	    val qq = NumberTheory.approx IEEEReal.TO_POSINF q
	    val rr = NumberTheory.approx IEEEReal.TO_POSINF r
	    fun multEs 0 = ""
	      | multEs i = "* E " ^ (multEs (Int.-(i, 1)))
	in
	    "("^(Real.toString qq)^" + "^(Real.toString rr)^"*E) "^(multEs exp)
	end

 
(* Old code. May still be useful for debugging  *)

(*
    fun toString p = 
	let val (q, r, exp) = approx2 p
	    open IntInf
	in
	    (toString q)^"E"^(Int.toString exp)^" + "^
	    (toString r)^"E"^(Int.toString (Int.+(exp, 1)))
	end



    val N = 5
    fun approxN p =
	let val {num=a, denom=b} = toFrac p
	    fun approx' (a, b, exp, acc) =
		if (exp = N) orelse (a = Izero) then acc else
		    let open IntInf
		    in
			if a >= b then 
			    let val (q, r) = quotrem (a, b)
			    in
				if r = Izero then {coef=q, exp=exp}::acc
				else approx' (r * e, b, Int.+(exp, 1), {coef=q, exp=exp}::acc)
			    end
			else approx' (a * e, b, Int.+(exp, 1), acc)
		    end
	in
	    IntInfPoly.fromCoefExpList (approx' (a, b, 0, nil))
	end


    fun toString p =
	IntInfPoly.toString (approxN p)
*)

end