(******************************************* ** RoundDownFloatInterval.sml ** sml ** ** Aleksandar Nanevski ** ** Implementation of interval arithmetic. ** Requires rounding down. *******************************************) structure RoundDownFloatInterval :> FLOAT_INTERVAL = struct exception Unordered structure R = Real (* assumes rounding to NegInf *) (* second argument is stored inverted, so that it *) (* gets rounded up *) (* a bound of the interval can be a NaN or an infinity *) datatype interval = Interval of real * real type t = interval val zero = Interval(0.0, Real.~ 0.0) val one = Interval(1.0, Real.~ 1.0) fun fromNumber x = Interval(x, Real.~ x) fun toNumber (Interval(x, y)) = 0.5 * (x - y) fun fromBounds {lo=x, hi=y} = if Real.<=(x, y) then Interval(x, Real.~ y) else raise Domain fun fromMidErr {mid=x, err=e} = Interval (x-e, ~x - e) fun toBounds (Interval(x, y)) = {lo=x, hi=Real.~ y} (* partial ordering (pairwise) *) fun op == (Interval(a1,a2), Interval(b1,b2)) = Real.==(a1, b1) andalso Real.==(a2, b2) fun op < (Interval(a1, a2), Interval(b1, b2)) = Real.<(a1, b1) andalso Real.>(a2, b2) fun op <= (Interval(a1, a2), Interval(b1, b2)) = Real.<= (a1, b1) andalso Real.>=(a2, b2) fun op >= (Interval(a1, a2), Interval(b1, b2)) = Real.>= (a1, b1) andalso Real.<=(a2, b2) fun op > (Interval(a1, a2), Interval(b1, b2)) = Real.>(a1, b1) andalso Real.<(a2, b2) fun op < (Interval(a1, a2), Interval(b1, b2)) = Real.<(a1, b1) andalso Real.>(a2, b2) fun op != (Interval(a1, a2), Interval(b1, b2)) = Real.!=(a1, b1) orelse Real.!=(a2, b2) fun inf(Interval(a1, a2), Interval(b1, b2)) = Interval(Real.min(a1, b1), Real.max(a2, b2)) fun sup(Interval(a1, a2), Interval(b1, b2)) = Interval(Real.max(a1, b1), Real.min(a2, b2)) (* raise Unordered if contains 0, but not [0,0], *) (* or cannot decide because of NaNs *) fun sign (a as Interval(a1, a2)) = if Real.>(a1, 0.0) then 1 else if Real.>(a2, 0.0) then ~1 else if (Real.==(a1, 0.0) andalso Real.==(a2, 0.0)) then 0 else raise Unordered fun compare (Interval(a1, a2), Interval(b1, b2)) = case Real.compareReal(a1, b1) of IEEEReal.LESS => if Real.>(a2, b2) then LESS else raise Unordered | IEEEReal.EQUAL => if Real.==(a2, b2) then EQUAL else raise Unordered | IEEEReal.GREATER => if Real.<(a2, b2) then GREATER else raise Unordered | _ => raise Unordered (* set related operations *) fun << (Interval(_, a2), Interval(b1, _)) = Real.<(Real.~ a2, b1) fun >> (Interval(a1, _), Interval(_, b2)) = Real.>(a1, Real.~ b2) (* seems to work correctly in all the cases *) (* including when one of the bounds is inf or NaN *) fun disjoint (Interval(a1, a2), Interval(b1, b2)) = Real.<(Real.~ a2, b1) orelse Real.<(Real.~ b2, a1) (* raise Unorder if disjoint or can't conclude *) (* because of NaNs *) fun meet (Interval(a1, a2), Interval(b1, b2)) = let val lo = Real.max(a1, b1) val hi = Real.max(a2, b2) in if Real.<=(lo, hi) then Interval(lo, hi) else raise Unordered end fun join (Interval(a1, a2), Interval(b1, b2)) = let val lo = Real.min(a1, b1) val hi = Real.min(a2, b2) in if Real.<=(lo, hi) then Interval(lo, hi) else raise Unordered end fun subset (Interval(a1, a2), Interval(b1, b2)) = Real.>=(a1, b1) andalso Real.>=(a2, b2) fun member (a, Interval(a1, a2)) = Real.<=(a1, a) andalso Real.<=(a, Real.~ a2) fun singleton (Interval(a1, a2)) = Real.==(a1, Real.~ a2) (* threaded continuous monotone functions *) val threadInf = inf val threadSup = sup fun threadAbs(a as Interval(a1, a2)) = if (Real.>(a1, 0.0)) then a else if (Real.>(a2, 0.0)) then Interval(a2, a1) else Interval(0.0, Real.max(Real.abs(a1), Real.abs(a2))) (* arithmetic operations *) fun op + (Interval(a1, a2), Interval(b1, b2)) = Interval(Real.+(a1, b1), Real.+(a2, b2)) fun op - (Interval(a1, a2), Interval(b1, b2)) = Interval(Real.+(a1, b2), Real.+(a2, b1)) fun op * (Interval(a1, a2), Interval(b1, b2)) = if Real.>=(a1, 0.0) then (* A >= [0,0] *) if Real.>=(b1, 0.0) then (* B >= [0,0] *) Interval(Real.*(a1, b1), Real.*(Real.~ a2, b2)) else if Real.>=(b2, 0.0) then (* B <= [0,0] *) Interval(Real.*(Real.~ a2, b1), Real.*(a1, b2)) else (* 0 \in B , or B = [NaN, NaN] *) let val a2' = Real.~ a2 in Interval(Real.*(a2', b1), Real.*(a2', b2)) end else if Real.>=(a2, 0.0) then (* A <= [0,0] *) if Real.>=(b1, 0.0) then (* B >= [0,0] *) Interval(Real.*(a1, Real.~ b2), Real.*(a2, b1)) else if Real.>=(b2, 0.0) then (* B <= [0,0] *) Interval(Real.*(a2, b2), Real.*(Real.~ a1, b1)) else (* 0 \in B, or B = [NaN, NaN] *) let val a1' = Real.~ a1 in Interval(Real.*(a1', b2), Real.*(a1',b1)) end else (* 0 \in A, or A = [NaN, NaN] *) if Real.>=(b1, 0.0) then (* B >= [0,0] *) let val b2' = Real.~ b2 in Interval(Real.*(a1, b2'), Real.*(a2, b2')) end else if Real.>=(b2, 0.0) then (* B <= [0,0] *) let val b1' = Real.~ b1 in Interval(Real.*(a2, b1'), Real.*(a1, b1')) end else (* 0 \in B, or B = [NaN, NaN] *) let val (a1', a2') = (Real.~ a1, Real.~ a2) val (l1, l2) = (Real.*(a1', b2), Real.*(a2', b1)) val (u1, u2) = (Real.*(a1', b1), Real.*(a2', b2)) in if Real.<(l1, l2) then if Real.<(u1, u2) then Interval(l1, u1) else Interval(l1, u2) else if Real.<(u1, u2) then Interval(l2, u1) else Interval(l2, u2) end fun op / (Interval(a1, a2), Interval(b1, b2)) = if Real.>=(a1, 0.0) then (* A >= [0,0] *) if Real.>(b1, 0.0) then (* B > [0,0] *) Interval(Real./(a1, Real.~ b2), Real./(a2, b1)) else if Real.>(b2, 0.0) then (* B < [0,0] *) Interval(Real./(a2, b2), Real./(a1,Real.~ b1)) else (* 0 \in B, or B = [NaN, NaN] *) Interval(Real.negInf, Real.negInf) (* raise Div *) else if Real.>=(a2, 0.0) then (* A <= [0,0] *) if Real.>(b1, 0.0) then (* B > [0,0] *) Interval(Real./(a1, b1), Real./(a2, Real.~ b2)) else if Real.>(b2, 0.0) then (* B < [0,0] *) Interval(Real./(Real.~ a2, b1), Real./(a1, b2)) else (* 0 \in B *) Interval(Real.negInf, Real.negInf) (* raise Div *) else (* 0 \in A *) if Real.>(b1, 0.0) then (* B > [0,0] *) Interval(Real./(a1, b1), Real./(a2, b1)) else if Real.>(b2, 0.0) then (* B < [0,0] *) Interval(Real./(a2, b2), Real./(a1, b2)) else (* 0 \in B *) Interval(Real.negInf, Real.negInf) (* raise Div *) fun scale (a, Interval(b1, b2)) = if Real.>=(a, 0.0) then Interval(Real.*(a, b1), Real.*(a, b2)) else let val a' = Real.~ a in Interval(Real.*(a', b2), Real.*(a', b1)) end fun op ~ (Interval(a1, a2)) = Interval(a2, a1) fun inv a = one / a (* structure-specific functions *) fun width (Interval(a1, a2)) = let open Real in {lo = ~a2 - a1, hi = ~ (a1 + a2)} end fun distance(Interval(a1, a2), Interval(b1, b2)) = let open Real val (l1, h1) = if a1 >= b1 then (a1 - b1, b1 - a1) else (b1 - a1, a1 - b1) val (l2, h2) = if a2 >= b2 then (a2 - b2, b2 - a2) else (b2 - a2, a2 - b2) in {lo = max(l1, l2), hi = ~ (min(h1, h2))} end fun abs(Interval(a1, a2)) = Real.max(Real.abs a1, Real.abs a2) fun isNan(Interval(a1, a2)) = (Real.isNan a1) orelse (Real.isNan a2) fun isFinite(Interval(a1, a2)) = (Real.isFinite a1) andalso (Real.isFinite a2) fun toString (Interval(a1, a2)) = "["^ (Real.toString a1) ^", "^ (Real.toString (Real.~ a2)) ^ "]" end