Computing Square Roots with the Newton-Raphson Method

See:

Newton–Raphson method
"Why Functional Programming Matters" by John Hughes

We can approximate the square root a of a number n by successive applications of this formula:

aᵢ₊₁ = (aᵢ + n/aᵢ)/2

We repeat until aᵢ is within some error tolerance of the actual root.

We use Rational Numbers to approximate reals. Given the number n and our initial guess for the root a as numerator denominator pairs together on the stack we can calculate the next aᵢ₊₁ like so:

                 nₙ dₙ nₐ dₐ
-----------------------------------------------
   nₙ dₙ nₐ dₐ div-frac nₐ dₐ add-frac 1 <<

(Note: to divide a ratio by two it suffices to double the denominator, this is achieved simply and efficiently by shifting it left by one bit.)

Let's start with:

dup2 [div-frac] dipd

Here's how it works:

          nₙ dₙ nₐ dₐ dup2 [div-frac] dipd
------------------------------------------------------
         nₙ dₙ nₐ dₐ nₐ dₐ [div-frac] dipd
------------------------------------------------------
             nₙ dₙ nₐ dₐ div-frac nₐ dₐ
------------------------------------------------------
              nₙ/ₐ dₙ/ₐ nₐ dₐ add-frac
------------------------------------------------------
                   nₙ/ₐ₊ₐ dₙ/ₐ₊ₐ

So this would be the "get next aᵢ₊₁" function:

dup2 [div-frac] dipd add-frac 1 <<

Note that this consumes the nₙ dₙ pair, but if we use it with the binary combinator we would only get the new denominator, the new numerator would be discarded! What to do?

I guess we just need a dup4?

[4 grabN] nullary i

This is elaborate, but it would compile to something simple, e.g. in Python:

@SimpleFunctionWrapper
def dup4(stack):
    a0, a1, a2, a3, _ = get_n_items(4, stack)
    return a0, (a1, (a2, (a3, stack)))

I like to call these sorts of functions "stack chatter". They serve to prove that we know what we're doing, then vanish in the compiler.

≡ `dup4`

[dup4 [4 grabN] nullary i] inscribe

Once we have that our function would be:

dup4 div-frac add-frac 1 <<

Like so:

       nₙ dₙ nₐ dₐ dup4 div-frac add-frac 1 <<
------------------------------------------------------
   nₙ dₙ nₐ dₐ nₙ dₙ nₐ dₐ div-frac add-frac 1 <<
------------------------------------------------------
        nₙ dₙ nₐ dₐ nₙ/ₐ dₙ/ₐ add-frac 1 <<
------------------------------------------------------
            nₙ dₙ nₙ/ₐ₊ₐ dₙ/ₐ₊ₐ 1 <<
------------------------------------------------------
             nₙ dₙ nₙ/ₐ₊ₐ 2dₙ/ₐ₊ₐ

That gives us our new ratio for aᵢ₊₁ and retains n on the stack below it.

But we will want to compare this aᵢ₊₁ to aᵢ to check if the (absolute value of) the difference between them is less than some error tolerance. (I was going to square the number and compare to n but I reread the paper and it's done the other way.) Keeping six numbers on the stack is a little much. What if we used a Generator Program?

nₙ dₙ nₐᵢ dₐᵢ [nₐᵢ₊₁ dₐᵢ₊₁ F] x
nₙ dₙ nₐᵢ dₐᵢ [nₐᵢ₊₁ dₐᵢ₊₁ F] nₐᵢ₊₁ dₐᵢ₊₁ F
nₙ dₙ nₐᵢ₊₁ dₐᵢ₊₁ [nₐᵢ₊₂ dₐᵢ₊₂ F]

Computing Square Roots with the Newton-Raphson Method

≡ dup4

≡ `dup4`