Square Root

Square root is implement using Newton's method.

Section SquareRoot.

Variable a : Q.

Hypothesis Ha : 1 <= a <= 4.
Square root is first defined on [1,4]. Iterating this formula will converge to the square root of a.
Definition root_step (b:Q) : Q := b/2 + a/(2*b).

Definition root_has_error (e:Qpos) (b:Q) := a <= (b+e)^2 /\ (b-e)^2 <= a.
Lemma root_has_error_le : forall (e1 e2:Qpos) (b:Q), e2 <= b -> e1 <= e2 -> root_has_error e1 b -> root_has_error e2 b.

Lemma root_error_bnd : forall (e:Qpos) b, e <= 1 -> 1 <= b -> (root_has_error e b) -> b <= 2 + e.

Lemma root_has_error_ball : forall (e1 e2:Qpos) (b1 b2:Q), (e1 + e2<=1) -> (1 <= b1) -> (1 <= b2) -> root_has_error e1 b1 -> Qball e2 b1 b2 -> root_has_error (e1+e2) b2.

Lemma ball_root_has_error : forall (e1 e2:Qpos) (b1 b2:Q), ((e1 + e2)<=1) -> (1<=b1) -> (1<=b2) -> root_has_error e1 b1 -> root_has_error e2 b2 -> Qball (e1+e2) b1 b2.

Lemma root_step_error : forall b (e:Qpos), (1 <= b) -> (e <= 1) -> root_has_error e b -> root_has_error ((1#2)*(e*e)) (root_step b).

Our initial estimate is (a+1)/2 with an error of 1/2
Definition initial_root :Q := ((1#2)*(a+1)).

Lemma initial_root_error : root_has_error (1#2) initial_root.

Lemma root_step_one_le : forall b, (1 <= b)-> (1 <= root_step b).

Lemma initial_root_one_le : (1 <= initial_root).

Each step squares the error
Fixpoint root_loop (e:Qpos) (n:nat) (b:Q) (err:positive) {struct n} : Q :=
match n with
| O => b
| S n' => if (Qlt_le_dec_fast e (1#err))
          then let err':= (err*err)%positive in root_loop e n' (approximateQ (root_step b) (2*err')) err'
          else b
end.


Lemma root_loop_one_le : forall e n b err, (1 <= b)-> (1 <= root_loop e n b err).

Lemma root_loop_error : forall (e:Qpos) n b err, (1 <= b) -> root_has_error (1#err) b -> (1#(iter_nat n _ (fun x => (x * x)%positive) err))<=e ->
 root_has_error (Qpos_min (1 # err) e) (root_loop e n b err).

Find a bound on the number of iterations we need to take.
Lemma root_max_steps : forall (n d:positive), (1#(iter_nat (S (Psize d)) _ (fun x => (x * x)%positive) 2%positive))<=(n#d)%Qpos.

Square root on [1,4].
Definition sqrt_raw (e:QposInf) : Q :=
match e with
| ∞ => 1
| Qpos2QposInf e' => root_loop e' (S (Psize (Qden e'))) initial_root 2
end.

Lemma sqrt_regular : is_RegularFunction sqrt_raw.

Definition rational_sqrt_mid : CR := Build_RegularFunction sqrt_regular.

Lemma rational_sqrt_mid_err : forall (e:Qpos), (e <= 1) -> root_has_error e (approximate rational_sqrt_mid e).

Lemma rational_sqrt_mid_one_le : forall (e:QposInf), 1 <= (approximate rational_sqrt_mid e).

Lemma rational_sqrt_mid_le_3 : forall (e:QposInf), (approximate rational_sqrt_mid e) <= 3.


Lemma rational_sqrt_mid_correct0 : (CRpower_positive 2 rational_sqrt_mid == ' a)%CR.

Lemma rational_sqrt_mid_correct1 : ('0 <= rational_sqrt_mid)%CR.

Lemma rational_sqrt_mid_correct : forall H, (rational_sqrt_mid == IRasCR (sqrt (inj_Q IR a) H))%CR.

End SquareRoot.

By scaling the input the range of square root can be extened upto 4^n.
Definition rational_sqrt_big_bounded (n:nat) a (Ha:1 <= a <= (4^n)%Z) : CR.
Defined.

Lemma rational_sqrt_big_bounded_correct : forall n a Ha H,
 (@rational_sqrt_big_bounded n a Ha == IRasCR (sqrt (inj_Q IR a) H))%CR.

By scaling the other direction we can extend the range down to 4^(-n).
Definition rational_sqrt_small_bounded (n:nat) a (Ha:/(4^n)%Z <= a <= 4) : CR.
Defined.

Lemma rational_sqrt_small_bounded_correct : forall n a Ha H,
 (@rational_sqrt_small_bounded n a Ha == IRasCR (sqrt (inj_Q IR a) H))%CR.

And hence it is defined for all postive numbers.
Definition rational_sqrt_pos a (Ha:0<a) : CR.
Defined.

Lemma rational_sqrt_pos_correct : forall a Ha H,
 (@rational_sqrt_pos a Ha == IRasCR (sqrt (inj_Q IR a) H))%CR.

We define the square root for all nonpositive numbers.
Definition rational_sqrt a : CR :=
match (Qlt_le_dec_fast 0 a) with
|left H => rational_sqrt_pos a H
|right _ => (' 0)%CR
end.

Lemma rational_sqrt_correct : forall a H,
 (@rational_sqrt a == IRasCR (sqrt (inj_Q IR a) H))%CR.

Square root is uniformly continuous everywhere.
Definition sqrt_modulus (e:Qpos) : Q+∞ := Qpos2QposInf (e*e).

Lemma sqrt_uc_prf : is_UniformlyContinuousFunction rational_sqrt sqrt_modulus.


Definition sqrt_uc : Q_as_MetricSpace --> CR :=
Build_UniformlyContinuousFunction sqrt_uc_prf.

Definition CRsqrt : CR --> CR := Cbind QPrelengthSpace sqrt_uc.

Lemma CRsqrt_correct : forall x H,
 (IRasCR (sqrt x H) == CRsqrt (IRasCR x))%CR.