Zsqrt_compat.v

Built with Alectryon, running Coq+SerAPI v8.10.0+0.7.0. Coq sources are in this panel; goals and messages will appear in the other. Bubbles () indicate interactive fragments: hover for details, tap to reveal contents. Use Ctrl+↑ Ctrl+↓ to navigate, Ctrl+🖱️ to focus.

(************************************************************************)
(*         *   The Coq Proof Assistant / The Coq Development Team       *)
(*  v      *   INRIA, CNRS and contributors - Copyright 1999-2018       *)
(* <O___,, *       (see CREDITS file for the list of authors)           *)
(*   \VV/  **************************************************************)
(*    //   *    This file is distributed under the terms of the         *)
(*         *     GNU Lesser General Public License Version 2.1          *)
(*         *     (see LICENSE file for the text of the license)         *)
(************************************************************************)

Require Import ZArithRing.
Require Import Omega.
Require Export ZArith_base.
Local Open Scope Z_scope.

THIS FILE IS DEPRECATED

Instead of the various Zsqrt defined here, please use rather Z.sqrt (or Z.sqrtrem). The latter are pure functions without proof parts, and more results are available about them. Some equivalence proofs between the old and the new versions can be found below. Importing ZArith will provides by default the new versions.

(**********************************************************************)

Definition and properties of square root on Z

The following tactic replaces all instances of (POS (xI ...)) by `2*(POS ...)+1`, but only when ... is not made only with xO, XI, or xH.

Ltac compute_POS :=
  match goal with
    |  |- context [(Zpos (xI ?X1))] =>
      match constr:(X1) with
	| context [1%positive] => fail 1
	| _ => rewrite (Pos2Z.inj_xI X1)
      end
    |  |- context [(Zpos (xO ?X1))] =>
      match constr:(X1) with
	| context [1%positive] => fail 1
	| _ => rewrite (Pos2Z.inj_xO X1)
      end
  end.

Inductive sqrt_data (n:Z) : Set :=
  c_sqrt : forall s r:Z, n = s * s + r -> 0 <= r <= 2 * s -> sqrt_data n.

Definition sqrtrempos : forall p:positive, sqrt_data (Zpos p).forall p : positive, sqrt_data (Z.pos p)
  refine
    (fix sqrtrempos (p:positive) : sqrt_data (Zpos p) :=
      match p return sqrt_data (Zpos p) with
	| xH => c_sqrt 1 1 0 _ _
	| xO xH => c_sqrt 2 1 1 _ _
	| xI xH => c_sqrt 3 1 2 _ _
	| xO (xO p') =>
          match sqrtrempos p' with
            | c_sqrt _ s' r' Heq Hint =>
              match Z_le_gt_dec (4 * s' + 1) (4 * r') with
		| left Hle =>
                  c_sqrt (Zpos (xO (xO p'))) (2 * s' + 1)
                  (4 * r' - (4 * s' + 1)) _ _
		| right Hgt => c_sqrt (Zpos (xO (xO p'))) (2 * s') (4 * r') _ _
              end
          end
	| xO (xI p') =>
          match sqrtrempos p' with
            | c_sqrt _ s' r' Heq Hint =>
              match Z_le_gt_dec (4 * s' + 1) (4 * r' + 2) with
		| left Hle =>
                  c_sqrt (Zpos (xO (xI p'))) (2 * s' + 1)
                  (4 * r' + 2 - (4 * s' + 1)) _ _
		| right Hgt =>
                  c_sqrt (Zpos (xO (xI p'))) (2 * s') (4 * r' + 2) _ _
              end
          end
	| xI (xO p') =>
          match sqrtrempos p' with
            | c_sqrt _ s' r' Heq Hint =>
              match Z_le_gt_dec (4 * s' + 1) (4 * r' + 1) with
		| left Hle =>
                  c_sqrt (Zpos (xI (xO p'))) (2 * s' + 1)
                  (4 * r' + 1 - (4 * s' + 1)) _ _
		| right Hgt =>
                  c_sqrt (Zpos (xI (xO p'))) (2 * s') (4 * r' + 1) _ _
              end
          end
	| xI (xI p') =>
          match sqrtrempos p' with
            | c_sqrt _ s' r' Heq Hint =>
              match Z_le_gt_dec (4 * s' + 1) (4 * r' + 3) with
		| left Hle =>
                  c_sqrt (Zpos (xI (xI p'))) (2 * s' + 1)
                  (4 * r' + 3 - (4 * s' + 1)) _ _
            | right Hgt =>
                c_sqrt (Zpos (xI (xI p'))) (2 * s') (4 * r' + 3) _ _
            end
        end
    end); clear sqrtrempos; repeat compute_POS;
 try (try rewrite Heq; ring); try omega.
Defined.

Define with integer input, but with a strong (readable) specification.

Definition Zsqrt :
  forall x:Z,
    0 <= x ->
    {s : Z &  {r : Z | x = s * s + r /\ s * s <= x < (s + 1) * (s + 1)}}.forall x : Z,
0 <= x ->
{s : Z &
{r : Z |
x = s * s + r /\ s * s <= x < (s + 1) * (s + 1)}}
  refine
    (fun x =>
      match
	x
	return
        0 <= x ->
        {s : Z &  {r : Z | x = s * s + r /\ s * s <= x < (s + 1) * (s + 1)}}
	with
	| Zpos p =>
          fun h =>
            match sqrtrempos p with
              | c_sqrt _ s r Heq Hint =>
		existT
                (fun s:Z =>
                  {r : Z |
                    Zpos p = s * s + r /\ s * s <= Zpos p < (s + 1) * (s + 1)})
                s
                (exist
                  (fun r:Z =>
                    Zpos p = s * s + r /\
                    s * s <= Zpos p < (s + 1) * (s + 1)) r _)
            end
	| Zneg p =>
          fun h =>
            False_rec
            {s : Z &
              {r : Z |
		Zneg p = s * s + r /\ s * s <= Zneg p < (s + 1) * (s + 1)}}
            (h (eq_refl Datatypes.Gt))
	| Z0 =>
          fun h =>
            existT
            (fun s:Z =>
              {r : Z | 0 = s * s + r /\ s * s <= 0 < (s + 1) * (s + 1)}) 0
            (exist
               (fun r:Z => 0 = 0 * 0 + r /\ 0 * 0 <= 0 < (0 + 1) * (0 + 1)) 0
               _)
    end); try omega.x:Z
p:positive
h:0 <= Z.pos p
s, r:Z
Heq:Z.pos p = s * s + r
Hint:0 <= r <= 2 * s
Z.pos p = s * s + r /\
s * s <= Z.pos p < (s + 1) * (s + 1)
 split; [ omega | rewrite Heq; ring_simplify (s*s) ((s + 1) * (s + 1)); omega ].
Defined.

Define a function of type Z->Z that computes the integer square root, but only for positive numbers, and 0 for others.

Definition Zsqrt_plain (x:Z) : Z :=
  match x with
    | Zpos p =>
      match Zsqrt (Zpos p) (Pos2Z.is_nonneg p) with
	| existT _ s _ => s
      end
    | Zneg p => 0
    | Z0 => 0
  end.

A basic theorem about Zsqrt_plain

Theorem Zsqrt_interval :
  forall n:Z,
    0 <= n ->
    Zsqrt_plain n * Zsqrt_plain n <= n <
    (Zsqrt_plain n + 1) * (Zsqrt_plain n + 1).forall n : Z,
0 <= n ->
Zsqrt_plain n * Zsqrt_plain n <= n <
(Zsqrt_plain n + 1) * (Zsqrt_plain n + 1)
Proof.forall n : Z,
0 <= n ->
Zsqrt_plain n * Zsqrt_plain n <= n <
(Zsqrt_plain n + 1) * (Zsqrt_plain n + 1)
  intros [|p|p] Hp.Hp:0 <= 0
Zsqrt_plain 0 * Zsqrt_plain 0 <= 0 <
(Zsqrt_plain 0 + 1) * (Zsqrt_plain 0 + 1)
p:positive
Hp:0 <= Z.pos p
Zsqrt_plain (Z.pos p) * Zsqrt_plain (Z.pos p) <=
Z.pos p <
(Zsqrt_plain (Z.pos p) + 1) *
(Zsqrt_plain (Z.pos p) + 1)
p:positive
Hp:0 <= Z.neg p
Zsqrt_plain (Z.neg p) * Zsqrt_plain (Z.neg p) <=
Z.neg p <
(Zsqrt_plain (Z.neg p) + 1) *
(Zsqrt_plain (Z.neg p) + 1)
  -Hp:0 <= 0
Zsqrt_plain 0 * Zsqrt_plain 0 <= 0 <
(Zsqrt_plain 0 + 1) * (Zsqrt_plain 0 + 1) now compute.
  -p:positive
Hp:0 <= Z.pos p
Zsqrt_plain (Z.pos p) * Zsqrt_plain (Z.pos p) <=
Z.pos p <
(Zsqrt_plain (Z.pos p) + 1) *
(Zsqrt_plain (Z.pos p) + 1) unfold Zsqrt_plain.p:positive
Hp:0 <= Z.pos p
(let (s, _) := Zsqrt (Z.pos p) (Pos2Z.is_nonneg p) in
 s) *
(let (s, _) := Zsqrt (Z.pos p) (Pos2Z.is_nonneg p) in
 s) <= Z.pos p <
((let (s, _) := Zsqrt (Z.pos p) (Pos2Z.is_nonneg p) in
  s) + 1) *
((let (s, _) := Zsqrt (Z.pos p) (Pos2Z.is_nonneg p) in
  s) + 1)
    now destruct Zsqrt as (s & r & Heq & Hint).
  -p:positive
Hp:0 <= Z.neg p
Zsqrt_plain (Z.neg p) * Zsqrt_plain (Z.neg p) <=
Z.neg p <
(Zsqrt_plain (Z.neg p) + 1) *
(Zsqrt_plain (Z.neg p) + 1) now elim Hp.
Qed.

Positivity

Theorem Zsqrt_plain_is_pos: forall n, 0 <= n ->  0 <= Zsqrt_plain n.forall n : Z, 0 <= n -> 0 <= Zsqrt_plain n
Proof.forall n : Z, 0 <= n -> 0 <= Zsqrt_plain n
  intros n m; case (Zsqrt_interval n); auto with zarith.n:Z
m:0 <= n
Zsqrt_plain n * Zsqrt_plain n <= n ->
n < (Zsqrt_plain n + 1) * (Zsqrt_plain n + 1) ->
0 <= Zsqrt_plain n
  intros H1 H2; case (Z.le_gt_cases 0 (Zsqrt_plain n)); auto.n:Z
m:0 <= n
H1:Zsqrt_plain n * Zsqrt_plain n <= n
H2:n < (Zsqrt_plain n + 1) * (Zsqrt_plain n + 1)
Zsqrt_plain n < 0 -> 0 <= Zsqrt_plain n
  intros H3; contradict H2; auto; apply Z.le_ngt.n:Z
m:0 <= n
H1:Zsqrt_plain n * Zsqrt_plain n <= n
H3:Zsqrt_plain n < 0
(Zsqrt_plain n + 1) * (Zsqrt_plain n + 1) <= n
  apply Z.le_trans with ( 2 := H1 ).n:Z
m:0 <= n
H1:Zsqrt_plain n * Zsqrt_plain n <= n
H3:Zsqrt_plain n < 0
(Zsqrt_plain n + 1) * (Zsqrt_plain n + 1) <=
Zsqrt_plain n * Zsqrt_plain n
  replace ((Zsqrt_plain n + 1) * (Zsqrt_plain n + 1))
     with (Zsqrt_plain n * Zsqrt_plain n + (2 * Zsqrt_plain n + 1));
  auto with zarith.n:Z
m:0 <= n
H1:Zsqrt_plain n * Zsqrt_plain n <= n
H3:Zsqrt_plain n < 0
Zsqrt_plain n * Zsqrt_plain n +
(2 * Zsqrt_plain n + 1) =
(Zsqrt_plain n + 1) * (Zsqrt_plain n + 1)
  ring.
Qed.

Direct correctness on squares.

Theorem Zsqrt_square_id: forall a, 0 <= a ->  Zsqrt_plain (a * a) = a.forall a : Z, 0 <= a -> Zsqrt_plain (a * a) = a
Proof.forall a : Z, 0 <= a -> Zsqrt_plain (a * a) = a
  intros a H.a:Z
H:0 <= a
Zsqrt_plain (a * a) = a
  generalize (Zsqrt_plain_is_pos (a * a)); auto with zarith; intros Haa.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
Zsqrt_plain (a * a) = a
  case (Zsqrt_interval (a * a)); auto with zarith.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <= a * a ->
a * a <
(Zsqrt_plain (a * a) + 1) * (Zsqrt_plain (a * a) + 1) ->
Zsqrt_plain (a * a) = a
  intros H1 H2.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
Zsqrt_plain (a * a) = a
  case (Z.le_gt_cases a (Zsqrt_plain (a * a))); intros H3.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:a <= Zsqrt_plain (a * a)
Zsqrt_plain (a * a) = a
a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:Zsqrt_plain (a * a) < a
Zsqrt_plain (a * a) = a
  -a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:a <= Zsqrt_plain (a * a)
Zsqrt_plain (a * a) = a Z.le_elim H3; auto.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:a < Zsqrt_plain (a * a)
Zsqrt_plain (a * a) = a
    contradict H1; auto; apply Z.lt_nge; auto with zarith.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:a < Zsqrt_plain (a * a)
a * a < Zsqrt_plain (a * a) * Zsqrt_plain (a * a)
    apply Z.le_lt_trans with (a * Zsqrt_plain (a * a)); auto with zarith.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:a < Zsqrt_plain (a * a)
a * Zsqrt_plain (a * a) <
Zsqrt_plain (a * a) * Zsqrt_plain (a * a)
    apply Z.mul_lt_mono_pos_r; auto with zarith.
  -a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H2:a * a <
(Zsqrt_plain (a * a) + 1) *
(Zsqrt_plain (a * a) + 1)
H3:Zsqrt_plain (a * a) < a
Zsqrt_plain (a * a) = a contradict H2; auto; apply Z.le_ngt; auto with zarith.a:Z
H:0 <= a
Haa:0 <= a * a -> 0 <= Zsqrt_plain (a * a)
H1:Zsqrt_plain (a * a) * Zsqrt_plain (a * a) <=
a * a
H3:Zsqrt_plain (a * a) < a
(Zsqrt_plain (a * a) + 1) * (Zsqrt_plain (a * a) + 1) <=
a * a
    apply Z.mul_le_mono_nonneg; auto with zarith.
Qed.

Zsqrt_plain is increasing

Theorem Zsqrt_le:
 forall p q, 0 <= p <= q  ->  Zsqrt_plain p <= Zsqrt_plain q.forall p q : Z,
0 <= p <= q -> Zsqrt_plain p <= Zsqrt_plain q
Proof.forall p q : Z,
0 <= p <= q -> Zsqrt_plain p <= Zsqrt_plain q
  intros p q [H1 H2].p, q:Z
H1:0 <= p
H2:p <= q
Zsqrt_plain p <= Zsqrt_plain q
  Z.le_elim H2; [ | subst q; auto with zarith].p, q:Z
H1:0 <= p
H2:p < q
Zsqrt_plain p <= Zsqrt_plain q
  case (Z.le_gt_cases (Zsqrt_plain p) (Zsqrt_plain q)); auto; intros H3.p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Zsqrt_plain p <= Zsqrt_plain q
  assert (Hp: (0 <= Zsqrt_plain q)).p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
0 <= Zsqrt_plain q
p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
Zsqrt_plain p <= Zsqrt_plain q
  {p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
0 <= Zsqrt_plain q apply Zsqrt_plain_is_pos; auto with zarith. }p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
Zsqrt_plain p <= Zsqrt_plain q
  absurd (q <= p); auto with zarith.p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
q <= p
  apply Z.le_trans with ((Zsqrt_plain q + 1) * (Zsqrt_plain q + 1)).p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
q <= (Zsqrt_plain q + 1) * (Zsqrt_plain q + 1)
p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
(Zsqrt_plain q + 1) * (Zsqrt_plain q + 1) <= p
  case (Zsqrt_interval q); auto with zarith.p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
(Zsqrt_plain q + 1) * (Zsqrt_plain q + 1) <= p
  apply Z.le_trans with (Zsqrt_plain p * Zsqrt_plain p); auto with zarith.p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
(Zsqrt_plain q + 1) * (Zsqrt_plain q + 1) <=
Zsqrt_plain p * Zsqrt_plain p
p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
Zsqrt_plain p * Zsqrt_plain p <= p
  apply Z.mul_le_mono_nonneg; auto with zarith.p, q:Z
H1:0 <= p
H2:p < q
H3:Zsqrt_plain q < Zsqrt_plain p
Hp:0 <= Zsqrt_plain q
Zsqrt_plain p * Zsqrt_plain p <= p
  case (Zsqrt_interval p); auto with zarith.
Qed.

Equivalence between Zsqrt_plain and Z.sqrt

Lemma Zsqrt_equiv : forall n, Zsqrt_plain n = Z.sqrt n.forall n : Z, Zsqrt_plain n = Z.sqrt n
Proof.forall n : Z, Zsqrt_plain n = Z.sqrt n
 intros.n:Z
Zsqrt_plain n = Z.sqrt n destruct (Z_le_gt_dec 0 n).n:Z
l:0 <= n
Zsqrt_plain n = Z.sqrt n
n:Z
g:0 > n
Zsqrt_plain n = Z.sqrt n
 symmetry.n:Z
l:0 <= n
Z.sqrt n = Zsqrt_plain n
n:Z
g:0 > n
Zsqrt_plain n = Z.sqrt n apply Z.sqrt_unique; trivial.n:Z
l:0 <= n
Zsqrt_plain n * Zsqrt_plain n <= n <
Z.succ (Zsqrt_plain n) * Z.succ (Zsqrt_plain n)
n:Z
g:0 > n
Zsqrt_plain n = Z.sqrt n
 now apply Zsqrt_interval.n:Z
g:0 > n
Zsqrt_plain n = Z.sqrt n
 now destruct n.
Qed.