Require Import ZArith.
Require Import BigNumPrelude.
Require Import DoubleType.
Require Import DoubleBase.
Open Local Scope Z_scope.
Section DoubleSqrt.
Variable w :
Type.
Variable w_is_even :
w ->
bool.
Variable w_compare :
w ->
w ->
comparison.
Variable w_0 :
w.
Variable w_1 :
w.
Variable w_Bm1 :
w.
Variable w_WW :
w ->
w ->
zn2z w.
Variable w_W0 :
w ->
zn2z w.
Variable w_0W :
w ->
zn2z w.
Variable w_sub :
w ->
w ->
w.
Variable w_sub_c :
w ->
w ->
carry w.
Variable w_square_c :
w ->
zn2z w.
Variable w_div21 :
w ->
w ->
w ->
w *
w.
Variable w_add_mul_div :
w ->
w ->
w ->
w.
Variable w_digits :
positive.
Variable w_zdigits :
w.
Variable ww_zdigits :
zn2z w.
Variable w_add_c :
w ->
w ->
carry w.
Variable w_sqrt2 :
w ->
w ->
w *
carry w.
Variable w_pred :
w ->
w.
Variable ww_pred_c :
zn2z w ->
carry (
zn2z w).
Variable ww_pred :
zn2z w ->
zn2z w.
Variable ww_add_c :
zn2z w ->
zn2z w ->
carry (
zn2z w).
Variable ww_add :
zn2z w ->
zn2z w ->
zn2z w.
Variable ww_sub_c :
zn2z w ->
zn2z w ->
carry (
zn2z w).
Variable ww_add_mul_div :
zn2z w ->
zn2z w ->
zn2z w ->
zn2z w.
Variable ww_head0 :
zn2z w ->
zn2z w.
Variable ww_compare :
zn2z w ->
zn2z w ->
comparison.
Variable low :
zn2z w ->
w.
Let wwBm1 :=
ww_Bm1 w_Bm1.
Definition ww_is_even x :=
match x with
|
W0 =>
true
|
WW xh xl =>
w_is_even xl
end.
Let w_div21c x y z :=
match w_compare x z with
|
Eq =>
match w_compare y z with
Eq => (
C1 w_1,
w_0)
|
Gt => (
C1 w_1,
w_sub y z)
|
Lt => (
C1 w_0,
y)
end
|
Gt =>
let x1 :=
w_sub x z in
let (
q,
r) :=
w_div21 x1 y z in
(
C1 q,
r)
|
Lt =>
let (
q,
r) :=
w_div21 x y z in
(
C0 q,
r)
end.
Let w_div2s x y s :=
match x with
C1 x1 =>
let x2 :=
w_sub x1 s in
let (
q,
r) :=
w_div21c x2 y s in
match q with
C0 q1 =>
if w_is_even q1 then
(
C0 (
w_add_mul_div (
w_pred w_zdigits)
w_1 q1),
C0 r)
else
(
C0 (
w_add_mul_div (
w_pred w_zdigits)
w_1 q1),
w_add_c r s)
|
C1 q1 =>
if w_is_even q1 then
(
C1 (
w_add_mul_div (
w_pred w_zdigits)
w_0 q1),
C0 r)
else
(
C1 (
w_add_mul_div (
w_pred w_zdigits)
w_0 q1),
w_add_c r s)
end
|
C0 x1 =>
let (
q,
r) :=
w_div21c x1 y s in
match q with
C0 q1 =>
if w_is_even q1 then
(
C0 (
w_add_mul_div (
w_pred w_zdigits)
w_0 q1),
C0 r)
else
(
C0 (
w_add_mul_div (
w_pred w_zdigits)
w_0 q1),
w_add_c r s)
|
C1 q1 =>
if w_is_even q1 then
(
C0 (
w_add_mul_div (
w_pred w_zdigits)
w_1 q1),
C0 r)
else
(
C0 (
w_add_mul_div (
w_pred w_zdigits)
w_1 q1),
w_add_c r s)
end
end.
Definition split x :=
match x with
|
W0 => (
w_0,w_0)
|
WW h l => (
h,l)
end.
Definition ww_sqrt2 x y :=
let (
x1,
x2) :=
split x in
let (
y1,
y2) :=
split y in
let (
q,
r) :=
w_sqrt2 x1 x2 in
let (
q1,
r1) :=
w_div2s r y1 q in
match q1 with
C0 q1 =>
let q2 :=
w_square_c q1 in
let a :=
WW q q1 in
match r1 with
C1 r2 =>
match ww_sub_c (
WW r2 y2)
q2 with
C0 r3 => (
a,
C1 r3)
|
C1 r3 => (
a,
C0 r3)
end
|
C0 r2 =>
match ww_sub_c (
WW r2 y2)
q2 with
C0 r3 => (
a,
C0 r3)
|
C1 r3 =>
let a2 :=
ww_add_mul_div (
w_0W w_1)
a W0 in
match ww_pred_c a2 with
C0 a3 =>
(
ww_pred a,
ww_add_c a3 r3)
|
C1 a3 =>
(
ww_pred a,
C0 (
ww_add a3 r3))
end
end
end
|
C1 q1 =>
let a1 :=
WW q w_Bm1 in
let a2 :=
ww_add_mul_div (
w_0W w_1)
a1 wwBm1 in
(
a1,
ww_add_c a2 y)
end.
Definition ww_is_zero x :=
match ww_compare W0 x with
Eq =>
true
|
_ =>
false
end.
Definition ww_head1 x :=
let p :=
ww_head0 x in
if (
ww_is_even p)
then p else ww_pred p.
Definition ww_sqrt x :=
if (
ww_is_zero x)
then W0
else
let p :=
ww_head1 x in
match ww_compare p W0 with
|
Gt =>
match ww_add_mul_div p x W0 with
W0 =>
W0
|
WW x1 x2 =>
let (
r,
_) :=
w_sqrt2 x1 x2 in
WW w_0 (
w_add_mul_div
(
w_sub w_zdigits
(
low (
ww_add_mul_div (
ww_pred ww_zdigits)
W0 p)))
w_0 r)
end
|
_ =>
match x with
W0 =>
W0
|
WW x1 x2 =>
WW w_0 (
fst (
w_sqrt2 x1 x2))
end
end.
Variable w_to_Z :
w ->
Z.
Notation wB := (
base w_digits).
Notation wwB := (
base (
ww_digits w_digits)).
Notation "[| x |]" := (
w_to_Z x) (
at level 0,
x at level 99).
Notation "[+| c |]" :=
(
interp_carry 1
wB w_to_Z c) (
at level 0,
x at level 99).
Notation "[-| c |]" :=
(
interp_carry (-1)
wB w_to_Z c) (
at level 0,
x at level 99).
Notation "[[ x ]]" := (
ww_to_Z w_digits w_to_Z x)(
at level 0,
x at level 99).
Notation "[+[ c ]]" :=
(
interp_carry 1
wwB (
ww_to_Z w_digits w_to_Z)
c)
(
at level 0,
x at level 99).
Notation "[-[ c ]]" :=
(
interp_carry (-1)
wwB (
ww_to_Z w_digits w_to_Z)
c)
(
at level 0,
x at level 99).
Notation "[|| x ||]" :=
(
zn2z_to_Z wwB (
ww_to_Z w_digits w_to_Z)
x) (
at level 0,
x at level 99).
Notation "[! n | x !]" := (
double_to_Z w_digits w_to_Z n x)
(
at level 0,
x at level 99).
Variable spec_w_0 : [|w_0|] = 0.
Variable spec_w_1 : [|w_1|] = 1.
Variable spec_w_Bm1 : [|w_Bm1|] =
wB - 1.
Variable spec_w_zdigits : [|w_zdigits|] =
Zpos w_digits.
Variable spec_more_than_1_digit: 1 <
Zpos w_digits.
Variable spec_ww_zdigits : [[
ww_zdigits]] =
Zpos (
xO w_digits).
Variable spec_to_Z :
forall x, 0 <= [|x|] <
wB.
Variable spec_to_w_Z :
forall x, 0 <= [[
x]] <
wwB.
Variable spec_w_WW :
forall h l, [[
w_WW h l]] = [|h|] *
wB + [|l|].
Variable spec_w_W0 :
forall h, [[
w_W0 h]] = [|h|] *
wB.
Variable spec_w_0W :
forall l, [[
w_0W l]] = [|l|].
Variable spec_w_is_even :
forall x,
if w_is_even x then [|x|]
mod 2 = 0
else [|x|]
mod 2 = 1.
Variable spec_w_compare :
forall x y,
match w_compare x y with
|
Eq => [|x|] = [|y|]
|
Lt => [|x|] < [|y|]
|
Gt => [|x|] > [|y|]
end.
Variable spec_w_sub :
forall x y, [|w_sub
x y|] = ([|x|] - [|y|])
mod wB.
Variable spec_w_square_c :
forall x, [[
w_square_c x]] = [|x|] * [|x|].
Variable spec_w_div21 :
forall a1 a2 b,
wB/2 <= [|b|] ->
[|a1|] < [|b|] ->
let (
q,r) :=
w_div21 a1 a2 b in
[|a1|] *wB+ [|a2|] = [|q|] * [|b|] + [|r|] /\
0 <= [|r|] < [|b|].
Variable spec_w_add_mul_div :
forall x y p,
[|p|] <=
Zpos w_digits ->
[|
w_add_mul_div p x y |] =
([|x|] * (2 ^ [|p|]) +
[|y|] / (
Zpower 2 ((
Zpos w_digits) - [|p|])))
mod wB.
Variable spec_ww_add_mul_div :
forall x y p,
[[
p]] <=
Zpos (
xO w_digits) ->
[[
ww_add_mul_div p x y ]] =
([[
x]] * (2^ [[
p]]) +
[[
y]] / (2^ (
Zpos (
xO w_digits) - [[
p]])))
mod wwB.
Variable spec_w_add_c :
forall x y, [+|w_add_c
x y|] = [|x|] + [|y|].
Variable spec_ww_add :
forall x y, [[
ww_add x y]] = ([[
x]] + [[
y]])
mod wwB.
Variable spec_w_sqrt2 :
forall x y,
wB/ 4 <= [|x|] ->
let (
s,r) :=
w_sqrt2 x y in
[[
WW x y]] = [|s|] ^ 2 + [+|r|] /\
[+|r|] <= 2 * [|s|].
Variable spec_ww_sub_c :
forall x y, [-[
ww_sub_c x y]] = [[
x]] - [[
y]].
Variable spec_ww_pred_c :
forall x, [-[
ww_pred_c x]] = [[
x]] - 1.
Variable spec_pred :
forall x, [|w_pred
x|] = ([|x|] - 1)
mod wB.
Variable spec_ww_pred :
forall x, [[
ww_pred x]] = ([[
x]] - 1)
mod wwB.
Variable spec_ww_add_c :
forall x y, [+[
ww_add_c x y]] = [[
x]] + [[
y]].
Variable spec_ww_compare :
forall x y,
match ww_compare x y with
|
Eq => [[
x]] = [[
y]]
|
Lt => [[
x]] < [[
y]]
|
Gt => [[
x]] > [[
y]]
end.
Variable spec_ww_head0 :
forall x, 0 < [[
x]] ->
wwB/ 2 <= 2 ^ [[
ww_head0 x]] * [[
x]] <
wwB.
Variable spec_low:
forall x, [|low
x|] = [[
x]]
mod wB.
Let spec_ww_Bm1 : [[
wwBm1]] =
wwB - 1.
Hint Rewrite spec_w_0 spec_w_1 w_Bm1 spec_w_WW spec_w_sub
spec_w_div21 spec_w_add_mul_div spec_ww_Bm1
spec_w_add_c spec_w_sqrt2:
w_rewrite.
Lemma spec_ww_is_even :
forall x,
if ww_is_even x then [[
x]]
mod 2 = 0
else [[
x]]
mod 2 = 1.
Theorem spec_w_div21c :
forall a1 a2 b,
wB/2 <= [|b|] ->
let (
q,r) :=
w_div21c a1 a2 b in
[|a1|] *
wB + [|a2|] = [+|q|] * [|b|] + [|r|] /\ 0 <= [|r|] < [|b|].
Theorem C0_id:
forall p, [+|C0
p|] = [|p|].
Theorem add_mult_div_2:
forall w,
[|w_add_mul_div (
w_pred w_zdigits)
w_0 w|] = [|w|] / 2.
Theorem add_mult_div_2_plus_1:
forall w,
[|w_add_mul_div (
w_pred w_zdigits)
w_1 w|] =
[|w|] / 2 + 2 ^
Zpos (
w_digits - 1).
Theorem add_mult_mult_2:
forall w,
[|w_add_mul_div
w_1 w w_0|] = 2 * [|w|]
mod wB.
Theorem ww_add_mult_mult_2:
forall w,
[[
ww_add_mul_div (
w_0W w_1)
w W0]] = 2 * [[
w]]
mod wwB.
Theorem ww_add_mult_mult_2_plus_1:
forall w,
[[
ww_add_mul_div (
w_0W w_1)
w wwBm1]] =
(2 * [[
w]] + 1)
mod wwB.
Theorem Zplus_mod_one:
forall a1 b1, 0 <
b1 -> (
a1 +
b1)
mod b1 =
a1 mod b1.
Lemma C1_plus_wB:
forall x, [+|C1
x|] =
wB + [|x|].
Theorem spec_w_div2s :
forall a1 a2 b,
wB/2 <= [|b|] -> [+|a1|] <= 2 * [|b|] ->
let (
q,r) :=
w_div2s a1 a2 b in
[+|a1|] *
wB + [|a2|] = [+|q|] * (2 * [|b|]) + [+|r|] /\ 0 <= [+|r|] < 2 * [|b|].
Theorem wB_div_4: 4 * (
wB / 4) =
wB.
Theorem Zsquare_mult:
forall p,
p ^ 2 =
p *
p.
Theorem Zsquare_pos:
forall p, 0 <=
p ^ 2.
Lemma spec_split:
forall x,
[|fst (
split x)|] *
wB + [|snd (
split x)|] = [[
x]].
Theorem mult_wwB:
forall x y, [|x|] * [|y|] <
wwB.
Hint Resolve mult_wwB.
Lemma spec_ww_sqrt2 :
forall x y,
wwB/ 4 <= [[
x]] ->
let (
s,r) :=
ww_sqrt2 x y in
[||WW
x y||] = [[
s]] ^ 2 + [+[
r]] /\
[+[
r]] <= 2 * [[
s]].
Lemma spec_ww_is_zero:
forall x,
if ww_is_zero x then [[
x]] = 0
else 0 < [[
x]].
Lemma wwB_4_2: 2 * (
wwB / 4) =
wwB/ 2.
Lemma spec_ww_head1
:
forall x :
zn2z w,
(
ww_is_even (
ww_head1 x) =
true) /\
(0 < [[
x]] ->
wwB / 4 <= 2 ^ [[
ww_head1 x]] * [[
x]] <
wwB).
Theorem wwB_4_wB_4:
wwB / 4 =
wB / 4 *
wB.
Lemma spec_ww_sqrt :
forall x,
[[
ww_sqrt x]] ^ 2 <= [[
x]] < ([[
ww_sqrt x]] + 1) ^ 2.
End DoubleSqrt.