(******************************************** ** NormalizedRational.sml ** sml ** ** Aleksandar Nanevski ** ** Implementation of infinite precision ** rational numbers, with normalization ** performed at every step. ********************************************) structure NormalizedRational :> RATIONAL = struct structure I = IntInf structure N = NumberTheory (* assumes denom > 0 *) type rat = {num: I.int, denom: I.int} type t = rat fun fromIntInf x = {num = x, denom = N.one} fun fromInt x = {num = I.fromInt x, denom = N.one} val zero = fromInt 0 val one = fromInt 1 fun normalize (x as {num=p, denom=q}) = let val d = N.gcd(p, q) in if d = N.one then x else {num = I.quot(p, d), denom = I.quot(q, d)} end fun fromFrac {num=x, denom=y} = if y = N.one then {num=x, denom=N.one} else normalize {num=x, denom=y} fun fromIntFrac {num=x, denom=y} = if y = 1 then {num=I.fromInt x, denom=N.one} else normalize {num=I.fromInt x, denom=I.fromInt y} fun fromManExp {man=man, exp=exp} = if Int.>(exp, 0) then {num = I.<<(man, Word.fromInt exp), denom = N.one} else normalize {num = man, denom = I.<<(N.one, Word.fromInt (Int.~ exp))} val bound = Real.precision + 1 val p960 = Real.fromManExp{man=1.0, exp=960} val q1020 = Real.fromManExp{man=1.0, exp= ~1020} val halfMax = Real.*(0.5, Real.maxFinite) val doubleMin = Real.*(2.0, Real.minPos) (* return the tightest possible enclosing interval of *) (* the fraction. the interval is computed as if the *) (* lower bound is obtained from the fraction by rounding *) (* down and the upper bound by rounding up *) fun intervalApprox ({num=a, denom=b}) = let val s = I.sign a in if s = 0 then {lo=0.0, hi=0.0} else let val a = I.abs a val (al, bl) = (I.log2 a, I.log2 b) val d = al - bl val (a', b') = if d > bound then (a, I.<<(b, Word.fromInt (d - bound))) else (I.<<(a, Word.fromInt(bound - d)), b) val (q, r) = I.quotrem(a', b') val (q', exp) = if (r = N.zero) then (q, d - bound) else (I.+(N.one, I.<<(q, 0w30)), (d - bound) - 30) val {lo=lo, hi=hi} = N.intervalApprox q' in if exp < ~1170 then if (s = 1) then {lo=0.0, hi=Real.minPos} else {lo=Real.~ Real.minPos, hi=0.0} else if exp < ~1020 then let open Real fun roundup (lo, hi, ~1) = if hi < doubleMin then {lo=0.0, hi=Real.minPos} else if lo < doubleMin then {lo=0.5*lo, hi=Real.minPos} else {lo=0.5*lo, hi=0.5*hi} | roundup (lo, hi, i) = if lo < doubleMin then {lo=0.0, hi=Real.minPos} else roundup(0.5*lo, 0.5*hi, Int.+(i, 1)) val (r as {lo=lo, hi=hi}) = roundup (lo*q1020, hi*q1020, Int.+(1020, exp)) in if (s = 1) then r else {lo=Real.~ hi, hi=Real.~ lo} end else if exp <= 960 then let open Real val p = fromManExp{man=1.0, exp=exp} in if (s = 1) then {lo=lo*p, hi=hi*p} else {lo=p * Real.~ hi, hi=p * Real.~ lo} end else if exp > 1023 then if (s = 1) then {lo=Real.maxFinite, hi=Real.posInf} else {lo=Real.negInf, hi=Real.~ Real.maxFinite} else let open Real fun roundup (lo, hi, 1) = if lo > halfMax then {lo=Real.maxFinite, hi=Real.posInf} else if hi > halfMax then {lo=2.0*lo, hi=Real.posInf} else {lo=2.0*lo, hi=2.0*hi} | roundup (lo, hi, i) = if lo > halfMax then {lo=Real.maxFinite, hi=Real.posInf} else roundup(2.0*lo, 2.0*hi, Int.-(i, 1)) val (r as {lo=lo, hi=hi}) = roundup (lo*p960, hi*p960, Int.-(exp, 960)) in if (s = 1) then r else {lo=Real.~ hi, hi=Real.~ lo} end end end fun approx rm ({num=a, denom=b}) = let val s = I.sign a in if s = 0 then 0.0 else let val rm = case (s, rm) of (~1, IEEEReal.TO_POSINF) => IEEEReal.TO_NEGINF | (~1, IEEEReal.TO_NEGINF) => IEEEReal.TO_POSINF | (_, _) => rm val a = I.abs a val (al, bl) = (I.log2 a, I.log2 b) val bound = Real.precision + 1 val d = al - bl val (a', b') = if d > bound then (a, I.<<(b, Word.fromInt (d - bound))) else (I.<<(a, Word.fromInt(bound - d)), b) val (q, r) = I.quotrem(a', b') val (q', exp) = if (r = N.zero) then (q, d - bound) else (I.+(N.one, I.<<(q, 0w30)), (d - bound) - 30) val app = Real.fromManExp{man = N.approx rm q', exp = exp} in if s = 1 then app else Real.~ app end end fun toFrac x = x fun op + (xx, yy) = let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy val d1 = N.gcd (a, b) in if d1 = N.one then {num = I.+(I.*(x, b), I.*(y, a)), denom = I.*(a, b)} else let val (a', b') = (I.quot(a, d1), I.quot(b, d1)) val t = I.+(I.*(x, b'), I.*(y, a')) val d2 = N.gcd (t, d1) in {num = I.quot(t, d2), denom = I.*(a', I.quot(b, d2))} end end fun op ~ ({num=x, denom=a}) = {num=I.~ x, denom=a} fun op - (xx, yy) = op + (xx, op ~ yy) fun op * (xx, yy) = let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy val (d1, d2) = (N.gcd(x, b), N.gcd(y, a)) in {num = I.*(I.quot(x,d1), I.quot(y,d2)), denom = I.*(I.quot(a,d2), I.quot(b,d1))} end fun sign ({num=p, denom=_}) = I.sign p fun max (xx, yy) = case Int.compare(sign xx, sign yy) of LESS => yy | GREATER => xx | EQUAL => let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy in if I.>(I.*(x, b), I.*(y, a)) then xx else yy end fun min (xx, yy) = case Int.compare(sign xx, sign yy) of LESS => xx | GREATER => yy | EQUAL => let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy in if I.>(I.*(x, b), I.*(y, a)) then yy else xx end fun compare (xx, yy) = case Int.compare(sign xx, sign yy) of EQUAL => let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy in I.compare(I.*(x, b), I.*(y, a)) end | s => s fun inv ({num=x, denom=a}) = if x = N.zero then raise Div else {num = a, denom = x} fun op / (xx, yy) = op * (xx, inv yy) fun == (xx, yy) = case Int.compare(sign xx, sign yy) of EQUAL => let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy in I.*(x, b) = I.*(y, a) end | _ => false fun != (xx, yy) = not (==(xx, yy)) fun op > (xx, yy) = case Int.compare(sign xx, sign yy) of EQUAL => let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy in I.>(I.*(x, b), I.*(y, a)) end | LESS => false | GREATER => true fun op < (xx, yy) = yy > xx fun op <= (xx, yy) = case Int.compare(sign xx, sign yy) of EQUAL => let val {num=x, denom=a} = toFrac xx val {num=y, denom=b} = toFrac yy in I.<=(I.*(x, b), I.*(y, a)) end | LESS => true | GREATER => false fun op >= (xx, yy) = yy <= xx fun abs ({num=x, denom=a}) = {num=I.abs x, denom=a} fun fromReal x = let val {odd=m, exp=e} = N.fromReal x in if Int.>=(e, 0) then {num = I.<<(m, Word.fromInt e), denom = N.one} else {num = m, denom = I.<<(N.one, Word.fromInt (Int.~ e))} end fun toString xx = let val {num=x, denom=a} = toFrac xx in (I.toString x) ^ " / " ^ (I.toString a) end end