# Iterative Square Root

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

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
##  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
##
## ──────────────────────────────────────────────────────────────────────────────