Invalid results from exp2 function (C math library)

[pcauchon]

When I try to run this code on my machine (iMac with MacOS 10.5), I
get very strange results. I am using this compiler:
i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5465)

The first program loop gives the results I expect (the value of 2**x
when x is smaller than -1024 should be smaller than 10**-308).
However the second loop gives values that seem to be corrupted
somehow. Is it specific to iMac? Is there any way to detect when the
returned value is corrupted?

Many thanks
Philippe

---- begin code ----
#include <math.h>
#include <stdio.h>

int main(int argc, const char * argv[]) {
double x, y;

printf("first loop\n");
for (x = 2040; x < 2060; ++x) {
y = exp2(-x);
printf("%.4f %.4f\n", x, y);
}

printf("second loop\n");
for (x = 2040.1313; x < 2060; ++x) {
y = exp2(-x);
printf("%.4f %.4f\n", x, y);
}

return 0;
}
--- end code ----

Program output on my machine:
first loop
2040.0000 0.0000
2041.0000 0.0000
2042.0000 0.0000
2043.0000 0.0000
2044.0000 0.0000
2045.0000 0.0000
2046.0000 0.0000
2047.0000 0.0000
2048.0000 0.0000
2049.0000 0.0000
2050.0000 0.0000
2051.0000 0.0000
2052.0000 0.0000
2053.0000 0.0000
2054.0000 0.0000
2055.0000 0.0000
2056.0000 0.0000
2057.0000 0.0000
2058.0000 0.0000
2059.0000 0.0000
second loop
2040.1313 0.0000
2041.1313 0.0000
2042.1313 0.0000
2043.1313 0.0000
2044.1313 0.0000
2045.1313 -3.6482
2046.1313 -2.0000
2047.1313 -2.0000
2048.1313 -2.0000
2049.1313 -2.0000
2050.1313 -2.0000
2051.1313 -233.7301
2052.1313 -467.4603
2053.1313 -934.9206
2054.1313 -1869.8412
2055.1313 -3739.6823
2056.1313 -7479.3646
2057.1313 -14958.7292
2058.1313 -29917.4584
2059.1313 -59834.9169

[bigcaterpillar]

The value x implies exp2(x) has overflow. You can use exp(x) to try.

[pcauchon]

On Jun 7, 2:52*pm, "[e-mail]"
<[e-mail]:

The value x implies exp2(x) has overflow. You can use exp(x) to try.

When I try with the standard exp function I do get normal results:
exp(-x) is always equal to 0 when x is larger than about 750.

I get the issue I described earlier only with the exp2 function.

[user923005]

On Jun 6, 10:13*pm, [e-mail] wrote:

When I try to run this code on my machine (iMac with MacOS 10.5), I
get very strange results. I am using this compiler:
i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5465)

The first program loop gives the results I expect (the value of 2**x
when x is smaller than -1024 should be smaller than 10**-308).
However the second loop gives values that seem to be corrupted
somehow. Is it specific to iMac? Is there any way to detect when the
returned value is corrupted?

Many thanks
Philippe

---- begin code ----
#include <math.h>
#include <stdio.h>

int main(int argc, const char * argv[]) {
* double x, y;

* printf("first loop\n");
* for (x = 2040; x < 2060; ++x) {
* * y = exp2(-x);
* * printf("%.4f %.4f\n", x, y);
* }

* printf("second loop\n");
* for (x = 2040.1313; x < 2060; ++x) {
* * y = exp2(-x);
* * printf("%.4f %.4f\n", x, y);
* }

* return 0;}

--- end code ----

Program output on my machine:
first loop
2040.0000 0.0000
2041.0000 0.0000
2042.0000 0.0000
2043.0000 0.0000
2044.0000 0.0000
2045.0000 0.0000
2046.0000 0.0000
2047.0000 0.0000
2048.0000 0.0000
2049.0000 0.0000
2050.0000 0.0000
2051.0000 0.0000
2052.0000 0.0000
2053.0000 0.0000
2054.0000 0.0000
2055.0000 0.0000
2056.0000 0.0000
2057.0000 0.0000
2058.0000 0.0000
2059.0000 0.0000
second loop
2040.1313 0.0000
2041.1313 0.0000
2042.1313 0.0000
2043.1313 0.0000
2044.1313 0.0000
2045.1313 -3.6482
2046.1313 -2.0000
2047.1313 -2.0000
2048.1313 -2.0000
2049.1313 -2.0000
2050.1313 -2.0000
2051.1313 -233.7301
2052.1313 -467.4603
2053.1313 -934.9206
2054.1313 -1869.8412
2055.1313 -3739.6823
2056.1313 -7479.3646
2057.1313 -14958.7292
2058.1313 -29917.4584
2059.1313 -59834.9169

If your exponent exceeds DBL_MAX_10_EXP (decimal) or DBL_MAX_EXP
(binary) or if your exponent is smaller than DBL_MIN_10_EXP (decimal)
or DBL_MIN_EXP (binary) then you are sure to overflow / underflow.
After that happens, you have a total loss of significance.

Here is what happens when I run your program:
C:\tmp>mtest
first loop

exp2 underflow error
2040.0000 0.0000

exp2 underflow error
2041.0000 0.0000

exp2 underflow error
2042.0000 0.0000

exp2 underflow error
2043.0000 0.0000

exp2 underflow error
2044.0000 0.0000

exp2 underflow error
2045.0000 0.0000

exp2 underflow error
2046.0000 0.0000

exp2 underflow error
2047.0000 -1.#IND

exp2 underflow error
2048.0000 -1.#IND

exp2 underflow error
2049.0000 -1.#IND

exp2 underflow error
2050.0000 -1.#IND

exp2 underflow error
2051.0000 -1.#IND

exp2 underflow error
2052.0000 -1.#IND

exp2 underflow error
2053.0000 -1.#IND

exp2 underflow error
2054.0000 -1.#IND

exp2 underflow error
2055.0000 -1.#IND

exp2 underflow error
2056.0000 -1.#IND

exp2 underflow error
2057.0000 -1.#IND

exp2 underflow error
2058.0000 -1.#IND

exp2 underflow error
2059.0000 -1.#IND
second loop

exp2 underflow error
2040.1313 -1.#IND

exp2 underflow error
2041.1313 -1.#IND

exp2 underflow error
2042.1313 -1.#IND

exp2 underflow error
2043.1313 -1.#IND

exp2 underflow error
2044.1313 -1.#IND

exp2 underflow error
2045.1313 -1.#IND

exp2 underflow error
2046.1313 -1.#IND

exp2 underflow error
2047.1313 -1.#IND

exp2 underflow error
2048.1313 -1.#IND

exp2 underflow error
2049.1313 -1.#IND

exp2 underflow error
2050.1313 -1.#IND

exp2 underflow error
2051.1313 -1.#IND

exp2 underflow error
2052.1313 -1.#IND

exp2 underflow error
2053.1313 -1.#IND

exp2 underflow error
2054.1313 -1.#IND

exp2 underflow error
2055.1313 -1.#IND

exp2 underflow error
2056.1313 -1.#IND

exp2 underflow error
2057.1313 -1.#IND

exp2 underflow error
2058.1313 -1.#IND

exp2 underflow error
2059.1313 -1.#IND

[user923005]

On Jun 6, 11:09*pm, pcauchon <[e-mail]:

On Jun 7, 2:52*pm, "[e-mail]"

<[e-mail]:

The value x implies exp2(x) has overflow. You can use exp(x) to try.

When I try with the standard exp function I do get normal results:
exp(-x) is always equal to 0 when x is larger than about 750.

I get the issue I described earlier only with the exp2 function.

GIGO.
Has it occurred to you that exp(-x) is not equal to zero when x is
larger than 750?
It seems strange that it should make you happy.
You are asking your floating point library to do something it isn't
able to do.

#include <math.h>
#include <stdio.h>
int main(void)
{
double dexp = exp(-750);
double dback = log(dexp);
printf("%.20g should be 750.0 but is it?\n", dback);
return 0;
}

[pcauchon]

On Jun 7, 3:19*pm, user923005 <[e-mail]:

On Jun 6, 11:09*pm, pcauchon <[e-mail]:

On Jun 7, 2:52*pm, "[e-mail]"

<[e-mail]:

The value x implies exp2(x) has overflow. You can use exp(x) to try.

When I try with the standard exp function I do get normal results:
exp(-x) is always equal to 0 when x is larger than about 750.

I get the issue I described earlier only with the exp2 function.

GIGO.
Has it occurred to you that exp(-x) is not equal to zero when x is
larger than 750?
It seems strange that it should make you happy.
You are asking your floating point library to do something it isn't
able to do.

#include <math.h>
#include <stdio.h>
int * * * * * * main(void)
{
* * double * * * * *dexp = exp(-750);
* * double * * * * *dback = log(dexp);
* * printf("%.20g should be 750.0 but is it?\n", dback);
* * return 0;

}

I agree with you it's not symetric. Anyway very few floating point
operation have exact symetric inverse and that's a well-known
limitation.

In my specific case I only take exponent of negative values. I would
expect the result of exp(x) or exp2(x) for x <= 0 to be between 0 and
1. I do not really care if I have an absolute error less than say
10**-10, so computing 0 instead of the exact 10**-483 or whatever the
value is doesn't bother me.

I have problems when I get a negative result from exp2 or when the
result magnitude is larger than 1. From what I see I could get both
even if it clearly doesn't make any sense.

Any clue on how to detect when a floating point underflow occurs?

[user923005]

On Jun 6, 11:30*pm, pcauchon <[e-mail]:

On Jun 7, 3:19*pm, user923005 <[e-mail]:

On Jun 6, 11:09*pm, pcauchon <[e-mail]:

On Jun 7, 2:52*pm, "[e-mail]"

<[e-mail]:
The value x implies exp2(x) has overflow. You can use exp(x) to try.

When I try with the standard exp function I do get normal results:
exp(-x) is always equal to 0 when x is larger than about 750.

I get the issue I described earlier only with the exp2 function.

GIGO.
Has it occurred to you that exp(-x) is not equal to zero when x is
larger than 750?
It seems strange that it should make you happy.
You are asking your floating point library to do something it isn't
able to do.

#include <math.h>
#include <stdio.h>
int * * * * * * main(void)
{
* * double * * * * *dexp = exp(-750);
* * double * * * * *dback = log(dexp);
* * printf("%.20g should be 750.0 but is it?\n", dback);
* * return 0;

}

I agree with you it's not symetric. Anyway very few floating point
operation have exact symetric inverse and that's a well-known
limitation.

In my specific case I only take exponent of negative values. I would
expect the result of exp(x) or exp2(x) for x <= 0 to be between 0 and
1. I do not really care if I have an absolute error less than say
10**-10, so computing 0 instead of the exact 10**-483 or whatever the
value is doesn't bother me.

I have problems when I get a negative result from exp2 or when the
result magnitude is larger than 1. From what I see I could get both
even if it clearly doesn't make any sense.

I worry that you are inserting values larger than -2000 into exp2() in
the first place. That is not too different from inserting -1 into
sqrt() or tan(pi/2) or any other value that is not going to be
representible because the domain does not make sense for the given
hardware. If you have values that extreme, why are you taking
negative exponentials of them?

Once intermediate calculations underflow, you should not expect
subsequent calculation to make sense.
Consider my example of exp(log(x)). Mathematically, the answer is x,
but if I insert a huge or tiny x, then the answer will no longer have
any meaning. I expect that some intermediate calculation hurled
chunks when you told it to take the exponential of something smaller
than is possible for the hardware.

Any clue on how to detect when a floating point underflow occurs?

It depends on your compiler. On some compilers you can set floating
point exceptions to have implementation defined behavior.

[pete]

pc******@gmail.com wrote:

When I try to run this code on my machine (iMac with MacOS 10.5), I
get very strange results. I am using this compiler:
i686-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5465)

The first program loop gives the results I expect (the value of 2**x
when x is smaller than -1024 should be smaller than 10**-308).
However the second loop gives values that seem to be corrupted
somehow. Is it specific to iMac? Is there any way to detect when the
returned value is corrupted?

Many thanks
Philippe

---- begin code ----
#include <math.h>
#include <stdio.h>

int main(int argc, const char * argv[]) {
double x, y;

printf("first loop\n");
for (x = 2040; x < 2060; ++x) {
y = exp2(-x);
printf("%.4f %.4f\n", x, y);
}

printf("second loop\n");
for (x = 2040.1313; x < 2060; ++x) {
y = exp2(-x);
printf("%.4f %.4f\n", x, y);
}

return 0;
}
--- end code ----

Program output on my machine:
first loop
2040.0000 0.0000
2041.0000 0.0000
2042.0000 0.0000
2043.0000 0.0000
2044.0000 0.0000
2045.0000 0.0000
2046.0000 0.0000
2047.0000 0.0000
2048.0000 0.0000
2049.0000 0.0000
2050.0000 0.0000
2051.0000 0.0000
2052.0000 0.0000
2053.0000 0.0000
2054.0000 0.0000
2055.0000 0.0000
2056.0000 0.0000
2057.0000 0.0000
2058.0000 0.0000
2059.0000 0.0000
second loop
2040.1313 0.0000
2041.1313 0.0000
2042.1313 0.0000
2043.1313 0.0000
2044.1313 0.0000
2045.1313 -3.6482
2046.1313 -2.0000
2047.1313 -2.0000
2048.1313 -2.0000
2049.1313 -2.0000
2050.1313 -2.0000
2051.1313 -233.7301
2052.1313 -467.4603
2053.1313 -934.9206
2054.1313 -1869.8412
2055.1313 -3739.6823
2056.1313 -7479.3646
2057.1313 -14958.7292
2058.1313 -29917.4584
2059.1313 -59834.9169

It looks like a bug to me.

ISO/IEC 9899:1999 (E)
7.12.1 Treatment of error conditions

5 The result underflows if the magnitude of the mathematical result
is so small that the mathematical result cannot be represented,
without extraordinary roundoff error,
in an object of the specified type.195) If the result underflows,
the function returns an implementation-defined value whose magnitude
is no greater than the smallest normalized positive number in the
specified type;
if the integer expression math_errhandling & MATH_ERRNO is nonzero,
whether errno acquires the value ERANGE is implementation-defined;
if the integer expression math_errhandling & MATH_ERREXCEPT
is nonzero, whether the ‘‘underflow’’ floating-point exception is
raised is implementation-defined.

195) The term underflow here is intended to encompass both
‘‘gradual underflow’’ as in IEC 60559 and also
‘‘flush-to-zero’’ underflow.

--
pete