Commit 0dd92f60 authored by Raphaël Rieu-Helft's avatar Raphaël Rieu-Helft

Merge branch 'fxp_sqrt2' into 'master'

GMP square root

See merge request !106
parents 8ecd13e9 bac70072
......@@ -16,6 +16,10 @@ module ref.Ref
syntax val (:=) "%1 = %2" prec 14 13 14
end
module mach.int.Bounded_int
syntax val of_int "%1" prec 0
end
module mach.int.Int32
syntax type int32 "int32_t"
......@@ -66,7 +70,7 @@ module mach.int.UInt32GMP
prelude
"
#define LOW_MASK 0x00000000FFFFFFFFULL
#define LOW_MASK 0x00000000FFFFFFFFUL
struct __add32_with_carry_result
{ uint32_t __field_0;
......@@ -205,8 +209,8 @@ struct __lsld32_result lsld32(uint32_t x, uint32_t cnt);
"(uint32_t)((((uint64_t)%1) | (((uint64_t)%2) << 32))/(uint64_t)(%3))"
prec 2
syntax val lsl "%1 << %2" prec 5 5 4
syntax val lsr "%1 >> %2" prec 5 5 4
syntax val lsl "%1 << %2" prec 5 5 2
syntax val lsr "%1 >> %2" prec 5 5 2
syntax val is_msb_set "%1 & 0x80000000U" prec 8 8
......@@ -219,7 +223,7 @@ end
module mach.int.Int64
syntax type int64 "int64_t"
syntax literal int64 "%dLL"
syntax literal int64 "%dL"
syntax val (+) "%1 + %2" prec 4 4 3
syntax val (-) "%1 - %2" prec 4 4 3
......@@ -246,7 +250,7 @@ end
module mach.int.UInt64
syntax literal uint64 "0x%16xULL"
syntax literal uint64 "0x%16xUL"
syntax val (+) "%1 + %2" prec 4 4 3
syntax val (-) "%1 - %2" prec 4 4 3
......@@ -501,9 +505,9 @@ static struct __lsld64_result lsld64(uint64_t x, uint64_t cnt)
return result;
}
"
syntax literal uint64 "0x%16xULL"
syntax literal uint64 "0x%16xUL"
syntax val uint64_max "0xffffffffffffffffULL" prec 0
syntax val uint64_max "0xffffffffffffffffUL" prec 0
syntax val (+) "%1 + %2" prec 4 4 3
syntax val (-) "%1 - %2" prec 4 4 3
......@@ -534,17 +538,22 @@ static struct __lsld64_result lsld64(uint64_t x, uint64_t cnt)
syntax val sub_mod "%1 - %2" prec 4 4 3
syntax val mul_mod "%1 * %2" prec 3 3 2
syntax val lsl "%1 << %2" prec 5 5 4
syntax val lsr "%1 >> %2" prec 5 5 4
syntax val lsl "%1 << %2" prec 5 5 2
syntax val lsr "%1 >> %2" prec 5 5 2
syntax val lsl_mod "%1 << %2" prec 5 5 2
syntax val lsr_mod "%1 >> %2" prec 5 5 2
syntax val is_msb_set "%1 & 0x8000000000000000ULL" prec 8 7
syntax val is_msb_set "%1 & 0x8000000000000000UL" prec 8 7
syntax val count_leading_zeros "__builtin_clzll(%1)" prec 1 15
syntax val to_int32 "(int32_t)%1" prec 2 2
syntax val of_int32 "(uint64_t)%1" prec 2 2
syntax val to_int64 "(int64_t)%1" prec 2 2
syntax val of_int64 "(uint64_t)%1" prec 2 2
syntax val of_int "%1" prec 0
end
module mach.array.Array32
......@@ -593,3 +602,15 @@ module mach.c.C
syntax val print_uint32 "printf(\"%#010x\",%1)" prec 1 15
end
module mach.fxp.Fxp
syntax val fxp_add "%1 + %2" prec 4 4 3
syntax val fxp_sub "%1 - %2" prec 4 4 3
syntax val fxp_mul "%1 * %2" prec 3 3 2
syntax val fxp_lsl "%1 << %2" prec 5 5 2
syntax val fxp_lsr "%1 >> %2" prec 5 5 2
syntax val fxp_asr "(uint64_t)((int64_t)%1 >> %2)" prec 2 1 2
syntax val fxp_asr' "(uint64_t)((int64_t)%1 >> %2)" prec 2 1 2
end
\ No newline at end of file
......@@ -23,6 +23,7 @@ transformation "eliminate_let"
transformation "simplify_formula"
transformation "simplify_unknown_lsymbols"
transformation "simplify_trivial_quantification"
transformation "simplify_computations"
transformation "introduce_premises"
transformation "instantiate_predicate"
transformation "abstract_unknown_lsymbols"
......@@ -43,7 +44,7 @@ theory int.Int
syntax function (+) "(%1 + %2)"
syntax function (-) "(%1 - %2)"
syntax function ( * ) "(%1 * %2)"
syntax function (*) "(%1 * %2)"
syntax function (-_) "(-%1)"
syntax predicate (<=) "dummy"
......@@ -56,6 +57,7 @@ theory int.Int
meta "gappa arith" predicate (<), "not ", ">=", "<="
meta "gappa arith" predicate (>), "not ", "<=", ">="
meta "inline:no" function (-)
meta "inline:no" predicate (<=)
meta "inline:no" predicate (>=)
meta "inline:no" predicate (>)
......@@ -87,7 +89,7 @@ theory real.Real
syntax function (+) "(%1 + %2)"
syntax function (-) "(%1 - %2)"
syntax function ( * ) "(%1 * %2)"
syntax function (*) "(%1 * %2)"
syntax function (/) "(%1 / %2)"
syntax function (-_) "(-%1)"
syntax function inv "(1.0 / %1)"
......@@ -102,6 +104,7 @@ theory real.Real
meta "gappa arith" predicate (<), "not ", ">=", "<="
meta "gappa arith" predicate (>), "not ", "<=", ">="
meta "inline:no" function (-)
meta "inline:no" predicate (<=)
meta "inline:no" predicate (>=)
meta "inline:no" predicate (>)
......@@ -215,3 +218,9 @@ theory floating_point.DoubleMultiRounding
meta "instantiate:auto" prop Bounded_value
end
theory mach.fxp.Fxp
syntax function trunc_at "fixed<%2,dn>(%1)"
end
......@@ -29,35 +29,35 @@ why3:
dir:
mkdir -p build
MLWFILES= $(addsuffix .mlw, toom logical div mul sub add compare util)
MLWFILES= $(addsuffix .mlw, sqrtrem sqrt toom logical div mul sub add compare util)
cfiles: why3 dir
$(WHY3) extract -D c -L . --recursive --modular --interface -o build/ \
$(WHY3) extract -D wmpn.drv -D c -L . --recursive --modular --interface -o build/ \
wmpn.mlw
#$(MLWFILES)
extract: why3 dir cfiles
CFILES = build/uint64gmp.c build/toom.c build/div.c build/logical.c build/mul.c build/sub.c build/add.c build/compare.c build/util.c build/int32.c
CFILES = build/uint64gmp.c build/fxp.c build/sqrt.c build/sqrt1.c build/toom.c build/div.c build/logical.c build/mul.c build/sub.c build/add.c build/compare.c build/util.c build/int32.c
tests: extract check-gmp
gcc $(CFLAGS) tests.c $(CFILES) -I$(GMP_DIR) -Irandom -L$(GMP_LIB) -lgmp -o build/tests
gcc $(CFLAGS) -DCOMPARE_MINI tests.c $(CFILES) -I$(GMP_DIR) -Irandom -Imini-gmp -o build/minitests
gcc $(CFLAGS) tests.c $(CFILES) -Iinclude -I$(GMP_DIR) -Irandom -L$(GMP_LIB) -lgmp -o build/tests
gcc $(CFLAGS) -DCOMPARE_MINI tests.c $(CFILES) -Iinclude -I$(GMP_DIR) -Irandom -Imini-gmp -o build/minitests
./build/tests
./build/minitests
bench-tests: extract
gcc $(CFLAGS) tests.c $(CFILES) -Ibench-include -Irandom -lgmp -o build/bench-tests
gcc $(CFLAGS) tests.c $(CFILES) -Iinclude -Ibench-include -Irandom -lgmp -o build/bench-tests
build/why3%bench: extract check-gmp
gcc $(CFLAGS) -DTEST_WHY3 -DTEST_`echo $* | tr [:lower:] [:upper:]` tests.c $(CFILES) -I$(GMP_DIR) -Irandom -L$(GMP_LIB) -lgmp -o $@
gcc $(CFLAGS) -DTEST_WHY3 -DTEST_`echo $* | tr [:lower:] [:upper:]` tests.c $(CFILES) -Iinclude -I$(GMP_DIR) -Irandom -L$(GMP_LIB) -lgmp -o $@
build/gmp%bench: extract check-gmp
gcc $(CFLAGS) -DTEST_GMP -DTEST_`echo $* | tr [:lower:] [:upper:]` tests.c $(CFILES) -I$(GMP_DIR) -Irandom -L$(GMP_LIB) -lgmp -o $@
gcc $(CFLAGS) -DTEST_GMP -DTEST_`echo $* | tr [:lower:] [:upper:]` tests.c $(CFILES) -Iinclude -I$(GMP_DIR) -Irandom -L$(GMP_LIB) -lgmp -o $@
build/minigmp%bench: extract
gcc $(CFLAGS) -DTEST_MINIGMP -DTEST_`echo $* | tr [:lower:] [:upper:]` tests.c $(CFILES) -I$(GMP_DIR) -Imini-gmp -Irandom -o $@
gcc $(CFLAGS) -DTEST_MINIGMP -DTEST_`echo $* | tr [:lower:] [:upper:]` tests.c $(CFILES) -Iinclude -I$(GMP_DIR) -Imini-gmp -Irandom -o $@
alltests: tests build/why3addbench build/why3mulbench build/why3toombench build/why3divbench build/gmpaddbench build/gmpmulbench build/gmptoombench build/gmpdivbench build/minigmpaddbench build/minigmpmulbench build/minigmpdivbench build/minigmptoombench
......
This diff is collapsed.
......@@ -76,7 +76,7 @@
</transf>
</goal>
<goal name="VC wmpn_cmp.16" expl="postcondition" proved="true">
<proof prover="3"><result status="valid" time="0.03"/></proof>
<proof prover="0"><result status="valid" time="0.16"/></proof>
</goal>
<goal name="VC wmpn_cmp.17" expl="assertion" proved="true">
<proof prover="3"><result status="valid" time="0.03"/></proof>
......@@ -103,7 +103,7 @@
</transf>
</goal>
<goal name="VC wmpn_cmp.21" expl="postcondition" proved="true">
<proof prover="0"><result status="valid" time="0.16"/></proof>
<proof prover="3"><result status="valid" time="0.03"/></proof>
</goal>
<goal name="VC wmpn_cmp.22" expl="loop variant decrease" proved="true">
<proof prover="3"><result status="valid" time="0.03"/></proof>
......
This diff is collapsed.
This diff is collapsed.
/* TODO GMP header stuff */
static const unsigned char invsqrttab[384] = /* The common 0x100 was removed */
{
0xff,0xfd,0xfb,0xf9,0xf7,0xf5,0xf3,0xf2, /* sqrt(1/80)..sqrt(1/87) */
0xf0,0xee,0xec,0xea,0xe9,0xe7,0xe5,0xe4, /* sqrt(1/88)..sqrt(1/8f) */
0xe2,0xe0,0xdf,0xdd,0xdb,0xda,0xd8,0xd7, /* sqrt(1/90)..sqrt(1/97) */
0xd5,0xd4,0xd2,0xd1,0xcf,0xce,0xcc,0xcb, /* sqrt(1/98)..sqrt(1/9f) */
0xc9,0xc8,0xc6,0xc5,0xc4,0xc2,0xc1,0xc0, /* sqrt(1/a0)..sqrt(1/a7) */
0xbe,0xbd,0xbc,0xba,0xb9,0xb8,0xb7,0xb5, /* sqrt(1/a8)..sqrt(1/af) */
0xb4,0xb3,0xb2,0xb0,0xaf,0xae,0xad,0xac, /* sqrt(1/b0)..sqrt(1/b7) */
0xaa,0xa9,0xa8,0xa7,0xa6,0xa5,0xa4,0xa3, /* sqrt(1/b8)..sqrt(1/bf) */
0xa2,0xa0,0x9f,0x9e,0x9d,0x9c,0x9b,0x9a, /* sqrt(1/c0)..sqrt(1/c7) */
0x99,0x98,0x97,0x96,0x95,0x94,0x93,0x92, /* sqrt(1/c8)..sqrt(1/cf) */
0x91,0x90,0x8f,0x8e,0x8d,0x8c,0x8c,0x8b, /* sqrt(1/d0)..sqrt(1/d7) */
0x8a,0x89,0x88,0x87,0x86,0x85,0x84,0x83, /* sqrt(1/d8)..sqrt(1/df) */
0x83,0x82,0x81,0x80,0x7f,0x7e,0x7e,0x7d, /* sqrt(1/e0)..sqrt(1/e7) */
0x7c,0x7b,0x7a,0x79,0x79,0x78,0x77,0x76, /* sqrt(1/e8)..sqrt(1/ef) */
0x76,0x75,0x74,0x73,0x72,0x72,0x71,0x70, /* sqrt(1/f0)..sqrt(1/f7) */
0x6f,0x6f,0x6e,0x6d,0x6d,0x6c,0x6b,0x6a, /* sqrt(1/f8)..sqrt(1/ff) */
0x6a,0x69,0x68,0x68,0x67,0x66,0x66,0x65, /* sqrt(1/100)..sqrt(1/107) */
0x64,0x64,0x63,0x62,0x62,0x61,0x60,0x60, /* sqrt(1/108)..sqrt(1/10f) */
0x5f,0x5e,0x5e,0x5d,0x5c,0x5c,0x5b,0x5a, /* sqrt(1/110)..sqrt(1/117) */
0x5a,0x59,0x59,0x58,0x57,0x57,0x56,0x56, /* sqrt(1/118)..sqrt(1/11f) */
0x55,0x54,0x54,0x53,0x53,0x52,0x52,0x51, /* sqrt(1/120)..sqrt(1/127) */
0x50,0x50,0x4f,0x4f,0x4e,0x4e,0x4d,0x4d, /* sqrt(1/128)..sqrt(1/12f) */
0x4c,0x4b,0x4b,0x4a,0x4a,0x49,0x49,0x48, /* sqrt(1/130)..sqrt(1/137) */
0x48,0x47,0x47,0x46,0x46,0x45,0x45,0x44, /* sqrt(1/138)..sqrt(1/13f) */
0x44,0x43,0x43,0x42,0x42,0x41,0x41,0x40, /* sqrt(1/140)..sqrt(1/147) */
0x40,0x3f,0x3f,0x3e,0x3e,0x3d,0x3d,0x3c, /* sqrt(1/148)..sqrt(1/14f) */
0x3c,0x3b,0x3b,0x3a,0x3a,0x39,0x39,0x39, /* sqrt(1/150)..sqrt(1/157) */
0x38,0x38,0x37,0x37,0x36,0x36,0x35,0x35, /* sqrt(1/158)..sqrt(1/15f) */
0x35,0x34,0x34,0x33,0x33,0x32,0x32,0x32, /* sqrt(1/160)..sqrt(1/167) */
0x31,0x31,0x30,0x30,0x2f,0x2f,0x2f,0x2e, /* sqrt(1/168)..sqrt(1/16f) */
0x2e,0x2d,0x2d,0x2d,0x2c,0x2c,0x2b,0x2b, /* sqrt(1/170)..sqrt(1/177) */
0x2b,0x2a,0x2a,0x29,0x29,0x29,0x28,0x28, /* sqrt(1/178)..sqrt(1/17f) */
0x27,0x27,0x27,0x26,0x26,0x26,0x25,0x25, /* sqrt(1/180)..sqrt(1/187) */
0x24,0x24,0x24,0x23,0x23,0x23,0x22,0x22, /* sqrt(1/188)..sqrt(1/18f) */
0x21,0x21,0x21,0x20,0x20,0x20,0x1f,0x1f, /* sqrt(1/190)..sqrt(1/197) */
0x1f,0x1e,0x1e,0x1e,0x1d,0x1d,0x1d,0x1c, /* sqrt(1/198)..sqrt(1/19f) */
0x1c,0x1b,0x1b,0x1b,0x1a,0x1a,0x1a,0x19, /* sqrt(1/1a0)..sqrt(1/1a7) */
0x19,0x19,0x18,0x18,0x18,0x18,0x17,0x17, /* sqrt(1/1a8)..sqrt(1/1af) */
0x17,0x16,0x16,0x16,0x15,0x15,0x15,0x14, /* sqrt(1/1b0)..sqrt(1/1b7) */
0x14,0x14,0x13,0x13,0x13,0x12,0x12,0x12, /* sqrt(1/1b8)..sqrt(1/1bf) */
0x12,0x11,0x11,0x11,0x10,0x10,0x10,0x0f, /* sqrt(1/1c0)..sqrt(1/1c7) */
0x0f,0x0f,0x0f,0x0e,0x0e,0x0e,0x0d,0x0d, /* sqrt(1/1c8)..sqrt(1/1cf) */
0x0d,0x0c,0x0c,0x0c,0x0c,0x0b,0x0b,0x0b, /* sqrt(1/1d0)..sqrt(1/1d7) */
0x0a,0x0a,0x0a,0x0a,0x09,0x09,0x09,0x09, /* sqrt(1/1d8)..sqrt(1/1df) */
0x08,0x08,0x08,0x07,0x07,0x07,0x07,0x06, /* sqrt(1/1e0)..sqrt(1/1e7) */
0x06,0x06,0x06,0x05,0x05,0x05,0x04,0x04, /* sqrt(1/1e8)..sqrt(1/1ef) */
0x04,0x04,0x03,0x03,0x03,0x03,0x02,0x02, /* sqrt(1/1f0)..sqrt(1/1f7) */
0x02,0x02,0x01,0x01,0x01,0x01,0x00,0x00 /* sqrt(1/1f8)..sqrt(1/1ff) */
};
......@@ -70,6 +70,22 @@ module Lemmas
meta remove_prop axiom prod_compat_strict_lr
let lemma prod_compat_lr (a b c d:int)
requires { 0 <= a <= b }
requires { 0 <= c <= d }
ensures { a * c <= b * d }
= ()
meta remove_prop axiom prod_compat_lr
let lemma simp_compat_strict_l (a b c:int)
requires { 0 <= a * b < a * c }
requires { 0 < a }
ensures { b < c }
= ()
meta remove_prop axiom simp_compat_strict_l
(** {3 Integer value of a natural number} *)
(** `value_sub x n m` denotes the integer represented by
......
......@@ -6,6 +6,7 @@
<prover id="1" name="CVC3" version="2.4.1" timelimit="5" steplimit="0" memlimit="1000"/>
<prover id="2" name="CVC4" version="1.5" timelimit="1" steplimit="0" memlimit="1000"/>
<prover id="3" name="Z3" version="4.5.0" timelimit="5" steplimit="0" memlimit="1000"/>
<prover id="4" name="CVC4" version="1.6" timelimit="1" steplimit="0" memlimit="1000"/>
<prover id="5" name="Alt-Ergo" version="2.0.0" timelimit="5" steplimit="0" memlimit="1000"/>
<file proved="true">
<path name=".."/>
......@@ -61,6 +62,12 @@
<goal name="VC prod_compat_strict_lr" expl="VC for prod_compat_strict_lr" proved="true">
<proof prover="2"><result status="valid" time="0.04"/></proof>
</goal>
<goal name="VC prod_compat_lr" expl="VC for prod_compat_lr" proved="true">
<proof prover="4"><result status="valid" time="0.30"/></proof>
</goal>
<goal name="VC simp_compat_strict_l" expl="VC for simp_compat_strict_l" proved="true">
<proof prover="4"><result status="valid" time="0.24"/></proof>
</goal>
<goal name="VC value_sub" expl="VC for value_sub" proved="true">
<transf name="split_goal_right" proved="true" >
<goal name="VC value_sub.0" expl="variant decrease" proved="true">
......
......@@ -28,6 +28,35 @@ module Logical
ensures { mod (x * y + z) x = mod z x }
=
()
let lsl_mod_ext (x cnt: limb) : limb
requires { 0 <= cnt < Limb.length }
ensures { result = mod (x * power 2 cnt) radix }
ensures { result <= radix - power 2 cnt }
=
let r = lsl_mod x cnt in
let ghost p = power 2 (Limb.to_int cnt) in
let ghost q = power 2 (Limb.length - Limb.to_int cnt) in
assert { p * q = radix };
let ghost d = div (Limb.to_int x * p) radix in
assert { d * q >= 0 by d >= 0 so q >= 0 };
assert { mod r p = 0
by x * p = d * radix + r
so mod (x * p) p = mod (p * x + 0) p = mod 0 p = 0
so mod (d * radix + r) p = 0
so d * radix + r = p * (d * q) + r
so mod (d * radix + r) p = mod (p * (d * q) + r) p = mod r p };
assert { r <= radix - p
by mod r p = 0
so r < radix
so radix = p * power 2 (Limb.length - cnt)
so mod radix p = mod (p * q + 0) p = mod 0 p = 0
so let d1 = div r p in
let d2 = div radix p in
(r <= radix - p by
r = p * d1 so radix = p * d2 so p * d1 < p * d2 so p > 0
so d1 < d2 so d1 <= d2 - 1
so p * d1 <= p * (d2 - 1) = radix - p) };
r
let lsld_ext (x cnt:limb) : (limb,limb)
requires { 0 <= cnt < Limb.length }
......
This diff is collapsed.
module Sqrt1
use int.Int
use int.EuclideanDivision
use real.RealInfix
use real.Square
use real.FromInt
use real.Truncate as Trunc
use mach.int.UInt64
use int.Power
use map.Map
use mach.c.C
use lemmas.Lemmas
use mach.fxp.Fxp
meta coercion function rval
let lemma real_sqr_compat (x y: real)
requires { 0. <. x <. y }
ensures { x *. x <. y *. y }
= ()
meta remove_prop axiom real_sqr_compat
let lemma trunc_lower_bound (x: real) (k:int)
ensures { x <. trunc_at x k +. pow2 k }
= assert { x *. pow2 (-k)
<. from_int (Trunc.floor (x *. pow2 (-k)) + 1)
= from_int (Trunc.floor (x *. pow2 (-k))) +. 1.
= from_int (Trunc.floor (x *. pow2 (-k))) *. (pow2 k *. pow2 (-k))
+. 1.
= trunc_at x k *. pow2 (-k) +. 1.
= (trunc_at x k +. pow2 k) *. pow2 (-k) }
val rsa_estimate (a: fxp): fxp
requires { 0.25 <=. a <=. 0xffffffffffffffffp-64 }
requires { iexp a = - 64 }
ensures { iexp result = -8 }
ensures { 256 <= ival result <= 511 }
ensures { 1. <=. result <=. 2. }
ensures { let rsa = 1. /. sqrt a in
let e0 = (result -. rsa) /. rsa in -0.00281 <=. e0 <=. 0.002655 }
let sqrt1 (rp: ptr uint64) (a0: uint64): uint64
requires { valid rp 1 }
requires { 0x4000000000000000 <= a0 }
ensures { result * result <= a0 < (result + 1) * (result + 1) }
ensures { result * result + (pelts rp)[offset rp] = a0 }
ensures { (pelts rp)[offset rp] <= 2 * result }
=
let a = fxp_init a0 (-64) in
assert { 0.25 <=. a <=. 0xffffffffffffffffp-64 };
assert { 0. <. a };
let rsa = pure { 1. /. sqrt a } in
let x0 = rsa_estimate a in
let e0 = pure { (x0 -. rsa) /. rsa } in
let a1 = fxp_lsr a 31 in
let ea1 = pure { (a1 -. a) /. a } in
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
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 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
let a2 = fxp_lsr a 24 in
let ea2 = pure { (a2 -. a) /. a } in
let u1 = fxp_mul x1 a2 in
let eu1 = pure { (u1 -. sqrt a) /. sqrt a } in
assert { eu1 = (1. +. e1) *. (1. +. ea2) -. 1. };
let u2 = fxp_lsr u1 25 in
let eu2 = pure { (u2 -. u1) /. u1 } in
let m2 = fxp_init 0x24000000000 (-78) in
let t2' = fxp_sub (fxp_sub (fxp_lsl a 14) (fxp_mul u2 u2)) m2 in
assert { sqrt a *. sqrt a = a };
let t2 = fxp_asr t2' 24 in
let x2 = fxp_add u1 (fxp_asr' (fxp_mul x1 t2) 15 1) in
let mx2 = pure { u1 +. x1 *. t2' *. 0.5 } in
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
by x1 <> 0. /\ a2 <> 0. };
let x = fxp_lsr x2 32 in
let sa = pure { trunc_at (sqrt a) (-32) } in
assert { -0x1.p-32 <=. x -. sa <=. 0. };
let ref c = ival x in
assert { c = Trunc.floor (pow2 32 *. x)
by x <=. sa <=. 1.
so iexp x = -32
so pow2 32 *. x <. pow2 32
so 0 <= Trunc.floor (pow2 32 *. x)
< power 2 32
< radix
so c = mod (Trunc.floor (pow2 32 *. x)) radix
= Trunc.floor (pow2 32 *. x) };
assert { c * c <= a0 < (c+2) * (c+2)
by x *. x <=. a
so from_int a0 = pow2 64 *. a
so from_int c = pow2 32 *. x
so pow2 32 *. pow2 32 = pow2 64
so from_int (c * c) = from_int c *. from_int c
= pow2 64 *. (x *. x) <=. from_int a0
so c * c <= a0
so let x' = x +. 0x2.p-32 in
let sa' = sa +. pow2 (-32) in
pow2 (-32) = 0x1.p-32
so sa' <=. x'
so sqrt a <. trunc_at (sqrt a) (-32) +. pow2 (-32) = sa'
so a = sqrt a *. sqrt a <. sa' *. sa'<=. x' *. x'
so from_int (c + 2) = from_int c +. from_int 2
= from_int c +. 2.
= pow2 32 *. x +. 2.
= pow2 32 *. (x +. 0x2.p-32)
= pow2 32 *. x'
so from_int a0
= pow2 64 *. a
<. pow2 64 *. (x' *. x')
= from_int ((c+2) * (c+2))
};
let ref s = c * c in
assert { (c+1) * (c+1) <= radix };
assert { s + c <= s + c + c < (c+1) * (c+1) <= radix };
if (s + 2 * 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.set rp (a0 - s);
c
meta remove_prop axiom trunc_lower_bound
end
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -2092,5 +2092,22 @@ with wmpn_toom32_mul (r x y scratch: ptr limb) (sx sy:int32) (ghost k: int) : un
sfree scratch
end*)
meta remove_prop axiom no_borrow
meta remove_prop axiom no_borrow_ptr
meta remove_prop lemma power_ge_1
meta remove_prop axiom valuation_def
meta remove_prop axiom valuation_spec
meta remove_prop lemma valuation_mul
meta remove_prop axiom valuation_times_pow
meta remove_prop axiom valuation_lower_bound
meta remove_prop axiom valuation_monotonous
meta remove_prop lemma valuation_nondiv
meta remove_prop lemma valuation_zero_prod
meta remove_prop axiom valuation_times_nondiv
meta remove_prop lemma valuation_split
meta remove_prop lemma valuation_prod
meta remove_prop lemma valuation_mod
meta remove_prop lemma valuation_decomp
meta remove_prop lemma valuation_pow
end
\ No newline at end of file
......@@ -136,7 +136,7 @@
<proof prover="1"><result status="valid" time="0.02"/></proof>
</goal>
<goal name="VC wmpn_zero_p.8" expl="postcondition" proved="true">
<proof prover="2"><result status="valid" time="0.03" steps="45"/></proof>
<proof prover="2"><result status="valid" time="0.03" steps="43"/></proof>
</goal>
<goal name="VC wmpn_zero_p.9" expl="assertion" proved="true">
<proof prover="1"><result status="valid" time="0.01"/></proof>
......@@ -154,7 +154,7 @@
<proof prover="1"><result status="valid" time="0.01"/></proof>
</goal>
<goal name="VC wmpn_zero_p.14" expl="postcondition" proved="true">
<proof prover="0"><result status="valid" time="0.02" steps="28"/></proof>
<proof prover="0"><result status="valid" time="0.02" steps="27"/></proof>
</goal>
</transf>
</goal>
......
module sqrt.Sqrt1
prelude "
#include \"sqrtinit.h\"
uint64_t rsa_estimate (uint64_t a) {
uint64_t abits, x0;
abits = a >> 55;
x0 = 0x100 | invsqrttab[abits - 0x80];
return x0;
}
"
end
\ No newline at end of file
......@@ -10,5 +10,7 @@ module Wmpn
use export mul.Mul
use export div.Div
use export toom.Toom
use export sqrt.Sqrt1
use export sqrtrem.Sqrt
end
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -22,3 +22,5 @@ val normalize_hyp : int option -> Decl.prsymbol option -> Env.env
val normalize_hyp_few : int option -> Decl.prsymbol option -> Env.env
-> Task.task Trans.tlist
val simplify : (Term.lsymbol -> bool) -> Env.env -> Task.task Trans.trans
This diff is collapsed.
This diff is collapsed.
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