ansaurus

Question

Answer 1

+4 A:

You should not divide by i and then assign to ty, devide by i after the addition, i.e. logx = logx + ty / i.

Edit this should work:

i=1;   // (i.e. not 3)
logx = 0 ;
ty = (x-1)/(x+1) ;
do
{
    logx = logx + ty / i;
    tty = ty ;
    ty = (ty * ((x-1)/(x+1)) * ((x-1)/(x+1)));
    i = i + 2 ;
} while(tty - ty > 0.0000005 );

mvds 2010-07-27 12:08:27

Yes it is. mathoverflow is for reseach. This is very much an implementation question.

duffymo 2010-07-27 12:09:40

now there's a `c` tag on the question, I agree ;-)

mvds 2010-07-27 12:14:29

Yup. That looks better. And although its outside the scope of the question I personally would prefer to see the logx assignment happen at the end of the loop. Otherwise you waste your last loop (you get the next ty but then never use it).

Chris 2010-07-27 12:19:03

That, plus storing `(x-1)/(x+1) * (x-1)/(x+1)` once outside the loop. `gcc -O3` apparently doesn't notice that the factor is the same every iteration.

mvds 2010-07-27 12:25:51

thank you. I get it.

Barshan Das 2010-07-27 18:19:33

Answer 2

A:

you're approximating a function with a finite number of iterations (the stop condition is tty-ty > ...); decrease the 0.00000005 and you should obtain better precision

then you might run into numerical precision problems, so try using double instead of float if you're not already doing so

finally remember that generally logarithms are irrational numbers so they just can't be represented with a finite number of digits

Nicola Montecchio 2010-07-27 12:09:09

Answer 3

A:

ISTR that the rate of convergence of log series are notoriously slow.

Visage 2010-07-27 12:09:11

Answer 4

+3 A:

Actually what OP's computing is

$2\sum_{n=1,n\text{ odd}}^\infty\frac1{n!!}\left(\frac{m-1}{m+1}\right)^n = \sqrt{2\pi} \exp\left(\frac12\left(\frac{m-1}{m+1}\right)^2\right)\operatorname{erf}\left(\frac1{\sqrt2}\frac{m-1}{m+1}\right)$

instead of

$2\sum_{n=1,n\text{ odd}}^\infty\frac1{n}\left(\frac{m-1}{m+1}\right)^n = 2\tanh^{-1}\frac{m-1}{m+1} = \ln m$

BTW, it's better to factor out the ((x-1)/(x+1)) * ((x-1)/(x+1)) to avoid recomputing it. (The compiler may or may not treat this as a loop invariant and take it out.)

double ty = (x-1)/(x+1);
double p = ty * ty;
int i = 1;
double logx = 0;
do {
  logx += ty / i;
  tty = ty;
  ty *= p;
  i += 2;
} while (tty/(i-2) - ty/i > 5e-6);

KennyTM 2010-07-27 12:20:23

thank you very much.

Barshan Das 2010-07-27 15:31:07

Answer 5

+1 A:

Your C implementation is not correct for the algorithm.

log(m) base e = 2[ (m-1/m+1) + (1/3) * (m-1/m+1)^3 + (1/5) * (m-1/m+1)^5 + ... ]
                                ^^^ what is this      ^^^ and this

Your algorithm:

i = 3;
logx = 0 ;
ty = (x-1)/(x+1) ;
do
{
    logx = logx + ty ;
    tty = ty ;
    ty = (ty * ((x-1)/(x+1)) * ((x-1)/(x+1))) / i ; 
//  ^^ when i = 3 ty after this line = ((x-1)/(x+1) * ((x-1)/(x+1))^2) / 3
//     when i = 5 ty after this line = ((((x-1)/(x+1)^3) / 3) * ((x-1)/(x+1))^2) / 5 
//                                   = ((x-1)/(x+1)^5) / 15 <<< wrong divisor (see above, should be 5)
//     and it gets worse for each iteration of the loop
    i = i + 2 ;
} while(tty - ty > 0.0000005 );
logx = 2*logx ;
printf("\n ln (%g) = %g \n", x, logx);

JeremyP 2010-07-27 12:28:40

thank you JeremyP. You perfectly showed me where I am doing wrong.Thanks a lot.

Barshan Das 2010-07-27 15:25:08

ansaurus

tags:

views:

answers:

value of natural logarithm

related questions