# Iterative Square Root

**rstats on Irregularly Scheduled Programming**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I saw a toot celebrating a short, clean implementation of a square root finding algorithm and wanted to dig a bit deeper into how it works, with a diversion into some APL.

This was the toot from Jim Gardner

Doubly pleased with myself.

Been doing the Tour of Go. Got to the section where you make a square root function, which should return once the calculated value stops changing. Struggled for ages. Trimmed and trimmed. Until finally… this!

The calculation for z was given, and I don’t understand it at all. But I don’t care. It was a total mess when I started and has turned out quite neat. I’m very satisfied.

But why “doubly pleased”? Because I’ve been solely using Neovim so far for Go!!

package main import ( "fmt" ) func Sqrt(x float64) float64 { z := 1.0 for { y := z z -= (z*z - x) / (2 * z) if z == y { return z } } } func main() { fmt.Println(Sqrt(16)) }

It’s a nice, not-too-complicated algorithm to play with, and I agree it’s hard to
see *why* it works for this application, so I thought it would be neat to walk
through that.

What we’re trying to solve here is the function \(y = x^2\) which we could write as \(f(x) = x^2 – y\) for which we want the value \(x\) where \(f(x) = 0\).

Newton’s Method is an iterative
method for solving equations of this type (not all equations, mind you – I have an
entire chapter of my PhD thesis discussing exactly *why* it can’t be used to
solve the equations I was solving that my supervisor insisted it could). It works
by using the slope (derivative) at a point to guide towards a solution. The formula
for the updated value \(x_{n+1}\) given some guess \(x_n\) is

\[ x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)} \] where \(f'(x)\) is the derivative of the function \(f\) at the point \(x\). For \(f(x) = x^2 – y\) the derivative is \(f'(x) = 2x\) so we can substitute this and \(f(x)\) into the above formula

\[ x_{n+1}=x_n−\frac{x_n^2-y}{2x_n} \] This is what the Go code calculates; given an initial guess of \(x_n = 1\) it calculates the next value as

x = x - (x*x - y) / (2 * x)

where, here, `y`

is the value we’re finding the square root of.

In R this could be written as

SQRT <- function(x) { z <- 1 while (TRUE) { y <- z z <- z - (z*z - x)/(2*z) if (abs(y -z) < 0.0001) return(z) } }

(since `base::sqrt`

is already defined) where I’ve used a tolerance rather than
relying on exact numerical equality. The `while(TRUE)`

construct is equivalent
to Go’s `for {}`

syntax; an infinite loop.

R actually has another way to write that which is even closer; `repeat {}`

SQRT <- function(x) { z <- 1 repeat { y <- z z <- z - (z*z - x)/(2*z) if (abs(y -z) < 0.0001) return(z) } }

One might notice that this approach requires essentially squaring a value, which is hardly expensive, but we can simplify and cancel out \(x_n\), so

\[ x_{n+1}=\frac{x_n-\frac{y}{x_n}}{2} \] in which case we have

SQRT <- function(x) { z <- 1 repeat { y <- z z <- (z + x/z)/2 if (abs(y -z) < 0.0001) return(z) } }

One of the reasons I wanted to dig into this was the fact that it’s a convergence…

In APL the power operator (`⍣`

aplwi)
applies a function some specified number of times, so

f ⍣n x

applies `f`

to `x`

`n`

times, i.e. `(f⍣3)x`

produces `f(f(f(x)))`

.

It can also be used as `⍣=`

where it will continue to apply the function until
the output no longer changes (is equal). A classic example is the
golden ratio; take the reciprocal
then add 1 until it converges, i.e.

\[ x_{n+1} = 1+\frac{1}{x_n} \]

which you can try for yourself here

1+∘÷⍣=1 1.618033989

In this, `+∘÷`

is the (tacit) function created by
composing (`∘`

) ‘addition of 1’ (`1+`

, a partial application of a function) and
‘reciprocal’ (`÷`

), which is iterated until it no longer changes (within
machine precision).

Iterating until convergence is exactly what we want, since we’re looking for the value satisfying

\[
x_n = x_{n+1}=\frac{x_n-\frac{y}{x_n}}{2}
\]
APL uses `⍵`

as the right argument placeholder and `⍺`

as the left, so the
function we want to apply repeatedly to the right argument is

{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}

If we provide `1`

as the right argument (the start value) and `16`

as the left
argument, we get

16{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}⍣=1 4

You can try this out yourself at tryapl.org (link should load that expression).

We can turn that into a function, once again using the argument placeholder

sqrt←{⍵{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}⍣=1} sqrt 25 5 sqrt 81 9

Taking the simplification above, we can write this a bit shorter as

sqrt←{⍵{(⍵+(⍺÷⍵))÷2}⍣=1} sqrt 144 12

As clean as the Go code looks, I think there’s a certain beauty to being able to write this in just 20 characters. It’s not for everyone, I get that.

I love these opportunities to learn a bit more about languages.

If you have comments, suggestions, or improvements, as always, feel free to use the comment section below, or hit me up on Mastodon.

##
`devtools::session_info()`

## ─ Session info ─────────────────────────────────────────────────────────────── ## setting value ## version R version 4.3.3 (2024-02-29) ## os Pop!_OS 22.04 LTS ## system x86_64, linux-gnu ## ui X11 ## language (EN) ## collate en_AU.UTF-8 ## ctype en_AU.UTF-8 ## tz Australia/Adelaide ## date 2024-05-29 ## pandoc 3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown) ## ## ─ Packages ─────────────────────────────────────────────────────────────────── ## package * version date (UTC) lib source ## blogdown 1.18 2023-06-19 [1] CRAN (R 4.3.2) ## bookdown 0.36 2023-10-16 [1] CRAN (R 4.3.2) ## bslib 0.6.1 2023-11-28 [3] CRAN (R 4.3.2) ## cachem 1.0.8 2023-05-01 [3] CRAN (R 4.3.0) ## callr 3.7.3 2022-11-02 [3] CRAN (R 4.2.2) ## cli 3.6.2 2023-12-11 [3] CRAN (R 4.3.2) ## crayon 1.5.2 2022-09-29 [3] CRAN (R 4.2.1) ## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.2) ## digest 0.6.34 2024-01-11 [3] CRAN (R 4.3.2) ## ellipsis 0.3.2 2021-04-29 [3] CRAN (R 4.1.1) ## evaluate 0.23 2023-11-01 [3] CRAN (R 4.3.2) ## fastmap 1.1.1 2023-02-24 [3] CRAN (R 4.2.2) ## fs 1.6.3 2023-07-20 [3] CRAN (R 4.3.1) ## glue 1.7.0 2024-01-09 [3] CRAN (R 4.3.2) ## htmltools 0.5.7 2023-11-03 [3] CRAN (R 4.3.2) ## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.3.2) ## httpuv 1.6.12 2023-10-23 [1] CRAN (R 4.3.2) ## icecream 0.2.1 2023-09-27 [1] CRAN (R 4.3.2) ## jquerylib 0.1.4 2021-04-26 [3] CRAN (R 4.1.2) ## jsonlite 1.8.8 2023-12-04 [3] CRAN (R 4.3.2) ## knitr 1.45 2023-10-30 [3] CRAN (R 4.3.2) ## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.2) ## lifecycle 1.0.4 2023-11-07 [3] CRAN (R 4.3.2) ## magrittr 2.0.3 2022-03-30 [3] CRAN (R 4.2.0) ## memoise 2.0.1 2021-11-26 [3] CRAN (R 4.2.0) ## mime 0.12 2021-09-28 [3] CRAN (R 4.2.0) ## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2) ## pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.3.2) ## pkgload 1.3.3 2023-09-22 [1] CRAN (R 4.3.2) ## prettyunits 1.2.0 2023-09-24 [3] CRAN (R 4.3.1) ## processx 3.8.3 2023-12-10 [3] CRAN (R 4.3.2) ## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.2) ## promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.2) ## ps 1.7.6 2024-01-18 [3] CRAN (R 4.3.2) ## purrr 1.0.2 2023-08-10 [3] CRAN (R 4.3.1) ## R6 2.5.1 2021-08-19 [3] CRAN (R 4.2.0) ## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.2) ## remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.3.2) ## rlang 1.1.3 2024-01-10 [3] CRAN (R 4.3.2) ## rmarkdown 2.25 2023-09-18 [3] CRAN (R 4.3.1) ## rstudioapi 0.15.0 2023-07-07 [3] CRAN (R 4.3.1) ## sass 0.4.8 2023-12-06 [3] CRAN (R 4.3.2) ## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2) ## shiny 1.7.5.1 2023-10-14 [1] CRAN (R 4.3.2) ## stringi 1.8.3 2023-12-11 [3] CRAN (R 4.3.2) ## stringr 1.5.1 2023-11-14 [3] CRAN (R 4.3.2) ## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.2) ## usethis 2.2.2 2023-07-06 [1] CRAN (R 4.3.2) ## vctrs 0.6.5 2023-12-01 [3] CRAN (R 4.3.2) ## xfun 0.41 2023-11-01 [3] CRAN (R 4.3.2) ## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.2) ## yaml 2.3.8 2023-12-11 [3] CRAN (R 4.3.2) ## ## [1] /home/jono/R/x86_64-pc-linux-gnu-library/4.3 ## [2] /usr/local/lib/R/site-library ## [3] /usr/lib/R/site-library ## [4] /usr/lib/R/library ## ## ──────────────────────────────────────────────────────────────────────────────

**leave a comment**for the author, please follow the link and comment on their blog:

**rstats on Irregularly Scheduled Programming**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.