Small algorithms · Viv docs

Newton's square root

Repeat x = (x + a/x) / 2; each step roughly doubles the number of correct digits.

Viv [1]

def newton_sqrt(a) {
    x = a
    for i to 20 {
        x = 0.5 * (x + a / x)
    }
    return x
}
show newton_sqrt(2.0)
show sqrt(2.0)

Bisection

Root of cos(x) on [1, 2]. Each step halves the interval, so 40 halvings leave it 2^-40 ≈ 9e-13 wide and give pi/2 to about twelve digits; show prints the first six.

Viv [2]

def bisect_cos(lo, hi) {
    for i to 40 {
        mid = 0.5 * (lo + hi)
        if cos(lo) * cos(mid) < 0.0 {
            hi = mid
        } else {
            lo = mid
        }
    }
    return 0.5 * (lo + hi)
}
show bisect_cos(1.0, 2.0)

gcd via mod

Euclid's algorithm. The floats hold whole numbers exactly, so mod is exact.

Viv [3]

def gcd(a, b) {
    if b < 0.5 {
        return a
    }
    return gcd(b, mod(a, b))
}
show gcd(1071.0, 462.0)

Horner evaluation

A polynomial as an array of coefficients, highest power first.

Viv [4]

def horner(coeffs, x) {
    acc = 0.0
    for i to len(coeffs) {
        acc = acc * x + coeffs[i]
    }
    return acc
}
show horner([2.0, 0.0 - 3.0, 5.0], 2.0)

The power method

The dominant eigenvalue of a symmetric 2x2, found by repeated matrix-vector products; compare eig_sym.

Viv [5]

def normalize(v) {
    n = sqrt(v[0] * v[0] + v[1] * v[1])
    return v / n
}
def rayleigh(A, v) {
    w = A @ v
    return v[0] * w[0] + v[1] * w[1]
}
A = [[2.0, 1.0], [1.0, 3.0]]
v = [1.0, 0.0]
for i to 30 {
    v = normalize(A @ v)
}
show rayleigh(A, v)

Compensated summation

Kahan's method carries each addition's rounding error into the next one. Adding 0.1 a million times this way gives exactly 100000.

Viv [6]

def kahan_sum(n, val) {
    s = 0.0
    comp = 0.0
    for i to n {
        y = val - comp
        t = s + y
        comp = (t - s) - y
        s = t
    }
    return s
}
show kahan_sum(1000000.0, 0.1) - 100000.0

Fibonacci two ways

Naive recursion, whose number of calls grows exponentially with n, against a loop that takes n steps.

Viv [7]

def fib_rec(n) {
    if n < 1.5 {
        return n
    }
    return fib_rec(n - 1.0) + fib_rec(n - 2.0)
}
def fib_loop(n) {
    a = 0.0
    b = 1.0
    for i to n {
        t = a + b
        a = b
        b = t
    }
    return a
}
show fib_rec(15.0)
show fib_loop(15.0)