The standard <cmath>'s function pow() fails to produce some sub-normalized numbers if the version that takes an integer exponent is used (using G++ 4.4.3, x86_64-linux-gnu). For example, the following piece of code will cause 'p' to get a zero value, instead of the expected sub-normalized value 5.56268464626e-309. double p=pow(2.0, -1024); I dug the include files to try to discover the reason for this strange result. It seemed to me that the function that gets called when the exponent is integer first checks if the exponent is negative, and calls yet another overloaded version that takes an unsigned exponent, passing to it the absolute value of the original exponent, and dividing 1.0 by the received result. So, apparently my original call 'pow(2.0, -1024)' became translated to '1.0/pow(2.0, 1024u)'. However, as 2**1024 is too large for a double, it becomes the positive infinity, and the result of dividing 1.0 by this infinity is zero. As a sanity check, I tried to make the exponent a double (i.e. I did double p=pow(2.0, -1024.0); instead). Not surprisingly, now I got 5.56268464626e-309, which was the result that I expected from the beginning with an integer exponent. POSSIBLE FIX: Maybe the fix is simple (yet I didn't work it to a very long extent): instead of making 'pow(some_double, some_negative_int)' become '1.0/pow(some_double, unsigned(-some_negative_int))', change it to 'pow(1.0/some_double, unsigned(-some_negative_int))'.
Integer exponent powers are computed using TYPE NAME (TYPE x, int m) { unsigned int n = m < 0 ? -m : m; TYPE y = n % 2 ? x : 1; while (n >>= 1) { x = x * x; if (n % 2) y = y * x; } return m < 0 ? 1/y : y; } or optimal power series expansion for constant exponents starting with GCC 4.5. ISTR the standard allows this. You should use a double exponent to use the real power function, pow (2., -1024.)
(In reply to comment #1) > Integer exponent powers are computed using > > TYPE > NAME (TYPE x, int m) > { > unsigned int n = m < 0 ? -m : m; > TYPE y = n % 2 ? x : 1; > while (n >>= 1) > { > x = x * x; > if (n % 2) > y = y * x; > } > return m < 0 ? 1/y : y; > } > > or optimal power series expansion for constant exponents starting with > GCC 4.5. > > ISTR the standard allows this. You should use a double exponent to > use the real power function, pow (2., -1024.) I have no objection to using such optimizations for integer exponents, and I am aware that it . But perhaps it could be improved further, by widening the range over which the result for an integer exponent will be consistent with the corresponding double (or whatever) exponent. If the code above was changed to something like what goes below, would it be worse? TYPE NAME (TYPE x, int m) { unsigned int n; if(m < 0){ n = -m; x = TYPE(1) / x; } TYPE y = n % 2 ? x : 1; while (n >>= 1) { x = x * x; if (n % 2) y = y * x; } return y; }
(In reply to comment #2) > (In reply to comment #1) > > Integer exponent powers are computed using > > > > TYPE > > NAME (TYPE x, int m) > > { > > unsigned int n = m < 0 ? -m : m; > > TYPE y = n % 2 ? x : 1; > > while (n >>= 1) > > { > > x = x * x; > > if (n % 2) > > y = y * x; > > } > > return m < 0 ? 1/y : y; > > } > > > > or optimal power series expansion for constant exponents starting with > > GCC 4.5. > > > > ISTR the standard allows this. You should use a double exponent to > > use the real power function, pow (2., -1024.) > > I have no objection to using such optimizations for integer exponents, and I am > aware that it . But perhaps it could be improved further, by widening the > range over which the result for an integer exponent will be consistent with the > corresponding double (or whatever) exponent. > > If the code above was changed to something like what goes below, would it be > worse? > > TYPE > NAME (TYPE x, int m) > { > unsigned int n; > if(m < 0){ > n = -m; > x = TYPE(1) / x; > } > TYPE y = n % 2 ? x : 1; > while (n >>= 1) > { > x = x * x; > if (n % 2) > y = y * x; > } > return y; > } I'm sure you can create a similar issue for this variant. Fact is that there are multiple roundings involved in computing integer powers.
I'm following with interest the numerical analysis side of the discussion, but I have to add that, *as a libstc++-v3 proper* issue, as mentioned by Richard, the issue is moot, because in currently active release branches either the library uses a front-end builtin (__builtin_powi*) or just pow unconditionally (in C++0x mode). In other terms, *minimal* open-coding in the library for these functions (and we mean to stick to this general design choice, AFAIK)
I don't think this is a library issue. Maybe Bill is interested / has an opinion..
(In reply to comment #5) > I don't think this is a library issue. Maybe Bill is interested / has an > opinion.. Just to agree that it is WAD. The proposed solution just moves the values of x and m for which the rounding will cause edge-case anomalies. Either way, the compounded rounding error means the two versions of pow will never produce identical results. When precision is important, use the real-valued version.