2012-07-28 01:32:38 +02:00
|
|
|
|
|
|
|
|
(* Computing the square root (up to eps) using binary search
|
|
|
|
|
|
|
|
|
|
Note: precondition 0. < eps is not necessary for soundness
|
2014-01-23 22:21:40 +01:00
|
|
|
but needed to ensure termination.
|
|
|
|
|
*)
|
2012-07-28 01:32:38 +02:00
|
|
|
|
|
|
|
|
module BinarySqrt
|
|
|
|
|
|
2018-06-15 16:45:31 +02:00
|
|
|
use real.Real
|
2018-06-15 15:10:50 +02:00
|
|
|
use int.Int as Int
|
|
|
|
|
use real.MinMax as MinMax
|
|
|
|
|
use real.FromInt as FromInt
|
|
|
|
|
use real.Truncate as Truncate
|
2012-07-28 01:32:38 +02:00
|
|
|
|
2014-01-23 22:21:40 +01:00
|
|
|
let rec sqrt (r: real) (eps: real) (ghost n:int) (ghost eps0:real) : real
|
|
|
|
|
requires { 0.0 <= r }
|
|
|
|
|
requires { eps0 > 0.0 /\ Int.(>=) n 1 }
|
|
|
|
|
requires { eps = (FromInt.from_int n) * eps0 }
|
|
|
|
|
variant { Int.(-) (Truncate.ceil (MinMax.max r 1.0 / eps0)) n }
|
2012-10-13 00:23:29 +02:00
|
|
|
ensures { result * result <= r < (result + eps) * (result + eps) }
|
2014-01-23 22:21:40 +01:00
|
|
|
= if r < eps && 1.0 < eps then
|
|
|
|
|
0.0
|
2012-07-28 01:32:38 +02:00
|
|
|
else
|
2014-01-23 22:21:40 +01:00
|
|
|
begin
|
|
|
|
|
assert { FromInt.from_int n * eps0 <= MinMax.max r 1.0 };
|
|
|
|
|
assert { 1.0 / eps0 > 0.0 };
|
|
|
|
|
assert { FromInt.from_int n * eps0 * (1.0 / eps0) <= MinMax.max r 1.0 * (1.0 / eps0)};
|
|
|
|
|
assert { FromInt.from_int n * eps0 / eps0 <= MinMax.max r 1.0 / eps0};
|
|
|
|
|
assert { FromInt.from_int n <= MinMax.max r 1.0 / eps0 };
|
|
|
|
|
let r' = sqrt r (2.0 * eps) (Int.(*) 2 n) eps0 in
|
2012-07-28 01:32:38 +02:00
|
|
|
if (r' + eps) * (r' + eps) <= r then r' + eps else r'
|
2014-01-23 22:21:40 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
|
|
let sqrt_main (r:real) (eps:real) : real
|
|
|
|
|
requires { 0.0 <= r }
|
|
|
|
|
requires { eps > 0.0 }
|
|
|
|
|
ensures { result * result <= r < (result + eps) * (result + eps) }
|
|
|
|
|
= sqrt r eps 1 eps
|
2012-07-28 01:32:38 +02:00
|
|
|
|
|
|
|
|
end
|