### A Pluto.jl notebook ### # v0.15.1 using Markdown using InteractiveUtils # ╔═╡ 9d5ecb30-32bb-4380-9986-01d160743f49 md" ## Babylonian's or Heron's Square Root Method (written-by PCM2021/07/26) *local file:* C:\Users\claus\PCM_Pluto&Julia\Pluto2021\MIT2021_ComputationalThinking\PCM20210717_MySquareRootMethod.jl *UOL file:* https://uol.de/f/2/dept/informatik/ag/lks/download/Probabilistic_Programming/JULIA/Pluto.jl/PCM20210717_MySquareRootMethod.jl " # ╔═╡ e741df04-912c-4488-a38f-23a15664593a md" ##### Context The iterative procedure for the square root function is a popular entry in programming courses. This is at least true for [Scheme](https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-10.html#%_sec_1.1.7) and [Julia](https://computationalthinking.mit.edu/Spring21/hw0/). " # ╔═╡ 11fd0742-3dca-4ddb-ac81-87575cb0b61e md" ##### Definition The specification of the iterative procedure derives from its [definition](https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-10.html#%_sec_1.1.7): $√{a} := \text{ the } x \text{ such that } x ≥ 0 \text{ and } x^2 = a.$ " # ╔═╡ 487f6053-4215-4229-81de-808ada762e2b md" ##### Pseudocode In their [Exercise 1](https://computationalthinking.mit.edu/Spring21/hw0/) Edelman et al.(2021) provide a solution sketch in pseudocode (here we interchanged x and a to meet the definition above). This exercise is titled *Square root by Newton's method* though we see that the presented pseudocode is a result of the much older [Babylonian or Heron's method](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method). It is this ancient method because, unlike Newton's, no derivatives are used here in the derivation and algorithm. We quote here the slightly modified (to meet the definition above) instruction by Edelman et al.(2021): **(Optional) [Exercise 1](https://computationalthinking.mit.edu/Spring21/hw0/) - Square root by Newton's (?) method** Given: $a$ Output: $√{a} = x$ 1. Take a guess $x$ 2. Divide $a$ by $x$ 3. Set $x :=$ the average of $a/x$ and $x$ 4. Repeat until $a/x$ is roughly equal to $x$. Return $x$ as the square root. In general, you will never get to the point where $a/x$ is exactly equal to $x$. So if our algorithm keeps going until $a/x == x$, then it will get stuck. So instead, the algorithm takes a parameter *error_margin*, which is used to decide when $a/x$ and $x$ are close enough to halt. $\blacksquare$ " # ╔═╡ 9f4044fe-7716-4834-ba97-88c4ff7d8d96 md" ##### Derivation of the Pseudocode How do we get the solution by iteration ? We can make a *case analysis*: Case 1. If we have a *perfect* $x = √{a}$ or or *near perfect* guess $x \approx √{a}$, we rewrite the [definition](https://en.wikipedia.org/wiki/Square_root) $a = x^2$ to into $\frac{a}{x} = x$ with *vanishing* error $ϵ = (\frac{a}{x} - x) = 0$ or with *acceptible* error $ϵ = (\frac{a}{x} - x) \approx 0$. An error $ϵ \approx 0$ is *acceptible* if its absolute value is below a positive *error margin* $ϵ_{max}$: $|ϵ| < ϵ_{max} \approx 0$. In both subcases, we are done. Case 2. If our guess a is *imperfect* then $x \ne \sqrt{a}$ and the error is $\epsilon = (\frac{a}{x} - x) \ne 0$. Then a 50%-reduction of the error is possible. This results in our *momentary step-size*: $\frac{1}{2}\epsilon=\frac{1}{2}(\frac{a}{x}-x).$ We get a *new improved* guess $x_{i+1}$, when we add the step size to the *previous* guess $x_i$: $x_{i+1} = x_i + \frac{1}{2}ϵ = x_i+\frac{1}{2}(\frac{a}{x_i}-x_i).$ After some simple manipulations we get the iteration formula: $x_{i+1} = x_i+\frac{1}{2}(\frac{a}{x_i}-x_i) = x_i+\frac{a}{2x_i}-\frac{x_i}{2} = \frac{x_i}{2} + \frac{a}{2x_i} = \frac{x_i^2+a}{2x_i} = \frac{1}{2}(\frac{a}{x_i} + x_i)$ , and in the limit $√{a} = \lim_{i \rightarrow \infty} x_i\;\;\;\blacksquare$ " # ╔═╡ e71638ba-0833-4c73-9eb8-a7e57c6fd2c1 md" ##### Good Choice for Initial Guess $x$ The *improved* error $ϵ_{i+1} = |(\frac{1}{2}(\frac{a}{x}+x)-\sqrt{a})|$ should be smaller than the *initial* error $ϵ_i = |(x-√{a})|$. As shown below this is true if the *initial* guess $x$ is greater than the *true* but *unknown* square root (absolute values assumed). So in our Julia-code we set the *initial* value x to the known a. $|ϵ_{i+1}| < |ϵ_i|$ $\Rightarrow |x_{i+1} - √{a}| < |x_i - √{a}|$ $\Rightarrow |\frac{1}{2}(\frac{a}{x_i}+x_i)-√{a}|<|x_i-√{a}|$ $\Rightarrow |\frac{1}{2}(\frac{a}{x_i}+x_i)|<|x_i|$ $\Rightarrow |\frac{1}{2}(\frac{a+x_i^2}{x_i})|<|x_i|$ $\Rightarrow |a+x_i^2| < |2x_i^2|$ $\Rightarrow |a| < |2x_i^2 - x_i^2|$ $\Rightarrow |a| < |x_i^2|$ $\Rightarrow |√{a}| < |x_i|\;\; \blacksquare$ " # ╔═╡ 7df9abad-3795-4f65-ba09-5a12aa403c96 # ╔═╡ 43d08b6d-d92c-4f37-899c-1c8c171de4d7 md" ##### Julia/Pluto.jl-Code " # ╔═╡ e3f99d39-7f42-46c2-a7fd-5168d121bcf3 #------------------------------------------------------------------------------------- # Julia-code by *** PCM 2021/07/26 *** # # a is obligatory argument of square root funtion # errorMax is default argument and is the maximal not tolerated error # x is a default keyword argument for the initial guess x0 # verbose is a default keyword argument; # if true then intermediate results are collected #------------------------------------------------------------------------------------- mySquareRoot(a::Float64, errorMax = 1.0E-6; x0::Float64 = a, verbose::Bool=false) = begin #------------------------------------------------------------------ # local definitions error(x) = a/x - x # if x = √a then error(a, x) = 0 newGuess(x, stepSize) = x + stepSize #------------------------------------------------------------------ x = x0 it = 0 intermediateSteps = [] tempError = error(x) verbose ? push!(intermediateSteps, ("it=$it", "x=$x", "err=$tempError")) : " " while abs(tempError) > errorMax it = it + 1 stepSize = tempError / 2 x = newGuess(x, stepSize) # new improved guess tempError = error(x) # if a = √x then error(x, a) = 0 verbose ? push!(intermediateSteps, ("it=$it", "x=$x", "err=$tempError")) : " " end return verbose ? intermediateSteps : x end #------------------------------------------------------------------------------------- # ╔═╡ 01071894-b251-483c-a255-5a38d0fcad79 mySquareRoot(2.0) # ╔═╡ 60071ee7-acc5-4278-81bf-1400fab1d0fc mySquareRoot(2.0, x0=1.E-1) # ╔═╡ 1379cf59-a7b9-41af-b09f-2f53b9e56da7 mySquareRoot(2.0)^2 # ╔═╡ e883b914-b69d-435b-b235-f8163d01e298 mySquareRoot(125348.0) # ╔═╡ 8fd66d80-3b04-4888-a5aa-4d80cb51443b mySquareRoot(125348.0)^2 # ╔═╡ b0a6aec9-127b-4632-840a-fe82025d0af5 mySquareRoot(125348.0; verbose=true) # ╔═╡ 0067950f-8438-4f7a-a732-46ad9ef64c3f md" Compare the intermediate results of the call with x=600 (below) with those [here](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method) " # ╔═╡ 64952529-f2c5-4b70-8690-a071367a31e1 mySquareRoot(125348.0, x0 = 600.0; verbose=true) # ╔═╡ Cell order: # ╟─9d5ecb30-32bb-4380-9986-01d160743f49 # ╟─e741df04-912c-4488-a38f-23a15664593a # ╟─11fd0742-3dca-4ddb-ac81-87575cb0b61e # ╟─487f6053-4215-4229-81de-808ada762e2b # ╟─9f4044fe-7716-4834-ba97-88c4ff7d8d96 # ╟─e71638ba-0833-4c73-9eb8-a7e57c6fd2c1 # ╠═7df9abad-3795-4f65-ba09-5a12aa403c96 # ╟─43d08b6d-d92c-4f37-899c-1c8c171de4d7 # ╠═e3f99d39-7f42-46c2-a7fd-5168d121bcf3 # ╠═01071894-b251-483c-a255-5a38d0fcad79 # ╠═60071ee7-acc5-4278-81bf-1400fab1d0fc # ╠═1379cf59-a7b9-41af-b09f-2f53b9e56da7 # ╠═e883b914-b69d-435b-b235-f8163d01e298 # ╠═8fd66d80-3b04-4888-a5aa-4d80cb51443b # ╠═b0a6aec9-127b-4632-840a-fe82025d0af5 # ╟─0067950f-8438-4f7a-a732-46ad9ef64c3f # ╠═64952529-f2c5-4b70-8690-a071367a31e1