# Library Coq.Numbers.Cyclic.DoubleCyclic.DoubleSqrt

Set Implicit Arguments.
Require Import ZArith.
Require Import BigNumPrelude.
Require Import DoubleType.
Require Import DoubleBase.
Local Open 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,
w_compare x y = Z.compare [|x|] [|y|].
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|] / (Z.pow 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,
ww_compare x y = Z.compare [[x]] [[y]].
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 spec_w_WW spec_w_sub
spec_w_add_mul_div spec_ww_Bm1 spec_w_add_c : 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.