sqrt.mlw 2.89 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
use lemmas.Lemmas
use mach.fxp.Fxp

13 14
meta coercion function rval

15
val rsa_estimate (a: fxp): fxp
16
  requires { 0.25 <=. a <=. 0xffffffffffffffffp-64 }
17
  requires { iexp a = - 64 }
18 19
  ensures { iexp result = -8 }
  ensures { 256 <= ival result <= 511 }
20 21 22
  ensures { 1. <=. result <=. 2. }
  ensures { let rsa = 1. /. sqrt a in
            let e0 = (result -. rsa) /. rsa in -0.00281 <=. e0 <=. 0.002655 }
23

24 25
let sqrt1 (rp: ptr uint64) (a0: uint64): uint64
  requires { valid rp 1 }
26 27
  requires { 0x4000000000000000 <= a0 }
  ensures { result * result <= a0 < (result + 1) * (result + 1) }
28 29
  ensures { result * result + (pelts rp)[offset rp] = a0 }
  ensures { (pelts rp)[offset rp] <= 2 * result }
30 31
=
  let a = fxp_init a0 (-64) in
32 33 34
  assert { 0.25 <=. a <=. 0xffffffffffffffffp-64 };
  assert { 0. <. a };
  let rsa = pure { 1. /. sqrt a } in
35
  let x0 = rsa_estimate a in
36
  let e0 = pure { (x0 -. rsa) /. rsa } in
37
  let a1 = fxp_lsr a 31 in
38
  let ea1 = pure { (a1 -. a) /. a } in
39 40
  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
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
43 44 45
  let mx1 = pure { x0 +. x0 *. t1' *. 0.5 } in
  assert { (mx1 -. rsa) /. rsa = -0.5 *. (e0 *. e0 *. (3. +. e0) +. (1. +. e0) *. (1. -. m1 +. (1. +. e0) *. (1. +. e0) *. ea1)) };
  let e1 = pure { (x1 -. rsa) /. rsa } in
46
  let a2 = fxp_lsr a 24 in
47
  let ea2 = pure { (a2 -. a) /. a } in
48
  let u1 = fxp_mul x1 a2 in
49 50
  let eu1 = pure { (u1 -. sqrt a) /. sqrt a } in
  assert { eu1 = (1. +. e1) *. (1. +. ea2) -. 1. };
51
  let u2 = fxp_lsr u1 25 in
52
  let eu2 = pure { (u2 -. u1) /. u1 } in
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
55
  assert { sqrt a *. sqrt a = a };
56
  let t2 = fxp_asr t2' 24 in
57
  let x2 = fxp_add u1 (fxp_asr' (fxp_mul x1 t2) 15 1) in
58
  let mx2 = pure { u1 +. x1 *. t2' *. 0.5 } in
59 60
  assert { let ee = (1. +. eu1) *. (1. +. eu1) *. eu2 *. (2. +. eu2) +. eu1 *. eu1 in
           (mx2 -. sqrt a) /. sqrt a = -0.5 *. (ee +. m2 /. a) *. (1. +. e1) -. e1 *. eu1
61
           by x1 <> 0. /\ a2 <> 0. };
62
  let x = fxp_lsr x2 32 in
63 64
  let sa = pure { trunc_at (sqrt a) (-32) } in
  assert { -0x1.p-32 <=. x -. sa <=. 0. };
65 66 67 68
  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
69
           by x <=. sa <=. 1.
70 71 72 73 74 75 76 77 78 79 80 81
           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