(********************************************
 ** UnNormalizedRational.sml
 ** sml
 **
 ** Aleksandar Nanevski 
 **
 ** Implementation of infinite precision
 ** rational numbers, with changeable 
 ** normalization strategy
 ********************************************)

structure UnNormalizedRational :> 
sig 
    include RATIONAL

    val normalize : rat -> rat

    datatype normalization_mode = ALWAYS | NEVER 
    val setNormalizationMode : normalization_mode -> unit
    val getNormalizationMode : unit -> normalization_mode
end =

struct
    structure I = IntInf
    structure N = NumberTheory
	
    datatype normalization_mode = ALWAYS | NEVER 

    val mode = ref NEVER

    (* assumes denom > 0 *)
    type rat = {num: I.int, denom: I.int, norm: bool} ref
    type t = rat

    fun setNormalizationMode m = (mode := m)
    fun getNormalizationMode () = !mode


    fun fromIntInf x = 
	ref {num = x, denom = N.one, norm = true}

    fun fromInt x = 
	ref {num = I.fromInt x, denom = N.one, norm = true}

    val zero = fromInt 0
    val one = fromInt 1
	
    fun normalize (xx as ref {num=p, denom=q, norm=norm}) =
	if norm then xx
	else 
	    let val d = N.gcd(p, q)
	    in
		xx before
		xx := 
		(if d = N.one then {num=p, denom=q, norm=true}
		 else {num = I.quot(p, d),
		       denom = I.quot(q, d),
		       norm = true})
	    end
	

    fun fromFrac {num=x, denom=y} = 
	if y = N.one then 
	    ref {num=x, denom=N.one,  norm=true}
	else 
	    let val xx = ref {num=x, denom=y, norm=false}
	    in
		case !mode of
		    ALWAYS => normalize xx
		  | NEVER => xx
	    end


    fun fromIntFrac {num=x, denom=y} = 
	if y = 1 then
	    ref {num=I.fromInt x, denom=N.one, norm=true}
	else
	    let val xx = ref {num=I.fromInt x, 
			      denom=I.fromInt y, 
			      norm=false}
	    in
		case !mode of 
		    ALWAYS => normalize xx
		  | NEVER => xx
	    end


    fun fromManExp {man=man, exp=exp} =
	if Int.>(exp, 0) then
	    ref {num = I.<<(man, Word.fromInt exp),
		 denom = N.one,
		 norm = true}
	else
	    ref {num = man,
		 denom = I.<<(N.one, Word.fromInt (Int.~ exp)),
		 norm = true}



    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)

    fun intervalApprox (ref {num=a, denom=b, norm=_}) = 
	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 (ref {num=a, denom=b, norm=_}) =
	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 (ref {num=x, denom=a, norm=_}) = {num=x, denom=a}
	
	
    fun op + (xx, yy) = 
	case !mode of
	    ALWAYS =>
		let val {num=x, denom=a} = toFrac(normalize xx)
		    val {num=y, denom=b} = toFrac(normalize yy)
		    val d1 = N.gcd (a, b)
		in
		    if d1 = N.one then 
			ref {num = I.+(I.*(x, b), I.*(y, a)),
			     denom = I.*(a, b),
			     norm = true}
		    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
			    ref {num = I.quot(t, d2),
				 denom = I.*(a', I.quot(b, d2)),
				 norm= true}
			end
		end
	  | NEVER => 
		let val {num=x, denom=a} = toFrac xx
		    val {num=y, denom=b} = toFrac yy
		in
		    ref {num = I.+(I.*(x, b), I.*(y, a)),
			 denom = I.*(a, b),
			 norm = false}
		end


    fun op ~ (xx as ref {num=x, denom=a, norm=norm}) = 
	case !mode of
	    ALWAYS => 
		let val {num=x, denom=a} = toFrac(normalize xx)
		in
		    ref {num=I.~ x, denom=a, norm=true}
		end
	  | NEVER =>
		ref {num=I.~ x, denom=a, norm=norm}
	    
	    
	    
    fun op - (xx, yy) = op + (xx, op ~ yy)
	

    fun op * (xx, yy) = 
	case !mode of
	    ALWAYS => 
		let val {num=x, denom=a} = toFrac(normalize xx)
		    val {num=y, denom=b} = toFrac(normalize yy)
		    val (d1, d2) = (N.gcd(x, b), N.gcd(y, a))
		in
		    ref {num = I.*(I.quot(x,d1), I.quot(y,d2)),
			 denom = I.*(I.quot(a,d2), I.quot(b,d1)),
			 norm = true}
		end
	  | NEVER =>
		let val {num=x, denom=a} = toFrac xx
		    val {num=y, denom=b} = toFrac yy
		in
		    ref {num = I.*(x, y),
			 denom = I.*(a, b),
			 norm = false}
		end


    fun sign (ref {num=p, denom=_, norm=_}) = I.sign p
	

    fun max (xx, yy) =
	case Int.compare(sign xx, sign yy) of
	    LESS => (case !mode of ALWAYS => normalize yy | NEVER => yy)
	  | GREATER => (case !mode of ALWAYS => normalize xx | NEVER => xx)
	  | EQUAL => 
		(case !mode of
		     ALWAYS =>
			 let val {num=x, denom=a} = toFrac(normalize xx)
			     val {num=y, denom=b} = toFrac(normalize yy)
			 in
			     if I.>(I.*(x, b), I.*(y, a)) then xx
			     else yy
			 end
		   | NEVER => 
			 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 => (case !mode of ALWAYS => normalize xx | NEVER => xx)
	  | GREATER => (case !mode of ALWAYS => normalize yy | NEVER => yy)
	  | EQUAL => 	    
		(case !mode of
		     ALWAYS =>
			 let val {num=x, denom=a} = toFrac(normalize xx)
			     val {num=y, denom=b} = toFrac(normalize yy)
			 in
			     if I.>(I.*(x, b), I.*(y, a)) then yy
			     else xx
			 end
		   | NEVER => 
			 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 (xx as ref {num=x, denom=a, norm=norm}) =
	if x = N.zero then raise Div 
	else
	    case !mode of
		ALWAYS =>
		    let val {num=x, denom=a} = toFrac(normalize xx)
		    in
			ref {num = a, denom = x, norm = true}
		    end
	      | NEVER => 
		    ref {num = a, denom = x, norm = norm}
		    

    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 (xx as ref {num=x, denom=a, norm=norm}) =
	case !mode of
	    ALWAYS => 
		let val {num=x, denom=a} = toFrac(normalize xx)
		in
		    ref {num=I.abs x, denom=a, norm=true}
		end
	  | NEVER => 
		ref {num=I.abs x, denom=a, norm=norm}
		



    fun fromReal x =
	let val {odd=m, exp=e} =  N.fromReal x 
	in 
	    if Int.>=(e, 0) then 
		ref {num = I.<<(m, Word.fromInt e),
		     denom = N.one,
		     norm = true}
	    else
		ref {num = m,
		     denom = I.<<(N.one, Word.fromInt (Int.~ e)),
		     norm = true}
	end



    fun toString xx = 
	let val {num=x, denom=a} = toFrac(normalize xx)
	in
	    (I.toString x) ^ " / " ^ (I.toString a)
	end
end