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