# 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.