datatype Float = Float of (IntInf.int*IntInf.int*int);
fun fbase (Float (_,b,_)) = b;

fun pow b a = IntInf.pow (b,a);
(*
solves b^(x-1) < m + 1 <= b^x
*)
fun ilog b m = 
    let fun least f n = if (f n) then n else least f (n + 1) in
        least (fn k => (m div (pow b k) = 0)) 1
    end

fun fmulsign m1 m2 (Float (m,b,e)) = 
    if IntInf.sameSign (m1,m2) then (Float (m,b,e))
    else (Float (~m,b,e));

fun norm p (f as (Float (m,b,e))) = 
    let fun nn (f as (Float (m,b,e))) = 
            if m mod b = 0 then nn (Float (m div b,b,e+1)) else f
        val m' = IntInf.abs m in
        if m' < pow b p then nn f
        else let val x = ilog b m' in
                 fmulsign m 1 (nn (Float (m' div (pow b (x - p)), b, e + x - p)))
             end
    end

fun finv p (Float (m,b,e)) = 
    let val x = ilog b m in
        norm p (Float ((pow b (x+p)) div m, b, ~x-p-e))
    end

exception BaseNotEq;

fun fmul p (Float (m1',b1,e1)) (Float (m2',b2,e2)) = 
    let val m1 = IntInf.abs m1' val m2 = IntInf.abs m2' in
        if b1 <> b2 then raise BaseNotEq
        else fmulsign m1' m2' (norm p (Float (m1*m2, b1, e1+e2)))
    end

fun fdiv p (Float (m1',b1,e1)) (Float (m2',b2,e2)) = 
    let val m1 = IntInf.abs m1' val m2 = IntInf.abs m2' in
        if b1 <> b2 then raise BaseNotEq
        else
            let val x = ilog b1 m2 in
                fmulsign m1' m2' (norm p (Float (m1*(pow b1 (x+p)) div m2, b1, ~x-p-e2+e1)))
            end
    end

fun fadd p (Float (m1,b1,e1)) (Float (m2,b2,e2)) = 
    if b1 <> b2 then raise BaseNotEq
    else if e1 <= e2 then
        norm p (Float (m1*(pow b1 (e2-e1)) + m2*(pow b1 (2*(e2-e1))), b1, 2*e1-e2))
    else
        norm p (Float (m2*(pow b1 (e1-e2)) + m1*(pow b1 (2*(e1-e2))), b1, 2*e2-e1));

fun fsub p f1 (Float (m,b,e)) = fadd p f1 (Float (~m,b,e));

fun fabs (Float (m,b,e)) = Float (IntInf.abs m,b,e);
fun convBase p (Float (m,b,e)) b' = if e < 0 then
                                        fdiv p (Float (m,b',0)) (Float (pow b (~e),b',0))
                                    else norm p (Float (m*(pow b e),b',0));
fun fcompare (Float (m1,b1,e1)) (Float (m2,b2,e2)) = 
    if b1 <> b2 then raise BaseNotEq
    else if e1 <= e2 then IntInf.compare (m1*(pow b1 (e2-e1)), m2*(pow b1 (2*(e2-e1))))
    else IntInf.compare (m1*(pow b1 (2*(e1-e2))), m2*(pow b1 (e1-e2)));

fun fsqrtnext p (f as Float (m,b,e)) x = fdiv p (fadd p (fdiv p f x) x) (Float (2,b,0));

fun fsqrtnthapprox p (f as Float (m,b,e)) n = 
    let fun ff m x = if m >= n then x else ff (m+1) (fsqrtnext p f x) in
        ff 0 (Float (1,b,0))
    end

val test1 = map (fsqrtnthapprox 3 (Float (3,2,0))) [2,3,4];

fun fsqrtn a x e = 
    let fun fs x n = 
            let val x' = fsqrtnext (2*e) a x in
                case fcompare (fabs (fsub (2*e) x' x)) (Float (1,fbase a,~e)) of
                    LESS => (x',n)
                  | _ => fs x' (n+1)
            end in
        fs x 0
    end

fun fsqrt a x e = #1 (fsqrtn a x e);

val sqrt2 = fsqrt (Float (2,10,0)) (Float (1,10,0)) 100;
val test = fsub 200 (fmul 200 sqrt2 sqrt2) (Float (2,10,0));

fun ntimes' f x 0 = [x]
  | ntimes' f x n = let val nt = ntimes' f x (n-1) in
                        (f (List.hd nt)) :: nt
                    end

fun ntimes f x n = rev (ntimes' f x n);

let val x = norm 70 (fsqrt (Float (2,10,0)) (Float (1,10,0)) 70) in
    norm 1 (fsub 200 (fmul 200 x x) (Float (2,10,0)))
end