Calculate ln(x) in C

[Bad.boy!]

I want to calculate ln(x) in C and don't want to use the standard libraries. So I wrote my own function, but I have problems with big numbers (around 2.0 and bigger). Looping for longer doesn't make the answer correct. Because google gives the same answer as my calculator and does it faster. And when I use it in another algorithm the values suddenly get bigger until I reach an overflow (inf). Can anyone tell me what's wrong?

http://upload.wikimedia.org/math/1/b/e/1be...a968fa523f5.png

 

 

// macht = powerint macht_i(int a, unsigned int b){int i, resultaat = a;if (b > 1) { for (i = 1; i < b; i++) { 	resultaat *= a;  } return resultaat;}else { switch (b) { case 0: 	return 1; case 1: 	return a; }}}double macht_d(double a, int b){int i;double resultaat = a;if (b < 1) { b *= -1; for (i = 1; i < b; i++) { 	resultaat *= a;  } return (1/resultaat);}else if (b > 1) { for (i = 1; i < b; i++) { 	resultaat *= a;  } return resultaat;}else { switch (b) { case -1: 	return (1/a); case 0: 	return 1; case 1: 	return a; }}}double ln(double x) {int i, bufferI;double bufferDI, bufferD, resultaat = 0.0;// Eigenlijk moet dit oneindig vaak herhaald worden, maar bij 1000 keer heb je het wel redelijk precies.for (i = 1; i < 5000; i++) { // (-1)^n+1 bufferI = -1; bufferI = macht_i(bufferI, (i+1)); bufferDI = bufferI;  // ((-1)^n+1)/n bufferDI /= i;  // (x-1)^n bufferD = (x-1); bufferD = macht_d(bufferD, i); // Vermenigvuldigen en optellen (of aftrekken) bufferD *= bufferDI; resultaat += bufferD;}return resultaat;}

 

 

EDIT: Solved, the function I used was only working for -1<x<1. I fixed it using ln(x) = -ln(1/x).

Edited January 17, 2013 by Bad.boy!

[K^2]

Yup. Here is the math for it, if you are interested.

 

The function is derived via Taylor Series of ln(x) around x=1.

 

f(x) = ln(x), f'(x) = 1/x, f''(x) = -1/x². f'''(x) = 2/x³ ...

f(1) = 0, f'(1) = 1, f''(1) = -1, f'''(1) = 2

 

g(x) = 0 + 1(x-1) - 1(x-1)²/2 + 2(x-1)³/6 - ... = Sum (-1)^(n+1) (1/n) (x-1)^n

 

The x-1 comes from the fact that we expanded around x=1. Now, the theorem tells you that if the series converges for a specific value of x, then g(x) = f(x). However, the series is not guaranteed to converge for all x, or indeed, for any x.

 

This is where Radius of Convergence comes in. It basically is the question of "Which values of x does the series converge for?"

 

The series for logarithm alternates. (Changes sign.) That means that series converges only if the sequence converges to zero. In other words, (x-1)^n / n must go to zero for large n. So lets look at that. I'm going to substitute z = x-1.

 

lim n->inf z^n / n = lim n->inf n z^(n-1) = Infinity for |z| > 1 and 0 for |z<1|. (Using L'Hopital's Rule.)

 

So the sequence converges for -1 < z < 1, which translates to 0 < x < 2.

 

And that's precisely what you saw. Your function worked for x < 2, but did not for x > 2, simply because the sequence would not converge.

 

 

Do keep in mind that using inverse method is correct, but may come at a penalty for precision.

 

 

 

An alternative algorithm for logarithm computation is based on Newton's method. Observe that if x = ln(a), then a = e^x. In other words, the zero of exp(x) - a gives you ln(a). Because Newton Method converges well for this particular function, thanks to exp(x) function being monotonous, you can get arbitrary precision for any value of a. However, that relies on you computing exp(x) many times, and that has it's own collection of methods. Partly thanks to the fact that the series for exp(x) converges for all x.

 

Numerical computing often comes to a compromise between speed and precision. Depending on what you are trying to achieve, one may be more important than the other.

 

 

It seems that you are doing this as just an exercise, but if you have a specific purpose in mind, maybe I can help you find a better solution.

[Bad.boy!]

I'm just doing this as an exercise, but thanks for the background info.