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 xapplies 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.618033989In 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
4You 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
9Taking the simplification above, we can write this a bit shorter as
sqrt←{⍵{(⍵+(⍺÷⍵))÷2}⍣=1}
sqrt 144
12As 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
##
## ──────────────────────────────────────────────────────────────────────────────