Limits with Signed Digit Streams
Franziskus Wiesnet

FriesBScSdmult
BergerMiyamotoSchwichtenbergSeisenberger11
CiaffaglineGianantonio, Gianantonio, Tsuiki2002
Correcting symmetry imperfections 
in  linear multipole traps
    M. Knoop
    December 30, 2023
=============================================================


We work with the signed digit representation of abstract real numbers, which 
roughly is the binary representation enriched by the additional digit -1. 
The main objective of this paper is an algorithm which takes a sequence of signed 
digit representations of reals and returns the signed digit representation of their limit, if the sequence converges. 
As a first application we use this algorithm together with Heron's method to build up an algorithm which 
converts the signed digit representation of a non-negative real number into the signed digit 
representation of its square root. 

Instead of writing the algorithms first and proving their correctness afterwards, 
we work the other way round, in the tradition of program extraction from proofs. 
In fact we first give constructive proofs, and from these proofs we then compute the extracted terms, which is the desired algorithm. 
The correctness of the extracted term follows directly by the Soundness Theorem 
of program extraction. In order to get the extracted term from some proofs which are
often quite long, we use the proof assistant Minlog. However, to apply the extracted terms, the programming language Haskell is useful. Therefore after each proof we show a 
notation of the extracted term, which can be easily rewritten as a definition in Haskell.

Keywords: signed digit code, real number computation, coinductive definitions, corecursion, program extraction, realizability, convergence of reals, square root, Heron's method



§ INTRODUCTION AND MOTIVATION

There are several ways to define constructive real numbers. One of the best-known methods is to define them as Cauchy sequences of rational numbers with a Cauchy modulus. In our case, rather than we are not interested in a specific definition of real numbers, we are interested in their signed digit representation (SD code). Therefore, every quantifier on reals is non-computational, i.e. we write ∀_x^nc and ∃_x^r. Theses decorations just mean that the bounded variable does not appear in the computational content of the proof. Logically one can ignore the decorations. By x and y we denote variables of reals. Instead of using the concrete real computationally, we use their signed digit representation in the extracted term.

Attention should be paid to the equality between reals numbers. By equality = between two reals we mean the “real equality”, which is an equivalence relation and compatible with the usual operators and relations on the reals. The specific definition of the real equality depends on the definition of the real numbers. Generally, the real equality in not the same as the Leibniz equality. For instance, if one defines real numbers as Cauchy sequences with modulus, the real numbers ⟨(2^-n)_n∈ℕ,id ⟩ and ⟨(0)_n∈ℕ,id ⟩   are equal w.r.t. the real equality but they are not Leibniz equal.

One of the first paper <cit.> where  signed digits are used to represent real numbers, was published by Edwin Wiedmer in 1980.
The SD code of reals is similar to the binary code of reals but in addition to the digits 0 and 1, the SD code has the digit -1, which we also denote by 1.
Since every real x can be represented as x=k+x', where k is an integer and -1≤ x'≤ 1, we work on the interval [-1,1].

A real number x in [-1,1] has a binary representation if it can be written as

    x = s∑_i=1^∞a_i2^-i,
 
where s∈{-1,1} and a_i∈{0,1} for every i. 

If one reads one by one the binary representation of a concrete real number, in each step the interval in which the real number is located, is halved. Thus from the binary code one can determine the real number arbitrarily exactly.

On the other hand it is not always possible to compute the binary representation of a given real number or even to compute the binary representation of x+y/2 out of the binary representations of x and y. Hereby “compute” means to get an algorithm which takes the binary streams of x and y and gives back the digits of the binary stream x+y/2 one by one in finitely many steps for each digit. That is, to compute finitely many binary digits of x+y/2 can only use finitely many binary digits of x and y. This is not possible due to the “gaps” in the binary representation. They are illustrated in Figure <ref> at 0, 1/2, - 1/2, 1/4 and so on. Form the first digit (i.e. + or -) of a real x, one can decide 0≤ x or x≤ 0, with in general can not be done for constructive reals.

The signed digit code fills these gabs. For a real number x∈ [-1,1] it is given by

    x=∑_i=0^∞d_i2^-i,

where d_i∈{1,0,1} for every i.

As the illustration in Figure <ref> shows, to compute the first signed digit of a real number x∈ [-1,1] one has to decide which of the cases x≤ 0, -1/2≤ 0 ≤1/2 or 0≤ x holds. That this is possible, follows from the Comparability Theorem, which says that for reals x, y and z with x<y one has z≤ y ∨ x≤ y. 

Figure <ref> also shows that the SD code of a real number except -1 and 1 is not unique, whereas the binary code is “almost” everywhere unique.

A stream of signed digits is an infinite list d_1d_2d_3… of elements in 
    Sd := {1,0,1}.
 We will not prove something about signed digit streams directly, but we use the coinductively define predicate ^coI, which is given in the next section. For a real number x a realiser of ^coIx is a signed digit stream of x. With the Soundness Theorem of program extraction proven in <cit.> the proofs of the corresponding theorems for the signed digit streams are obtained.
 The idea to use coinductive algorithms to describe the operators on the reals goes back to Alberto Ciaffaglione and Pietro Di Gianantonio <cit.>. The idea to use coinductively defined predicates and the Soundness Theorem in this context is due to Ulrich Berger and Monika Seisenberger <cit.>. The notation and definitions in this paper are taken from Kenji Miyamoto and Helmut Schwichtenberg <cit.>.

For computing the extraced terms and verifying the correctness of the proofs, we have used the proof assistant Minlog <cit.> in some cases. After each proof we  state its computational content not in the notation of Minlog but in the notation of Haskell, since the runtime of the programs in Haskell is shorter, and one can define the terms in a more readable way. 
 

We now proceed as follows:

In Section <ref> we give a definition of ^coI and state the two axioms ^coI^- and ^coI^+ for this predicate. The two lemmas <ref> and <ref> are often used in the proofs below. The main result of this work is Theorem <ref>, the Convergence Theorem, at the end of Section <ref>. An application of this theorem is shown in Section <ref>. Here we use Heron's method together with the Convergence Theorem to get a signed digit representation of the square root of a non-negative real number from its  signed digit representation. Section <ref> is about potential extensions of this work.


§ FORMALISATION

In this section we use the theory of coinductively defined predicates given in <cit.> to formalise the statement that a real x has an SD representation.

We define ^coI as the greaterst fixpoint of 
    Φ(X):= {x|∃_d,x'^r(Sd d∧ Xx'∧ |x|≤ 1∧ x=d+x'/2)}

Here ∃^r_xAx := μ_X(∀_x^nc(Ax→ X)) is the existential quantifier where the quantified variable x does not appear in the computational content.

In a short form we have ^coI:=ν_X(Φ(X)). Therefore a realiser of ^coIx has the type  
    τ(^coI)=μ_τ(X)(τ(Φ(X))→τ(X))=μ_ξ(Sd→ξ→ξ ).
 Here we have identified τ(Sd)=μ_ξ(ξ,ξ,ξ) with Sd itself. We define :=τ(^coI) and with  we denote the only constructor of . Then in Haskell notation  is given by 

    ::=  Sd .

In this notation one easily sees that an element d v can be interpreted as an SD code with the first digit d and tail v. Sometimes we abbreviate d v by only writing dv.

From the definition of ^coI as greatest fixpoint of Φ we get two axioms for ^coI. 

    ^coI^-:    ^coI⊆Φ(^coI)
    ^coI^+:    X⊆Φ(^coI∪ X)→ X⊆^coI

where X is an unary predicate variable.
The first axiom ^coI^- says that ^coI is a fixpoint of Φ and for obvious reasons it is called the elimination rule. Expressed in elementary formulas it is given by

    ^coI^-: ∀_x^nc.^coIx→∃^r_d,x'(Sd d∧^coIx'∧ |x|≤ 1∧ x=d+x'/2).

The type of this axiom is τ(^coI^-)=→Sd× and a realiser of it is given by the destructor 𝒟 with computation rule

    𝒟   : →Sd×,
    𝒟(dv)    := ⟨ d,v⟩.

The destructor takes a stream and returns a pair consisting of its first digit and its tail. With the projectors π_0 and π_1 one gets the first digit and the tail, respectively.

The second axiom ^coI^+ is called the introduction axiom of ^coI and says that ^coI is the greatest fixpoint in the a strong sense. We use the following long version:

    ^coI^+: ∀_x^nc.   Xx→
       ∀_x^nc(Xx →∃_d,x'^r(Sd d∧(^coI∪ X)x'∧|x|≤ 1∧ x=d+x'2))
       →^coIx

The type of this axiom depends on the type of the predicate variable X:

    τ(^coI^+)= τ(X) →(τ(X)→Sd×(+τ(X)))→

A realiser of ^coI^+ is the corecursion operator ^coℛ. It is given by the computation rule

    ^coℛ   :τ(X) →(τ(X)→Sd×(+τ(X)))→,
    
    ^coℛt f    :=  (π_0(ft))[,λ_t'^coℛt' f]π_1(ft),

where [F_0,F_1]in_iT:=F_iT for i∈{0,1}. Here in_0 and in_1 are the two constructors of the type sum +τ(X).
If π_1(ft) has the form in_0v, the corecursion stops and we have (π_0(ft))v as signed digit representation. If it has the form in_1t', the corecursion goes further with the new argument t'. In both cases we have obtained at least the first digit of the stream. With iteration of the corecursion, if necessary, we generate each digit of the stream one by one.

Now we prove two lemmas, which are often used  in our proofs:

The predicate ^coI is compatible with the real equality, i.e.

    ∀_x,y^nc.^coIx→ x=y →^coIy.



We apply ^coI^+ to the predicate Px:=∃_y^r(^coIy∧ x=y):

    ∀_x^nc.   ∃_y^r(^coIy∧ x=y)→
       ∀_x^nc(Px →∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2))
       →^coIx

As one sees, it is sufficient to prove the second premise. Therefore we show

    ∀_x^nc(∃_y^r(^coIy∧ x=y) →∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2)).

Let x, y with ^coIy and x=y be given. From ^coIy we get e∈Sd and y'∈^coI with |y|≤ 1 and y=e+y'/2. We now define d:=e and x':=y' and our goal follows directly because of x=y.

In the following proofs this theorem is used tacitly. If one have a model of reals in which the real equality is the same as the Leibniz equality, one do not need this theorem.
The extracted term of it is given by

    λ_u.^coℛu(λ_v⟨π_0(𝒟v),in_0(π_1(𝒟v))⟩).

For stream of the form du this term is the identity function. Since we always deal only with such streams, we drop this term hereafter.


    ∀_x,d^nc.Sd d →^coIx →^coId+x/2



We apply ^coI^+ to the predicate 
    Px := ∃^r_d,x'(Sd d ∧^coIx' ∧ x=d+x'/2).

This leads to the formula

    ∀_x^nc.   ∃^r_d,x'(Sd d ∧^coIx' ∧ x=d+x'/2)→
       ∀_x^nc(Px →∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2))
       →^coIx.

In order to prove the goal formula, it is sufficient to prove the second premise. Therefore the new goal formula is

    ∀_x^nc.    ∃^r_d,x'(Sd d ∧^coIx' ∧ x=d+x'/2) →
       ∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2).

This is almost a tautology. The only part one has to consider is |x|≤ 1. From Sd d and ^coIx' we get |d|,|x'|≤ 1 and therefore |x|=|d+x'/2| ≤ 1. 

Formally the extracted term of this proof is given by

    λ_dλ_u. ^coℛuλ_v⟨ d,in_0v⟩.

If we use the computation rule of ^coℛ once, we get λ_dλ_u. du, which is identified with  itself. Therefore the constructor  is actually the computational content of this lemma. 


§ CONVERGENCE THEOREM

The aim of this chapter is to prove the convergence theorem in the SD code. It says that the limit of each convergent sequence in ^coI is also in ^coI. As extracted term we expect a function which takes a sequence of signed digit streams (i.e. a term of type ℕ→) and its modulus of convergence and returns a signed digit stream.

We do the proof step by step and therefore we first prove a few lemmas:


    ∀^nc_x.^coIx → x≤ 0 →^coI(x+1)
    ∀^nc_x.^coIx → 0≤ x →^coI(x-1)



Since both formulas are shown analogously we only show the first formula. Therefore we use the introduction axiom of ^coI, which is given by

    ∀_x^nc.   Px→
       ∀_x^nc(Px →∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2))
       →^coIx.

For Px we insert ∃_y( ^coIy∧ y≤ 0 ∧ y+1 = x ). Hence it is sufficient to prove the second premise. Let x, y, ^coIy, y≤ 0 and y+1=x be given. Our goal is 

    ∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2).

From ^coIy we get e and y' with Sd e, ^coIy', |y|≤ 1 and y=e+y'/2. Independent on d and x' we always get |x|≤ 1 out of |y|≤ 1, y≤ 0 and x=y+1.
In order to prove the remaining part of the formula, we do case distinction on Sd e:

If e = -1, we define d:=1 and x':=y'. Here Sd d and ^coIx' follow directly and we also have 
    x=y+1=-1+y'/2+1 = 1+y'/2=d+x'/2.

If e=0, we define d:=1 and x':=y'+1. We directly have Sd d. In this case we prove Px'. Hence we show ^coIy', y'≤ 0 and x'=y'+1. ^coIy' and x'=y'+1 are already given and y'≤ 0 follows directly from y≤ 0 and y=0+y'/2.

The last case is e=1. Because of y≤ 0, y=-1+y'/2 and |y'|≤ 1 it is only possible that y is equal to 0 and therefore x=1. Hence we define d:=1 and x':=1. Then Sd d and x=d+x'/2 are obviously true and ^coI 1 is true because 1 has the SD representation 111….

A realiser of the first formula in this lemma is a function , which takes a signed digit stream of a real number x and returns a signed digit stream of x+1 if x≤ 0. 

Using the formal definition of the extracted term we get for  the term

    = λ_u.^coℛu(λ_v.[⟨ 1,in_0(π_1(𝒟v))⟩,⟨ 1,in_1(π_1(𝒟v))⟩,⟨ 1,in_01⃗⟩]π_0(𝒟v)),

where [F_1,F_0,F_1]d:=F_d for d∈Sd and 
1⃗ is the infinite list with each entry equal to 1. 
Another way to characterise this function  is to give its computation rules:

    (1v)    := 1v 
     (0v)    := 1(v)  
     (1v)    := [1,1,…]

Analogously as extracted term of the second statement of this lemma, we get a function :→ which is characterised by the rules

    (1v)    := [1,1,…] 
     (0v)    := 1(v)  
     (1v)    := 1v.

It takes a signed digit stream of a real x and returns a signed digit stream of x-1 if 0≤ x.

With this lemma we can now prove the following lemma easily:


    ∀^nc_x.^coIx → |x|≤1/2→^coI(2x)



From ^coIx we get d, x' with Sd d, ^coIx', |x|≤ 1 and x=d+x'/2. Case distinction on Sd d gives three cases:

d=1: Here 2x=1+x' and with |x|≤1/2 it follows x'≤ 0. Therefore the first formula of Lemma <ref> gives ^coI(2x).

d=-1: The proof in this case is done analogously but we need the second statement of Lemma <ref>.

d=0: This leads to 2x=x' and ^coIx' is already given.

If we keep the definition of the functions  and  form the computational context of the privios lemma, the computational content :→ of this lemma is given by

    = λ_u.[(π_1(𝒟u)),π_1(𝒟u),(π_1(𝒟u))]π_0(𝒟u).

Again we give a more readable characterisation of  by the computation rules

    (1u)   :=u
    (0u)   :=u
    (1u)   :=u.



    ∀^nc_x.^coIx  →^coI( x/2±1/4)



^coIx gives x'∈^coI and d ∈Sd with x=d+x'/2. We show only ^coI( x/2 + 1/4) because the proof of ^coI( x/2 - 1/4) is done analogously. Case distinction on Sd d leads to the following three cases: 

d=1 gives x/2 +1/4 = 2+x'/4 = 1+x'/2/2, 

d=0 gives x/2 +1/4 = 1+x'/4 = 1+x'/2/2 and

d=-1 gives x/2 +1/4=x'/4=x'/2/2.

In each case we get ^coI(x/2 +1/4) by using Lemma <ref> twice.

We denote the extracted term of the proven statement by q^+. From the proof and the fact, that the extracted term of Lemma <ref> is given by , one easily sees that q^+ has the following computation rules:

    q^+(1u)   := 0 0 u
    
    q^+(0u)   := 0 1 u
    
    q^+(1u)   := 1 0 u

Analogously, the computational content q^- of the statement ∀^nc_x.^coIx  →^coI( x/2 - 1/4) is characterised by

    q^-(1u)   := 1 0 u
    
    q^-(0u)   := 0 1 u
    
    q^-(1u)   := 0 0 u.


Before we prove the Convergence Theorem, we have to give a definition of convergence. In this definition the witness of convergence is included:

Let a real x, a sequence f:ℕ→ℝ of reals and a modulus M:ℤ^+→ℕ be given. We say f convergences to x with modulus M if

    ∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p.



Let f:ℕ→ℝ be a sequence of reals in ^coI which converges to a real x with a modulus M, then also x is in ^coI.

Expressed in a formula:

    ∀^nc_x∀^nc_f∀_M.∀_n^coIf(n)∧∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p→^coIx



We show the equivalent formula

    ∀^nc_x.∃^r_f∃^d_M(∀_n^coIf(n)∧∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p)→^coIx.

The existence quantifier ∃^d is an inductively defined predicate and given by ∃^d_xA(x) := μ_X(∀_x(A(x)→∃_x^dA(x))). This means that in contrast to ∃^r the quantified variable is also part of the computational content.

To prove this formula we apply ^coI^+ to the predicate 

    Px := ∃^r_f∃^d_M(∀_n^coIf(n)∧∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p)
 
and show the second premise, which is given by

    ∀_x^nc.   ∃^r_f∃^d_M(∀_n^coIf(n)∧∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p) →
       ∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2).

Let x, f, M, ∀_n^coIf(n) and ∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p be given. We show

    ∃_d,x'^r.   Sd d∧∃^r_g∃^d_N(∀_n^coIg(n)∧∀_p∀_n≥ N(p)|g(n)-x'|≤ 2^-p)∧
       |x|≤ 1∧ x=d+x'2.

Regardless of the choice of d and x' we get ∀_n |f(n)|≤ 1 from ∀_n^coIf(n) and with ∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p it follows |x|≤ 1. Therefore in each case we consider |x|≤ 1 as proven.

Specializing ∀_n^coIf(n) to M(4) leads to ^coIf(M(4)) and triple application of ^coI^- gives d_1,d_2,d_3∈Sd and y'∈^coI such that f(M(4))=4d_1+2d_2+d_3+y'/8 or in short notation 
    f(M(4))=d_1d_2d_3y'.

Now we do case distinction on this representation of f(M(4)).

If f(M(4)) has one of the forms 11d_3y', 10d_3y', 111y', 110y' 011y' or 010y', it follows that 1/8≤ f(M(4)). In this case we choose 

    d:=1 and x':=2x-1,
 
then Sd d and x=d+x'/2 follow directly. Furthermore we define 
    g(n):= 2f(M(4)∨ n)-1
 for all n∈ℕ, where m∨ l := max{m,l}, and 
    N(p):=M(p+1)
 for all p∈ℤ^+. The formular ∀_p∀_n≥ N(p)|g(n)-x'|≤ 2^-p is a direct consequence of ∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p and it remains to show ∀_n^coIg(n). We calculate 
    g(n)=2f(M(4)∨ n)-1 =4( f(M(4)∨ n)/2 - 1/4).
 Lemma <ref> gives ^coI( f(M(4)∨ n)/2 - 1/4). Futhermore we have f(M(4))≥1/8 and therefore 
    f(M(4)∨ n)    = (f(M(4)∨ n)-x)+(x -f(M(4)))+f(M(4)) 
       ≥ -1/16 -1/16 + 1/8 =0.

Thus 0 ≤ f(M(4)∨ n)≤ 1 implies |f(M(4)∨ n)/2 - 1/4| ≤1/4 and Lemma <ref> applied twice gives ^coIg(n). 

If f(M(4)) has one of the forms 11d_3y', 10d_3y', 111y', 110y' 011y' or 010y', it follows f(M(4))≤ -1/8. Here we define 
    d:=-1, x':=2x+1,g:=λ_n (2f(M(4)∨ n)+1) and N:=λ_pM(p+1).
 The proof in this case is analogous to the proof of the first case. 

It remains to consider the case that f(M(4)) has one of the forms 00d_3y', 111y', 111y', 011y' or 011y'. Here we have -1/4≤ f(M(4)) ≤1/4 and we define 
    d:=0 and  x':=2x.

The formulars Sd d and x=d+x'/2 are obvious. In order to prove 
    ∃^r_g∃^d_N(∀_n^coIg(n)∧∀_p∀_n≥ N(p)|g(n)-x'|≤ 2^-p),
 we define 
    g:=λ_n 2f(M(4)∨ n) and N:=λ_pM(p+1).
 The second part of the conjunction follows from ∀_p∀_n≥ M(p)|f(n)-x|≤ 2^-p, which is given. And because of 
    |f(M(4)∨ n)|   ≤ |f(M(4)∨ n)-x| + |x - f(M(4))| +|f(M(4))|
       ≤1/16 + 1/16 +1/4≤1/2
 for every n, we have ∀_n ^coIg(n) by Lemma <ref>.

We denote the extracted term by . It has the type 
    :(ℤ^+→ℕ)→ (ℕ→)→.

It takes the modulus of convergence and the sequence of streams and returns the stream of the limit value. In order to give an at most readable characterisation of , we define the following sets:

    R   := {11v, 10v, 111v, 110v, 011v, 010v| v∈}
    M   :={00v, 111v, 111v, 011v, 011v|v∈}
    L   := {11v, 10v, 111v, 110v,011v, 010v|v∈}

These three sets correspond to the three cases in the proof and therefore we have the following rule for :

    M F :=
    
         1 ( λ_pM(p+1) λ_n(q^-F(M(4)∨ n)))   if  F(M(4))∈R
    
        0 ( λ_pM(p+1) λ_n(F(M(4)∨ n)))   if  F(M(4))∈M
    1 ( λ_pM(p+1) λ_n(q^+F(M(4)∨ n)))   if  F(M(4))∈L

Because of readability, we have omitted some brackets, for example F(M(4)∨ n)) shall be read as (F(M(4)∨ n))).
The functions , q^+ and q^- are the computational content of the lemmas above and, as in the proof, ∨ between natural numbers means their maximum.

One note that the definition of the new sequence is not unique. For reasons of efficiency one should be flexible with the choice of the new sequence, which is called g in the proof above. For example by choosing g one can replace M(4) ∨ n by M(4)+n. The efficiency depends on the concrete sequence. In the next section we define the Heron sequence and apply it to . In this case the definition with M(4) ∨ n is most efficient.

Since the definition of this function is not so easy, we give an example for an implementation as a Haskell program. First the sets R and L are realised as Boolean functions:

Then one can define  by case distinction:

In this implementation the constructor  is denoted by  and written as an infix. The tree elements in Sd are given by ,  and  and shall be interpreted by 1, 0 and 1, respectively.


§ SQUARE ROOT


We define H:ℝ→ℕ→ℝ by the computation rules

    H(x,0)    := 1
    H(x,n+1)   := 1/2( H(x,n)+x/H(x,n)).

For every non-negative x the sequence λ_nH(x,n) =: H(x):ℕ→ℝ is the sequence, we get from Heron's method with start value 1.

Please note that H is well-defined for non-negative x since one can easily prove H(x,n)≥ 2^-n by induction.

For x∈ [0,1] we have that H(x) converges to √(x) with modulus ι:ℤ^+→ℕ which is the inclusion from ℤ^+ to ℕ. Furthermore ∀_n√(x)≤ H(x,n) holds.


Let x∈ [0,1] be given. For each n∈ℕ  we define Δ(x,n) := H(x,n)-√(x). A short calculation gives

    Δ(x,n+1)   =1/2( H(x,n)+x/H(x,n)) -√(x) = (H(n,x))^2 -2H(x,n)√(x) +x/2H(x,n)
       = (Δ(x,n))^2/2H(x,n).

By induction on n one easily gets 0≤ H(x,n) and therefore 0≤Δ(x,n+1). Because of  Δ(x,0)=1-√(x)≥ 0 we have ∀_n √(x)≤ H(x,n).

Furthermore we calculate as follows:

    Δ(x,n+1)   = (Δ(x,n))^2/2H(x,n)=1/2Δ(x,n)Δ(x,n)/H(x,n)=1/2Δ(x,n)(1-√(x)/H(x,n))
       ≤1/2Δ(x,n)

Therefore by induction we have |H(x,n)-√(x)|=Δ(x,n)≤ 2^-n and this implies ∀_p∀_n≥ p|H(x,n)-√(x)|≤ 2^-p i.e. H(x) converges to √(x) with modulus ι.

This lemma does not have any computational content, but it says that ι is a witness for the convergence of H x to √(x) from above. In some special cases we can even improve the modulus:

If x∈[1/4,1] we can also choose :ℤ^+→ℕ as modulus of the converges from H(x) to √(x). For a positive integer p we define (p) as the least natural number n with p≤ 2^n.


Let x∈ [1/4,1] be given. From Lemma <ref> we know ∀_n√(x)≤ H(x,n) and therefore ∀_n1/2≤ H(x,n). In the proof of Lemma <ref> the formula

    Δ(x,n+1)= (Δ(x,n))^2/2H(x,n)

is proven. In total it follows Δ(x,n+1)≤ (Δ(x,n))^2. Because of 1/2≤√(x) we get by induction

    Δ(x,n)≤ 2^-2^n

for each natural number n. Hence for given p and n≥(p) we have p≤ 2^n and therefore 

    |H(x,n)-√(x)|=Δ(x,n)≤ 2^-2^n≤2^-p


One possibility, to implement the function , is by defining an auxiliary function :ℤ^+→ℕ→ℕ with the computation rules

    p n := 
    
         n   if  p≤ 2^n 
     p  (n+1)   otherwise

and then setting (p) :=  p 0.

We do not use this lemma and the function  in the main theorem but we use it for concrete calculations to reduce their duration.

For all x ∈^coI with 1/16≤ x we have ∀_n^coI(H(x,n)). Expressed as a formula this means

    ∀_x.^coIx→1/16≤ x →∀_n^coI(H(x,n)).



We use the results  of <cit.> and of Section 3.3 from <cit.>. In both scripts there are the following statements proven:

    ∀_x,y^nc.^coIx →^coIy →^coIx+y/2


    ∀_x,y^nc.^coIx →^coIy → |x| ≤ y →1/4≤ y →^coIx/y

With those formulas the proof of this lemma is done by induction on n:

For n=0 it is easy since H(x,0)=1.

For arbitrary n we have H(x,n+1)=1/2( H(x,n)+x/H(x,n)).
By Lemma <ref> we have √(x)≤ H(n,x) and therefore √(1/16) =1/4≤ H(x,n) and x≤√(x)≤ H(x,n). Additionally the induction hypothesis  claims ^coI(H(x,n)). With (<ref>) we have ^coIx/H(x,n) and with (<ref>) we get ^coI(H(x,n+1))

With  and  we denote the computational content of (<ref>) and (<ref>). Each of these terms takes two streams of reals and returns a stream of their average and their quotient, respectively.
 As we have used induction over n to prove the lemma, the extracted term  is defined by recursion:

    v 0    := [1,1,…]
     v (n+1)   :=  ( v n) ( v ( v n))

This is actually Definition <ref> in the notation of streams.


    ∀^nc_x.^coIx → 0≤ x →^coI√(x)



We apply the introduction axiom of ^coI^+ to the predicate 
    Px := ∃_y(^coIy∧ 0≤ y∧√(y)=x):


    ∀_x^nc.   ∃_y(^coIy∧ 0≤ y∧√(y)=x)→
       ∀_x^nc(Px →∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2))
       →^coIx

In order to show the goal formula, we show the second premise:

    ∀_x^nc.   ∃_y(^coIy∧ 0≤ y∧√(y)=x) →
       ∃_d,x'^r(Sd d∧(^coI∪ P)x'∧|x|≤ 1∧ x=d+x'2)

Let x,y with ^coIy, 0≤ x and √(y) =x be given. |x|≤ 1 follows directly from 0≤ y≤ 1 and y=√(x). Form ^coIy we get d_1,d_2,d_3∈Sd and y'∈^co𝐈 with y=d_1d_2d_3y'. We do case distinction into three cases:

If y has one of the forms 1d_2d_3y', 01d_3y' or 001y', it follows y≤ 0 and therefore x=√(y)=0. Hence we define d:=0 and x':=0 then Sd d, ^coIx' and x=d+x'2 are obvious.

If y has one of the forms 000y', 001y', 011y' or 111y', we can rewrite y=00ey' for an e∈{0,1}. Here we define d:=0 and x':= √(e+y'/2). Then Sd d and 
    x=√(y)=√(e+y'/8)=√(e+y/2)/2=d+x'/2
 follow. Furthermore we have ^coI(e+y'/2) from Lemma <ref>, we have  0≤e+y'/2 since 0≤ y= e+y'/8 and we have √(e+y'/2) =x', therefore Px'.

The remaining case is that y has one the forms 010y', 011y', 111y', 110y', 10d_3y' or 11d_3y'. Here it follows 1/8≤ y. Therefore we have ∀_n ^coI(H(y,n)) by Lemma <ref>. Furthermore from Lemma <ref> we know that H(y) converges to √(y) with modulus ι:ℤ^+→ℕ. Thus by Theorem <ref> we have ^coIx. Applying x and ^coIx to ^coI^- leads to 
    ∃_d,x'^r(Sd d∧^coIx'∧|x|≤ 1∧ x=d+x'2),
 which proves the goal formular in this case.

With the definitions of  as the extracted term of Theorem <ref> and  as the extracted term of Lemma <ref> we have the following rules for the computational content :→ of this theorem:  

    (1u)    := [0,0,…]
    (01u)    := [0,0,…]
    (00u)    := 0  u
    (011u)    := 0  1u
    (111u)    := 0  1u
     u    :=  ι ( u)

The last rule shall only be applied if the other rules do not fit.
 
This algorithm has an inefficient runtime, which comes from the recursive definition of . In each step the function  is used twice and it has already a quadratic runtime. Therefore  has an at least exponential runtime and Haskell already takes quite a few minutes to compute even the first digit of [1,0,0,…]  10. By using  from Lemma <ref> instead of ι we get a bit more digits of  √(1/2).
If we enter 

in Haskell and wait about one minute, we get

as output. These are the first 29 digits of √(1/2)= √(2)/2. One can check that this result is indeed valid.

To improve the runtime it could be helpful to use another programming language and another data type but this is not the subject of this paper.


§ OUTLOOK AND SUGGESTIONS
 
It is possible to generalise the Heron sequence for roots of higher order. For a positive integer n and x∈ [0,1] we define

    G n x 0    := 1
    
    G n x (k+1)   := 1/n((n-1)(G n x k) + x/(G n x k)^n-1).

This sequence comes from Newton's method applied to the function y↦ y^n-x and one easily sees G 2 = H. With this formula and the modulus from Newton's method one could prove a general version of Theorem <ref> like

    ∀^nc_x.^coIx → 0≤ x →^coI√(x).

A difference of this generalisation is that one has to take a look at the first n+1 digits of the radicand to compute the first digit of the root. Another problem which increase the duration for higher roots is, that in the definition of G one divides by (G n x k)^n-1. If  (G n x k) is small, (G n x k)^n-1 is even smaller for large n and the smaller the divisor is, the longer is the duration of the 
division.

An generalisation of this is given in the context of functions with can be defined as power series. The main examples here are the trigonometric functions sin and cos. They are given by

    cos(x)    := ∑_k=0^∞(-1)^kx^2k/(2k)!
    sin(x)    := ∑_k=0^∞(-1)^kx^2k+1/(2k+1)!

From Chapter 8 Section 3 of <cit.> follows that for ^coIx i.e. x∈ [-1,1] one can choose ι:ℤ^+→ℕ as modulus of convergence for both sequences. To prove that every real in these sequences are in ^coI, one can use the formulas (<ref>) and (<ref>) and Lemma <ref>. 

As an even larger generalisation one could ask the question: How can we use Theorem <ref> to get an algorithm which takes a continuous function f from [-1,1]  to [-1,1] and returns a computable function which takes a signed digit stream of a real x and returns a signed digit stream of f(x)? Constructive continuous functions are for example defined in <cit.>.

Another direction in which one could extend this work is to replace the signed digit code by the Gray code. Similar to the signed digit code, the Gray code is also suitable to represent real numbers. In contrast to the signed digit code, the Gray code is unique and it has the property that a small change in the value of a real number effects only a small change in the Gray code of this real number. To implement Gray code, one needs simultaneously defined typs and simultaneously coinductivly defined predicates. In <cit.> there are the analogous statements of Lemma <ref> and Lemma <ref> and also analogous statements to the formulas (<ref>) and (<ref>) for Gray code proven. The main goal here would be to prove Theorem <ref> in terms of Gray code.



§ ACKNOWLEDGEMENTS

At the end of this paper the author would like to thank Peter Schuster and Helmut Schwichtenberg. Peter Schuster has supported him in writing this paper. Helmut Schwichtenberg has introduced him to the signed digit code and has given some background information. 

The author also thanks the Istituto Nazionale di Alta Matematica “F. Severi” (INdAM) for their financial support. The results presented in this paper were achieved during the authors  Marie Skłodowska-Curie fellowship by this institute.

