Calculating the function and checking its domain

Question

my task is to calculate the following function:

∛log₅(x² - √x)

for argument x = [-5 ; 10], with increment 0.2. So in order for x to satisfy the function's domain, I guess we have to test it against 2 cases:

double myPow(double x, double y)
{
    if (x < 0.0)
        return -pow(-x, y);
    else
        return pow(x, y);
}
//...
for(double x = -5.0 ; x <= 10.0 ; x += 0.2)
{
    //...
    if(x >= 0.0 && (pow(x, 2.0) - sqrt(x)) > 0.0)
    {
        //myPow correctly calculates the cube root if x < 0 (see def. of std. pow)
        double sum = myPow(log(pow(x, 2.0) - sqrt(x))/log(5), 1/3.0);
        printf("%f ", sum);
    }
    //...
}

However, when x == 1.0, the second condition should equal to:

pow(1.0, 2.0) - sqrt(1.0) == 0.0

And it does, when I write it by itself as above, but in my for loop it doesn't - it equals to something like that: 0.00000000000000310862...

Why is this happening and how can I rectify it to not pass the if condition (and thus satisfy the logarithm constraints)?

Can you please replace the math notation (eg; `[-5 ; 10]`) with more general/common terminology? For example, I'm not sure if this is a open or closed interval. Also I think you're just running into rounding errors in the way floating point works — Avery, Jun 23 '14 at 00:38
Accumulated roundoff. Loop from -50 to 100 by 2 and divide by 10 inside the loop. (This is a duplicate, but of what I can't find.) — tmyklebu, Jun 23 '14 at 00:39
Where is the cube root in your formula? The formula as written is only defined for x>1, no need to check complicated conditions. — n. m. could be an AI, Jun 23 '14 at 05:17
@tmyklebu from the math's standpoint, your solution is equivalent to mine. Is there any article that explains this floating-point arithmetic and round errors **intuitively**? - I find it confusing... — user3765753, Jun 23 '14 at 09:20
@n.m. in `myPow` function call - I raise the logarithm to the power of 1/3 which is a cube root. — user3765753, Jun 23 '14 at 09:21
@user3765753: I don't know if you can read Pascal. If so, I wrote an article explaining these things from my perspective (non-mathematician, I'm actually a dentist), but using Delphi as example: http://rvelthuis.de/articles/articles-floats.html. I don't know if that helps. — Rudy Velthuis, Jun 23 '14 at 10:08
@Avery: `[a ; b]` means that the "edges" are inclusive, so `a <= x <= b`. — Rudy Velthuis, Jun 23 '14 at 10:37
@user3765753: `printf("%.30f\n", 0.2);` ought to clear things up a little. — tmyklebu, Jun 23 '14 at 13:51
@tmyklebu I still don't get it, those outputs: `for(double x = -50, x2 = -5 ; x <= 100 ; x += 2, x2 += 0.2) { printf("%.30f\t", x/10.0); printf("%.30f\n", x2); }` are very similar (and still not accurate). — user3765753, Jun 23 '14 at 14:59
@user3765753: The first one gives you the closest `double` to `x / 10` because `x` is computed exactly. The second one gives you the sum of a bunch of closest-`double`s-to-0.2. — tmyklebu, Jun 23 '14 at 15:40
@tmyklebu I get it know. 1) So I should always check whether the increment can be represented exactly with `printf("%.30f\n", /*literal_increment_here*/);` and if not - look for such workarounds like yours? 2) Your soulution makes x not pass the second condition just like I want it to, so is it preferable to additionally include comparisons involving epsilon, e.g. `const double EPS = 0.000001;` for "just in case"? — user3765753, Jun 23 '14 at 15:52
Use `x * x` instead of `pow`. [`pow` is usually crap.](http://stackoverflow.com/questions/16881749/different-results-on-using-sqrt-and-pow-functions/16881975#16881975) — tmyklebu, Jun 23 '14 at 15:54
@tmyklebu Okey. And could you answer my two questions from the previous comment? :) — user3765753, Jun 23 '14 at 16:22
Why would you toss an epsilon into the mix? Stop using `pow` and `pow`'s crappiness will go away. — tmyklebu, Jun 23 '14 at 17:12
@tmyklebu OK, I'll stop using it. I recalled that the epsilon is generally needed for floating point comparisons, isn't it in my case? http://stackoverflow.com/a/253874/3765753 — user3765753, Jun 23 '14 at 17:18
No. All you care about is that the `log` of `x*x - sqrt(x)` makes sense. And all you need for that is `x*x - sqrt(x) > 0`. I'd also suggest using `cbrt` instead of your `pow` wrapper. — tmyklebu, Jun 23 '14 at 17:22
@tmyklebu thanks, I forgot that cbrt was already in C99. The thing with epsilon is that it isn't needed here because 0.0 is easily stored and compared with other `double`s? — user3765753, Jun 23 '14 at 17:27
Note: rather than twice calculating `(pow(x, 2.0) - sqrt(x))`, do it once and save as temp result, say `t`. Due to subtle rules in C about allowing to go to higher precision, the 2 calculations may result in _slightly_ different intermediate values. Best to insure the value used for `t < 0` has exactly the same sign used in `ln()`. — chux - Reinstate Monica, Jul 08 '14 at 04:10

score 2 · Answer 1 · answered Jun 23 '14 at 20:29

double x=0.2 is not exactly 1/5, it's exactly 0.200000000000000011102230246251565404236316680908203125

Thus, when you repeatedly add this number you cumulate small round off errors.

Assuming there is a strict application of IEEE 754 double precision operations (without intermediate extra precision), your loop would provide those double (I used the shortest decimal representation that would round to same double assuming correct rounding to nearest even):

-5.0 -4.8 -4.6 -4.3999999999999995 -4.199999999999999 -3.999999999999999 -3.799999999999999 -3.5999999999999988 -3.3999999999999986 -3.1999999999999984 -2.9999999999999982 -2.799999999999998 -2.599999999999998 -2.3999999999999977 -2.1999999999999975 -1.9999999999999976 -1.7999999999999976 -1.5999999999999976 -1.3999999999999977 -1.1999999999999977 -0.9999999999999978 -0.7999999999999978 -0.5999999999999979 -0.39999999999999786 -0.19999999999999785 2.1649348980190553e-15 0.20000000000000218 0.4000000000000022 0.6000000000000022 0.8000000000000023 1.0000000000000022 1.2000000000000022 1.4000000000000021 1.600000000000002 1.800000000000002 2.000000000000002 2.2000000000000024 2.4000000000000026 2.6000000000000028 2.800000000000003 3.000000000000003 3.2000000000000033 3.4000000000000035 3.6000000000000036 3.800000000000004 4.0000000000000036 4.200000000000004 4.400000000000004 4.600000000000004 4.800000000000004 5.000000000000004 5.200000000000005 5.400000000000005 5.600000000000005 5.800000000000005 6.000000000000005 6.2000000000000055 6.400000000000006 6.600000000000006 6.800000000000006 7.000000000000006 7.200000000000006 7.400000000000007 7.600000000000007 7.800000000000007 8.000000000000007 8.200000000000006 8.400000000000006 8.600000000000005 8.800000000000004 9.000000000000004 9.200000000000003 9.400000000000002 9.600000000000001 9.8 10.0

You can improve by using

for(int i=-25; i<=50; ++i) {
    double x = i * 0.2;

But it's not perfect, 0.2 is not represented exactly as we saw above, and you get:

-5.0 -4.800000000000001 -4.6000000000000005 -4.4 -4.2 -4.0 -3.8000000000000003 -3.6 -3.4000000000000004 -3.2 -3.0 -2.8000000000000003 -2.6 -2.4000000000000004 -2.2 -2.0 -1.8 -1.6 -1.4000000000000001 -1.2000000000000002 -1.0 -0.8 -0.6000000000000001 -0.4 -0.2 0.0 0.2 0.4 0.6000000000000001 0.8 1.0 1.2000000000000002 1.4000000000000001 1.6 1.8 2.0 2.2 2.4000000000000004 2.6 2.8000000000000003 3.0 3.2 3.4000000000000004 3.6 3.8000000000000003 4.0 4.2 4.4 4.6000000000000005 4.800000000000001 5.0 5.2 5.4 5.6000000000000005 5.800000000000001 6.0 6.2 6.4 6.6000000000000005 6.800000000000001 7.0 7.2 7.4 7.6000000000000005 7.800000000000001 8.0 8.200000000000001 8.4 8.6 8.8 9.0 9.200000000000001 9.4 9.600000000000001 9.8 10.0

If you try with

for(int i=-25; i<=50; ++i) {
    double x = i / 5.0;

It's the best that you can get with floating point, because you get the nearest floating point approximation to your multiples of 1/5:

-5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0

Calculating the function and checking its domain

1 Answers1

Linked