(***********************************************
 ** InfiniteInteger.sml
 ** sml
 **
 ** Aleksandar Nanevski
 **
 ** IntInf structure extended with operations
 ** for shifting and conversion to Reals. Not
 ** to be used directly in the programs, but
 ** rather accessed only through the IntInf
 ** and NumberTheory structures.
 ***********************************************)


structure InfiniteInteger :> 
sig
    include INT_INF
    
    val zero : int
    val one : int
    val intervalApprox : int -> {lo : real, hi : real}
    val approx : IEEEReal.rounding_mode -> int -> real
    val toOddExp : int -> {odd : int, exp : Int.int}
    val fromReal : real -> {odd : int, exp : Int.int} 
end =
struct
    open IntInfBasis

    val itow = BN.itow
    val wtoi = BN.wtoi
    val lgBase = BN.lgBase
	
    (* the code assumes lgBase > 2 *)
	
    val wlgBase = itow lgBase
	
    val correction = itow (Int.-(Word.wordSize, lgBase))
	
	
    (* assumes 0 < i < BN.lgBase *) 
    fun digitRightShift (xs, i:Word.word) = 
	let val j = Word.-(itow Word.wordSize, i)
	    fun shift' ([], 0) = []
	      | shift' ([], d) = [d]
	      | shift' (x::xt, d) = 
		let val xw = itow x (* 30 bits always *)
		    val xr = wtoi (Word.>> (xw, i))
		    val xl = wtoi (Word.>> ((Word.<< (xw, j), correction)))
		in 
		    Int.+(xl, d)::shift' (xt, xr)
		end
	in
	    case xs of
		nil => nil
	      | x::xt => 
		    shift' (xt, wtoi (Word.>> (itow x, i)))
	end
    
    
    fun bnRightShift (xs, i:Word.word) =
	let val (q, r) = (Word.div (i, wlgBase), 
			  Word.mod (i, wlgBase))
	in
	    (* safe to do wtoi; lgBase>2 so q cannot be negative *)
	    digitRightShift (List.drop (xs, wtoi q), r) handle Subscript => nil
	end
    
    
    (* assumes 0 < i < BN.lgBase *) 
    (* assumes xs <> []          *)

    fun digitLeftShift (xs, i:Word.word) = 
	let val k = Word.+(correction, i)
	    val j = Word.-(itow Word.wordSize, k)
	    fun shift' ([], 0) = []
	      |	shift' ([], d) = [d]
	      | shift' (x::xt, d) = 
		let val xw = itow x (* 30 bits always *)
		    val xr = wtoi (Word.>> (xw, j))
		    val xl = wtoi (Word.>> ((Word.<< (xw, k), correction)))
		in 
		    Int.+(xl, d)::shift' (xt, xr)
		end
	in
	    shift' (xs, 0)
	end
    
    
    fun bnLeftShift (xs, i:Word.word) =
	let fun pad' (l, 0w0) = l	
	      | pad' (l, k) = pad' (0::l, Word.-(k,0w1))
	    val (q, r) = (Word.div (i, wlgBase), 
			  Word.mod (i, wlgBase))
	in
	    case xs of
		nil => nil
	      | _ => pad'(digitLeftShift (xs, r), q)
	end
    
    
    fun ~>> (BI{sign=s, digits=xs}, i) =
	let val ys = bnRightShift(xs, i)
	in
	    if (ys = BN.zero) then 
		BI{sign=POS, digits=ys}
	    else
		BI{sign=s, digits=ys}
	end
    
    fun << (BI{sign=s, digits=xs}, i) =
	let val ys = bnLeftShift(xs, i)
	in
	    BI{sign=s, digits=ys}
	end
    
    
    
    fun bnApprox rm xs =
	(* sum' sums all the elements even *)
	(* though not all are always needed *)
	let 
	    fun sum' (nil, k, acc) = acc
	      | sum' (x::xt, k, acc) =  
		let val acc = 
		    if x = 0 then acc
		    else Real.+(acc, 
				Real.fromManExp{man=Real.fromInt x, 
						exp=k})
		in
		    sum' (xt, Int.+(k, lgBase), acc)
		end
	    val rm' = IEEEReal.getRoundingMode()
	in
	    (IEEEReal.setRoundingMode(rm);
	     sum'(xs, 0, 0.0) before
	     IEEEReal.setRoundingMode(rm'))
	end
    


    val e = Real.fromManExp{man=1.0, exp= ~53}
    val p30 = Real.fromManExp{man=1.0, exp= Int.~ lgBase}
    val p60 = Real.*(p30, p30)
    val p990 = Real.fromManExp{man=1.0, exp=990}
    val halfMax = Real.*(0.5, Real.maxFinite)

    fun correction xx =
	let fun lsb' (0w0, i) = Real.*(i, e)
	      | lsb' (x, i) = lsb' (Word.>>(x, 0w1), Real.*(2.0,i))
	in
	    lsb' (Word.fromInt xx, 1.0)
	end


    (* doesn't handle xs = nil *)
    fun bnIntervalApprox xs =
	let fun zeroScan ([x], n) = 
		let val xx = Real.fromManExp{man=Real.fromInt x, exp=n}
		in
		    {lo=xx, hi=xx, exp=n}
		end
	      | zeroScan ([y, x], n) = 
		let open Real
		    val n = Int.+(n, lgBase)
		    val xx = fromInt x
		    val yy = p30 * Real.fromInt y
		    val rr = xx + yy
		in
		    case Real.compare (yy, rr - xx) of
			  LESS => {lo=rr - correction x, hi=rr, exp=n}
			| EQUAL => {lo=rr, hi=rr, exp=n}
			| GREATER => {lo=rr, hi=rr + correction x, exp=n}
		end
	      | zeroScan ([z, y, x], n) = 
		let open Real
		    val n = Int.+(n, Int.*(2, lgBase))
		    val xx = fromInt x 
		    val yy = p30 * Real.fromInt y
		    val rr = xx + yy
		in
		     case (Real.compare (yy, rr - xx)) of
			  EQUAL => 
			      let open Real
				  val zz = p60 * fromInt z
				  val ss = rr + zz
			      in
				  case (Real.compare(zz, ss - rr)) of
				      EQUAL =>  {lo=ss, hi=ss, exp=n}
				    | LESS => {lo=ss - correction x, hi=ss, exp=n}
				    | GREATER => {lo=ss, hi=ss + correction x, exp=n}
			      end
			| LESS => {lo=rr - correction x, hi=rr, exp=n}
			| GREATER => {lo=rr, hi=rr + correction x, exp=n}
		end
	      | zeroScan (x::t, n) = 
		if x = 0 then zeroScan (t, Int.+(n, lgBase))
		else positiveScan (t, Int.+(n, lgBase))

	    and positiveScan ([z, y, x], n) =  
	        let open Real
		    val n = Int.+(n, Int.*(2, lgBase))
		    val xx = fromInt x
		    val yy = p30 * Real.fromInt y
		    val rr = xx + yy
		in
		    case Real.compare(yy, rr - xx) of
			EQUAL => 
			    let val zz = p60 * fromInt z
				val ss = rr + zz
			    in
				case Real.compare(zz, ss - rr) of
				    EQUAL => {lo=ss, hi=ss+correction x, exp=n}
				  | LESS => {lo=ss-correction x, hi=ss, exp=n}
				  | GREATER => {lo=ss, hi=ss+correction x, exp=n}
			    end
		      | LESS => {lo=rr - correction x, hi=rr, exp=n}
		      | GREATER => {lo=rr, hi=rr + correction x, exp=n}
		end
	      | positiveScan (x::t, n) = positiveScan(t, Int.+(n, lgBase))
	    val {lo=lo, hi=hi, exp=n} = zeroScan (xs, 0)
	in 
	    if Int.<=(n, 990) then 
		let open Real
		    val p = fromManExp{man=1.0, exp=n}
		in
		    {lo=lo*p, hi=hi*p}
		end
	    else
		if Int.>(n, 1023) then {lo=Real.maxFinite, hi=Real.posInf}
		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))
		    in
			roundup (lo*p990, hi*p990, Int.-(n, 990))
		    end
	end
    

    fun intervalApprox (BI{sign=POS, digits=xs}) =
	(case xs of
	     nil => {lo=0.0, hi=0.0}
	   | _ => bnIntervalApprox xs)
      | intervalApprox (BI{sign=NEG, digits=xs}) =
	let val {lo=lo, hi=hi} = bnIntervalApprox xs
	in
	    {lo=Real.~ hi, hi=Real.~ lo}
	end


    fun approx rm (BI{sign=POS, digits=xs}) = 
	bnApprox rm xs
      | approx rm (BI{sign=NEG, digits=xs}) =
	(case rm of
	     IEEEReal.TO_POSINF =>
		 Real.~(bnApprox IEEEReal.TO_NEGINF xs)
	   | IEEEReal.TO_NEGINF => 
		 Real.~(bnApprox IEEEReal.TO_POSINF xs)
	   | _ => 
		 Real.~(bnApprox rm xs))


	     
    (* assumes xs <> [] *)
    fun bnDiv2 xs = 
	let fun div2' (0, acc) = acc
	      | div2' (x, acc) = 
	        if (Int.mod(x, 2) <> 0) then acc
		else div2' (Int.div(x, 2), Int.+(acc, 1))
	    fun pow2' ([], acc) = ([], acc)
	      | pow2' (0::xt, acc) = pow2'(xt, Int.+(lgBase, acc))
	      | pow2' (xs as (x::xt), acc) = 
		let val k = div2'(x, 0)
		in
		    (bnRightShift(xs, itow k), Int.+(acc, k))
		end
	in
	    pow2' (xs, 0)
	end
    
    
    (* works even if xx = 0 *)
    fun toOddExp (xx as (BI{sign=s, digits=xs})) = 
	let val (xs', k) = bnDiv2 xs
	in
	    {odd=BI{sign=s, digits=xs'}, exp=k}
	end


    (* maxInt + 1 *)
    val q30 = Real.fromManExp {man=1.0, exp=30}

    fun fromReal x =
	(* assumes positive x *)
	let fun manexp'(m, k, acc) = 
	    let val {whole=w, frac=f} = 
		Real.split (Real.*(m, q30)) (* result < maxInt + 1 *)
	    in 
		if Real.==(f, 0.0) then 
		    ((Real.floor w)::acc, Int.-(k, lgBase))
		else manexp'(f, Int.-(k, lgBase), (Real.floor w)::acc)
	    end
	    val s = Real.sign x
	in
	    if Real.==(x, 0.0) then 
		{odd = BI{sign=POS, digits=BN.zero}, exp = 0}
	    else 
		if Real.>(x, 0.0) then
		    if Real.==(x, Real.posInf) then raise Overflow
		    else
			let val {man=m, exp=exp} = Real.toManExp x
			    val (w, k) = manexp' (m, exp, nil)
			    val (w',k') = bnDiv2 w
			in
			    {odd = BI{sign=POS, digits=w'}, 
			     exp = Int.+(k, k')}
			end
		else
		    if Real.<(x, 0.0) then
			if Real.==(x, Real.negInf) then raise Overflow
			else
			    let val {man=m, exp=exp} = Real.toManExp (Real.~ x)
				val (w, k) = manexp' (m, exp, nil)
				val (w', k') = bnDiv2 w
			    in
				{odd = BI{sign=NEG, digits=w'}, 
				 exp = Int.+(k, k')}
			    end
		    else raise Domain (* NaN *)
	end
end