Approximating Natural Log By Hand

March 15, 2012 7:54 AM

Calculating/approximating ln(x) 'by hand'?

I've got an application that only has a limited set of mathematical operators: +, -, parens (), multiplication/division, pow() or **, mod, and div.

I've got a request to calculate ln(x). How can I best approximate this using the operators available to me, if at all? (As far as I can tell, x will be between 0.00 and 100.00 - two digit precision.)

posted by unixrat to science & nature (18 answers total)

You can use a Taylor approximation.
posted by Tooty McTootsalot at 8:01 AM on March 15, 2012 [4 favorites]

Given the limitations of the application, could you create a log table in some way?
posted by ManInSuit at 8:05 AM on March 15, 2012 [1 favorite]

Tooty's answer is right.
posted by escabeche at 8:07 AM on March 15, 2012

Wikipedia lists several ways to calculate a logarithm, including some more efficient ways than the Taylor approximation.
posted by grouse at 8:15 AM on March 15, 2012 [1 favorite]

(Bear with me here, it's been a long time)

So looking at that Taylor approximation, would this be an accurate implementation?

If I start with x...

I can do
x**-1 = y
and then
y-((y**2)/2)+((y**3)/3)-((y**4)/4)... and so on?
posted by unixrat at 8:22 AM on March 15, 2012

No, the power series you give is computing ln(1+y) so you want to let y be (x**-1 - 1) and then apply it.
posted by escabeche at 8:36 AM on March 15, 2012

Great! So to recap:

Starting with x...

((x**-1)-1)=y

and then...

y-((y**2)/2)+((y**3)/3)-((y**4)/4)+((y**5)/5)...

...would give me an approximation of ln(x), correct?
posted by unixrat at 8:48 AM on March 15, 2012

The taylor series is very slowly convergent (the wiki article grouse links shows you the convergence after 100 iterations. You also need to include a check for whether or not x is greater than 1: you don't want to invert x unless x>1.

The wiki article's other suggestion, using the arctan formulation, is probably better (and is labelled as more efficient).

y=(x-1)/(x+1);
z=x+x**3/3+x**5/5+x**7/7+....;
z=2*z;

gives you z ~ log(x). If you want 2 digit precision in the range [0.01 to 100], you need 80 terms it looks like (i.e. the last term in the sum should be x**161/161). Using the taylor form, you would need 319 terms in order to get two digits for x=0.01 or x=100.

Can you perform a loop? (do, while, for, etc)? If so, you can compute z pretty quickly, without having to write this out by hand.
posted by bessel functions seem unnecessarily complicated at 9:01 AM on March 15, 2012

Yikes. I can't loop at all.

I'm not sure I follow the arctan as laid out in bessel's example immediately above. If x is my input, why do I calculate a 'y'? It doesn't seem to be reused.

"Can't be done with these tools" is also an acceptable answer even if the client will be unhappy. Better unhappy than inaccurate.
posted by unixrat at 9:10 AM on March 15, 2012

oops, sorry, that should be
y=(x-1)/(x+1);
z=y+y**3/3+y**5/5+y**7/7+....;
z=2*z;
Apologies for the confusing typo.

If you're ok hand-writing out all 80 terms, it will be accurate to two decimal digits over [0.01,100] it seems.
posted by bessel functions seem unnecessarily complicated at 9:14 AM on March 15, 2012

Okay, I understand that. One last clarification - it should be all addition, right? The Taylor alternated between addition and subtraction.

Many thanks to everyone for their help.
posted by unixrat at 9:51 AM on March 15, 2012

If it were me doing this, I would probably graph ln(x) as determined by my computer's math libraries minus ln(x) as determined by this method for many numbers in the expected domain to make sure I was doing it right, and also to show the client what the expected accuracy is. Of course, your application might introduce additional errors or imprecisions.
posted by grouse at 10:02 AM on March 15, 2012

Yep, all addition. One final point to mention: You're computing y**3, then y**5, then y**7, and so on up to y**161. This is actually a really slow thing to do. If you're actually going to use this, a faster implementation would be

y=(x-1)/(x+1);
w=y; // w is (x-1)/(x+1)
y=y*y; // now, y=(x-1)**2/(x+1)**2
z=0;

z=z+w;
w=w*y; // now w=(x-1)**3/(x+1)**3
z=z+w/3;
w=w*y; //now w=(x-1)**5/(x+1)**5
z=z+w/5;
w=w*y; //now w=(x-1)**7/(x+1)**7
z=z+w/7;
....
// 77 more terms following the same pattern
// constantly updating w.
....
z=z*2;


That is: use the previously computed w=(x-1)**5/(x+1)**5 to compute the new w=(x-1)**7/(x+1)**7. I'm not sure if the speed of evaluation is much of an issue, but this will be an improvement.
posted by bessel functions seem unnecessarily complicated at 10:21 AM on March 15, 2012

Can you make recursive functions or use a stack data structure to simulate recursion? If so, you might try using the Binary Logarithm Algorithm and then the logarithm change-of-base identity. ln(x) = lg2(x)/lg2(e). Store the constant 1/lg2(e) in your code.
posted by scose at 10:22 AM on March 15, 2012

Okay, so I was able to max out the application's mathematics capabilities after 15 iterations of Bessel's suggestions. The numbers feel a bit rough but I've got a QA/QC team reviewing the output against the suggested graph of ln() as per grouse's suggestion.

A million thanks to all for helping, no pun intended. :)
posted by unixrat at 10:33 AM on March 15, 2012

I had a look at the implementation of log and log2 in eglibc, a free software C library. Unfortunately, these implementations (they were several) relied on being able to find the base2 exponent of the floating-point number (like frexp in C), and then compute the log₂ or natural log of a number in the restricted domain [.5, 1), which leads to much smaller polynomial. This doesn't seem to be applicable to the restrictions you are working under.
posted by jepler at 1:37 PM on March 15, 2012

You have variables. Do you have arrays? If so, and the progam storage can handle it, I would think about just externally pre-computing ln(x) for 0.00 .. 100.00, and make it an array lookup ArrayName[x*100]. That's only 10001 elements in the array.
posted by fings at 9:22 PM on March 15, 2012

One very important trick is to reduce the range. All polynomial approximations have a range where they work well, and other ranges where they work poorly. Use the fact that ln(x) = log2(x)/log2(e) to work with base 2. Then use the fact that log2(x/2) = log2(x) - 1, and divide x by 2, N times, to change the problem to calculating log2(y), where y is within the correct range. When you're done, add N to the answer, divide by log2(e), and you're done.
posted by ubiquity at 9:33 AM on March 19, 2012