sqrt.mlw 3.07 KB
Newer Older
1 2
module Sqrt1

3 4 5 6 7
use int.Int
use real.RealInfix
use real.Square
use mach.int.UInt64
use int.Power
8 9
use map.Map
use mach.c.C
10 11 12 13
use lemmas.Lemmas
use mach.fxp.Fxp

val rsa_estimate (a: fxp): fxp
14 15
  requires { 0.25 <=. rval a <=. 0xffffffffffffffffp-64 }
  requires { iexp a = - 64 }
16 17 18 19 20 21
  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 }

22 23
let sqrt1 (rp: ptr uint64) (a0: uint64): uint64
  requires { valid rp 1 }
24 25
  requires { 0x4000000000000000 <= a0 }
  ensures { result * result <= a0 < (result + 1) * (result + 1) }
26 27
  ensures { result * result + (pelts rp)[offset rp] = a0 }
  ensures { (pelts rp)[offset rp] <= 2 * result }
28 29 30
=
  let a = fxp_init a0 (-64) in
  assert { 0.25 <=. rval a <=. 0xffffffffffffffffp-64 };
31
  assert { 0. <. rval a };
32 33 34 35 36 37
  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
38 39
  let m1 = fxp_sub (fxp_init 0x2000000000000 (-49)) (fxp_init 0x30000 (-49)) in
  let t1' = fxp_sub m1 (fxp_mul (fxp_mul x0 x0) a1) in
40 41 42
  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
43
  assert { (mx1 -. rsa) /. rsa = -0.5 *. (e0 *. e0 *. (3. +. e0) +. (1. +. e0) *. (1. -. rval m1 +. (1. +. e0) *. (1. +. e0) *. ea1)) };
44 45 46 47 48 49
  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
50
  let ghost eu2 = (rval u2 -. rval u1) /. rval u1 in
51
  assert { rval a -. ((rval x1 *. rval a) *. (rval x1 *. rval a)) = -. rval a *. e1 *. (2. +. e1) };
52 53 54
  let m2 = fxp_init 0x24000000000 (-78) in
  let t2' = fxp_sub (fxp_sub (fxp_lsl a 14) (fxp_mul u2 u2)) m2 in
  let t2 = fxp_asr t2' 24 in
55
  let x2 = fxp_add u1 (fxp_asr' (fxp_mul x1 t2) 15 1) in
56 57 58 59 60 61
  let ghost mx2 = rval u1 +. rval x1 *. rval t2' *. 0.5 in
  assert { let p1a2 = (1. +. e1) *. (1. +. ea2) in
           let e1a2 = p1a2 -. 1. in
           let ee = p1a2 *. p1a2 *. eu2 *. (eu2 +. 2.) +. e1a2 *. e1a2 in
           (mx2 -. sa) /. sa = -0.5 *. (ee +. rval m2 /. rval a) *. (1. +. e1) -. e1 *. e1a2
           by rval x1 <> 0. /\ rval a2 <> 0. };
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
  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;
82
  C.set rp (a0 - s);
83
  c
84 85

end