3

I want to implement the formula:

enter image description here

#include <iostream>
#include <cmath>
#include <limits>

using namespace std;

int main(){
    double a, eps = numeric_limits<double>::epsilon();
    cout << "a=";
    cin >> a;
    double w = sqrt(a), prod = 1 + w;
    int n = 1;
    do {
        w = sqrt(w);
        prod *= 1 + w;
        cout << "ln(" << a << ")=" << ldexp(a-1, n)/prod << endl;
        n++;
    } while(abs(w - 1) > eps);

    return 0;
}

But e.g. ln(2.78787)=0.512639 and that cannot be true. Where is the mistake?

ATW
  • 233
  • 2
  • 14
  • These links may be relevant: [Is floating point math broken?](https://stackoverflow.com/q/588004/5910058) and [What Every Computer Scientist Should Know About Floating-Point Arithmetic](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) – Jesper Juhl Aug 02 '20 at 16:21
  • You should start with `prod = 1`, not `prod = 1 + w`. You effectively compute the product to `n+1`, not to `n`. – Igor Tandetnik Aug 02 '20 at 16:23
  • Yes @JesperJuhl but that isn‘t just a tiny difference which I would ascribe to floating point arithmetic. The difference is ~0.5. – ATW Aug 02 '20 at 16:23
  • @ATW I don't have an answer to your question. I just wanted to point you at some, potentially relevant, resources. – Jesper Juhl Aug 02 '20 at 16:25

1 Answers1

4

Initial conditions for the loop are off. Start with

double w = a, prod = 1;

and the program would produce expected results

Igor Tandetnik
  • 50,461
  • 4
  • 56
  • 85