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, c at level 99).
  Notation "[-| c |]" :=
   (interp_carry (-1) wB w_to_Z c) (at level 0, c 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, c at level 99).
  Notation "[-[ c ]]" :=
   (interp_carry (-1) wwB (ww_to_Z w_digits w_to_Z) c)
   (at level 0, c 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.