Commit bb95fc9f authored by Guillaume Melquiond's avatar Guillaume Melquiond Committed by Raphael Rieu-Helft

Add early version of square root algorithm.

parent e2c99a74
use int.Int
use real.RealInfix
use real.Square
use mach.int.UInt64
use int.Power
use lemmas.Lemmas
use mach.fxp.Fxp
val rsa_estimate (a: fxp): fxp
ensures { iexp result = -8 }
ensures { 256 <= ival result <= 511 }
ensures { 1. <=. rval result <=. 2. }
ensures { let rsa = 1. /. sqrt(rval a) in
let e0 = (rval result -. rsa) /. rsa in -0.00281 <=. e0 <=. 0.002655 }
let sqrt1 (a0: uint64): uint64
requires { 0x4000000000000000 <= a0 }
ensures { result * result <= a0 < (result + 1) * (result + 1) }
=
let a = fxp_init a0 (-64) in
assert { 0.25 <=. rval a <=. 0xffffffffffffffffp-64 };
let ghost sa = sqrt (rval a) in
let ghost rsa = 1. /. sa in
let x0 = rsa_estimate a in
let ghost e0 = (rval x0 -. rsa) /. rsa in
let a1 = fxp_lsr a 31 in
let ghost ea1 = (rval a1 -. rval a) /. rval a in
let t1' =
fxp_sub
(fxp_sub (fxp_init 0x2000000000000 (-49)) (fxp_init 0x30000 (-49)))
(fxp_mul (fxp_mul x0 x0) a1) in
let ghost m1 = 0x3.p-33 in
let t1 = fxp_asr t1' 16 in
let x1 = fxp_add (fxp_lsl x0 16) (fxp_asr' (fxp_mul x0 t1) 18 1) in
let ghost mx1 = rval x0 +. rval x0 *. rval t1' *. 0.5 in
assert { (mx1 -. rsa) /. rsa = -0.5 *. (e0 *. e0 *. (3. +. e0) +. (1. +. e0) *. (m1 +. (1. +. e0) *. (1. +. e0) *. ea1)) };
let ghost e1 = (rval x1 -. rsa) /. rsa in
let a2 = fxp_lsr a 24 in
let ghost ea2 = (rval a2 -. rval a) /. rval a in
let a3 = fxp_lsl a 14 in
let ghost ea3 = (rval a3 -. rval a) /. rval a in
let u1 = fxp_mul x1 a2 in
assert { rsa *. rval a = sa };
let u2 = fxp_lsr u1 25 in
let t2' = fxp_sub (fxp_sub a3 (fxp_mul u2 u2)) (fxp_init 0x10000000000 (-78)) in
let ghost m2 = 0x1.p-38 in
let t2 = fxp_asr t2' 24 in
let x2 = fxp_add u1 (fxp_asr' (fxp_mul x1 t2) 15 1) in
let ghost mx2 = rval x1 *. rval a2 +. rval x1 *. rval t2' *. 0.5 in
assert { (mx2 -. sa) /. sa = 0.5 *. ea3 *. (1. +. e1) +. ea2 *. (1. +. e1) *. (-. e1 *. (2. +. e1) -. (1. +. e1) *. (1. +. e1) *. 0.5 *. ea2) -. 0.5 *. m2 /. rval a *. (1. +. e1) -. e1 *. e1 *. 0.5 *. (3. +. e1) };
let x = fxp_lsr x2 32 in
let ghost sa = trunc_at sa (-32) in
assert { -0x1.p-32 <=. rval x -. sa <=. 0. };
let ref c = ival x in
assert { c * c <= a0 < (c+2) * (c+2) };
let ref s = c * c in
assert { (c+1) * (c+1) <= radix
by rval x <=. sa <=. 1.
so iexp x = -32
so c < power 2 32
so c+1 <= power 2 32
so (c+1) * (c+1) <= power 2 32 * power 2 32 = radix };
assert { s + c <= s + c + c < (c+1) * (c+1) <= radix };
if (s + c + c <= a0 - 1)
then begin
assert { (c+1) * (c+1) <= a0 };
s <- s + 2 * c + 1;
c <- c + 1;
assert { s = c * c };
end;
c
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE why3session PUBLIC "-//Why3//proof session v5//EN"
"http://why3.lri.fr/why3session.dtd">
<why3session shape_version="5">
<prover id="0" name="Alt-Ergo" version="2.2.0" timelimit="5" steplimit="0" memlimit="1000"/>
<prover id="2" name="Gappa" version="1.3.2" timelimit="5" steplimit="0" memlimit="1000"/>
<file>
<path name=".."/>
<path name="sqrt.mlw"/>
<theory name="Top">
<goal name="VC sqrt1" expl="VC for sqrt1">
<transf name="split_vc" >
<goal name="VC sqrt1.0" expl="assertion">
</goal>
<goal name="VC sqrt1.1" expl="precondition" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.2" expl="precondition" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.3" expl="fxp overflow" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.4" expl="precondition" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.5" expl="fxp alignment" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.6" expl="fxp alignment" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.7" expl="fxp overflow" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.8" expl="fxp overflow" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.9" expl="fxp alignment" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.10" expl="assertion">
</goal>
<goal name="VC sqrt1.11" expl="precondition" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.12" expl="fxp overflow" proved="true">
<proof prover="2"><result status="valid" time="0.01"/></proof>
</goal>
<goal name="VC sqrt1.13" expl="precondition" proved="true">
<proof prover="2"><result status="valid" time="0.01"/></proof>
</goal>
<goal name="VC sqrt1.14" expl="precondition" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.15" expl="assertion">
</goal>
<goal name="VC sqrt1.16" expl="fxp overflow" proved="true">
<proof prover="2"><result status="valid" time="0.31"/></proof>
</goal>
<goal name="VC sqrt1.17" expl="fxp alignment" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.18" expl="fxp alignment" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.19" expl="fxp overflow">
<proof prover="2"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.20" expl="fxp overflow">
<proof prover="2"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.21" expl="fxp alignment" proved="true">
<proof prover="2"><result status="valid" time="0.00"/></proof>
</goal>
<goal name="VC sqrt1.22" expl="assertion">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.23" expl="fxp overflow" proved="true">
<proof prover="2"><result status="valid" time="1.10"/></proof>
</goal>
<goal name="VC sqrt1.24" expl="assertion">
<proof prover="2"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.25" expl="assertion">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.26" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.27" expl="assertion">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.28" expl="assertion" proved="true">
<proof prover="0"><result status="valid" time="0.52" steps="177"/></proof>
</goal>
<goal name="VC sqrt1.29" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.30" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.31" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.32" expl="assertion" proved="true">
<proof prover="0"><result status="valid" time="0.57" steps="187"/></proof>
</goal>
<goal name="VC sqrt1.33" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.34" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.35" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.36" expl="integer overflow">
<proof prover="0"><result status="timeout" time="5.00"/></proof>
</goal>
<goal name="VC sqrt1.37" expl="assertion" proved="true">
<proof prover="0"><result status="valid" time="0.60" steps="196"/></proof>
</goal>
<goal name="VC sqrt1.38" expl="postcondition" proved="true">
<proof prover="0"><result status="valid" time="0.73" steps="197"/></proof>
</goal>
<goal name="VC sqrt1.39" expl="postcondition" proved="true">
<proof prover="0"><result status="valid" time="0.41" steps="187"/></proof>
</goal>
</transf>
</goal>
</theory>
</file>
</why3session>
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment