C++ program which calculates ln for a given variable x without using any ready functions

Question

I've searched for the equation which calculates the ln of a number x and found out that this equation is:

and I've written this code to implement it:

double ln = x-1 ;
    for(int i=2;i<=5;i++)
    {
        double tmp = 1 ;
        for(int j=1;j<=i;j++)
            tmp *= (x-1) ;
        if(i%2==0)
            ln -= (tmp/i) ;
        else
            ln += (tmp/i) ;
    }
    cout << "ln: " << setprecision(10) << ln << endl ;

but unfortunately I'm getting outputs completely different from output on my calculator especially for large numbers, can anyone tell me where is the problem ?

Answer 1

The equation you link to is an infinite series as implied by the ellipsis following the main part of the equation and as indicated more explicitly by the previous formulation on the same page:

In your case, you are only computing the first four terms. Later terms will add small refinements to the result to come closer to the actual value, but ultimately to compute all infinite steps will require infinite time.

However, what you can do is approximate your response to something like:

double ln(double x) {
  // validate 0 < x < 2
  double threshold = 1e-5;  // set this to whatever threshold you want
  double base = x-1;        // Base of the numerator; exponent will be explicit
  int den = 1;              // Denominator of the nth term
  int sign = 1;             // Used to swap the sign of each term
  double term = base;       // First term
  double prev = 0;          // Previous sum
  double result = term;     // Kick it off

  while (fabs(prev - result) > threshold) {
      den++;
      sign *=- 1;
      term *= base;
      prev = result;
      result += sign * term / den;
  }

  return result;
}

Caution: I haven't actually tested this so it may need some tweaking.

What this does is compute each term until the absolute difference between two consecutive terms is less than some threshold you establish.

Now this is not a particularly efficient way to do this. It's better to work with the functions the language you're using (in this case C++) provides to compute the natural log (which another poster has, I believe already shown to you). But there may be some value in trying this for yourself to see how it works.

Also, as barak manos notes below, this Taylor series only converges on the range (0, 2), so you will need to validate the value of x lies in that range before trying to run actual computation.

Answer 2

I believe the natural log in C++ language is simply log

Answer 3

That formula won't work for large inputs, because it would require you to take in consideration the highest degree member, which you can't because they are infinity many. It will only work for small inputs, where only the first terms of your series are relevant.

You can find ways to do that here: http://en.wikipedia.or/wiki/Pollard%27s_rho_algorithm_for_logarithms

and here: http://www.netlib.org/cephes/qlibdoc.html#qlog

Answer 4

It wouldn't hurt to use long and long double instead of int and double. This may get a little more accuracy on some larger values. Also, your series only extending 5 levels deep is also limiting your accuracy.

Using a series like this is basically an approximation of the logarithmic answer.

Answer 5

This version should be somewhat faster:

double const scale = 1.5390959186233239e-16;
double const offset = -709.05401552996614;

double fast_ln(double x)
{
    uint64_t xbits;
    memcpy(&xbits, &x, 8);
    // if memcpy not allowed, use 
    //       for( i = 0; i < 8; ++i ) i[(char*)xbits] = i[(char*)x];
    return xbits * scale + offset;
}

The trick is that this uses a 64-bit integer * 64-bit floating-point multiply, which involves a conversion of the integer to floating-point. Said floating-point representation is similar to scientific notation and requires a logarithm to find the appropriate exponent... but it is done purely in hardware and is very fast.

However it is doing a linear approximation within each octave, which is not very accurate. Using a lookup table for those bits would be far better.

Answer 6

This should work. You just needed the part where if x>=2 shrink x by half and add 0.6931. The reason for 0.6931 is that is ln(2). If you wanted to you could add if (x >= 1024) return myLN(x/1024) + 6.9315 where 6.9315 is ln(1024). This will add speed for big values of x. The for loop with 100 could be much less like 20. I believe to get exact result for an integer its 17.

double myLN(double x) {

    if (x >= 2) {
        return myLN(x/2.0) + 0.6931;
    }

    x = x-1;

    double total = 0.0;
    double xToTheIPower = x;

    for (unsigned i = 1; i < 100; i++) {

        if (i%2 == 1) {
            total += xToTheIPower / (i);
        } else {
            total -= xToTheIPower / (i);
        }

        xToTheIPower *= x;
    }

    return total;
}