sqrt.mlw 2.72 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
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 };
22
  assert { 0. <. rval a };
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
  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 u1 = fxp_mul x1 a2 in
  assert { rsa *. rval a = sa };
  let u2 = fxp_lsr u1 25 in
44
  assert { rval a -. ((rval x1 *. rval a) *. (rval x1 *. rval a)) = -. rval a *. e1 *. (2. +. e1) };
45
  let t2 = fxp_asr (fxp_sub (fxp_sub (fxp_lsl a 14) (fxp_mul u2 u2)) (fxp_init 0x1000000000 (-78))) 24 in
46
  let ghost m2 = 0x1.p-38 in
47
  let ghost mt2 = rval a -. (rval x1 *. rval a2) *. (rval x1 *. rval a2) -. m2 in
48
  let x2 = fxp_add u1 (fxp_asr' (fxp_mul x1 t2) 15 1) in
49
  let ghost mx2 = rval x1 *. rval a2 +. rval x1 *. mt2 *. 0.5 in
50
  assert { (mx2 -. sa) /. sa = 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) };
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
  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