Aproximating LogBASE in nim

I've split this tutorial into two parts, one that explains the math and one that explains how to implement the math. This tutorial requires knowledge of some math and basic nim.

Math

As many know, logarithms generate irrational numbers, so they can only be approximated, and the exact value is infinitely long. To approximate these numbers, we will be using Newton's method, which can approximate the inverse of many operations really well. Generally Newton's method is $x_{n+1}=x_n-{{f(x_n)}\over{f'(x_n)}}$ where $f(x)$ is calculated using $$x=\log_b(n)$$ $$n=b^x$$ $$0=b^x-n$$ $$f(x)=b^x-n$$ (This is how I calculate it, and might not be the correct method. This method supports the inverse of any operator, and I test it with nthRoot, logBASE, arcsin and arcocs)

In the end we get $x_{n+1}=x_n-{{b^x-n}\over{\ln(b)\times b^x}}$ with $\ln(b)\times b^x$ being the first derivative of $f(x)$. But, right now, we still cannot calculate the log function without using std/math's ln() function, so we must make our own. Thankfully it is not a infinitely recursive loop, and the first derivative of $f(x)=e^x-n$ is simply $f'(x)=e^x$, and I am allowing the use of std/math's $e$ during this tutorial. (Generating $e$ is not that difficult, but i will not be going through it in this tutorial.)

So, in the end, the functions you will be implementing are: $$x_{n+1}=x_n-{{e^{x_n}-n}\over{e^{x_n}}}$$ $$x_{n+1}=x_n-{{b^{x_n}-n}\over{\ln(b)\times b^{x_n}}}$$

Code (ln)

We start this by importing std/math, as we need floatPOW and almostEqual procs:

from std/math import pow, almostEqual, E # Almost equal is used becuause of float realted bugs

We will start with the ln function:

proc ln(n: float64): float64 =
  var next: float64 = 1.0

We start with the next as 1.0, and that as our "initial guess" for newton's method. Keep in mind that no matter the inital guess it usally adjusts itself to the right number, the only thing that the guess affects is speed. The next var will usally be $x_{n+1}$.

The next thing is the main loop:

while not almostEqual(result, next, 1): # (1): Very accurate, but can still catch bugs
  result = next
  next = result - ((pow(E, result) - n) / pow(E, result))

First, we set $n_x = n_{x+1}$ using result = next, then we do our math. $$r_esult-({{((e ^ {r_esult}) - n)}\over{(e ^ {r_esult})}})$$ That equation above is equal to result - ((pow(E, result) - n) / pow(E, result)), and we set that to out next var. We repeat that until the float no longer changes do to the decimal limit of a float64, and we get our result! Test it out with:

echo ln(5.0) # Check with a good calculator

Code (logBASE)

We first start out with same structure of ln():

proc log(n, base: float64): float64 =
  var next: float64 = 1.0
  while not almostEqual(result, next, 1): 
    result = next

Next we will split the equation into top and bottom, with the top being: $b^x-n$ and the bottom: $b^x\times \ln(b)$. We can represent these in code like this:

# (Inside the while loop)
let top = pow(base, result) - n
let bottom = pow(base, result) * ln(base)

Then we simpliy finish it off with:

next = result - (top / bottom) # Finish newton's method

You can try it by doing log(5, 10), and you should get $\approx 0.6989700043$!

Result

from std/math import pow, almostEqual, E

proc ln(n: float64): float64 =
  var next: float64 = 1.0
  while not almostEqual(result, next, 1): # (1): Very accurate, but can still catch bugs
    result = next
    next = result - ((pow(E, result) - n) / pow(E, result))

echo ln(5)

proc log(n, base: float64): float64 =
  var next: float64 = 1.0
  while not almostEqual(result, next, 1): # (1): Very accurate, but can still catch bugs
    result = next
    let top = pow(base, result) - n
    let bottom = pow(base, result) * ln(base)
    next = result - (top / bottom)

echo log(5, 10)