Iterative Square Root

[This article was first published on 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
## 
## ──────────────────────────────────────────────────────────────────────────────


To 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.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)