Bug 14786 – The built-in exponentiation operator ^^ sometimes returns a value with the wrong sign.
Status
RESOLVED
Resolution
FIXED
Severity
minor
Priority
P1
Component
phobos
Product
D
Version
D2
Platform
x86_64
OS
Linux
Creation time
2015-07-07T23:18:00Z
Last change time
2016-01-03T14:15:09Z
Assigned to
nobody
Creator
thomas.bockman
Comments
Comment #0 by thomas.bockman — 2015-07-07T23:18:25Z
Negative one raised to any odd power should yield negative one.
However, when using 80-bit real arithmetic, -1 raised to a power greater than 2^^63 and less than 2^^64 yields +1.
This is despite the fact that integers in this range can still (just barely) be represented without loss of precision by an 80-bit real; even and odd can still be distinguished.
** Test case that currently fails **
void main(string[] args)
{
real b = -1;
real e = long.max;
e += 2;
assert(e % 2 == 1);
real r = b ^^ e;
assert(r == 1); // Passes, but shouldn't
assert(r == -1); // Fails, but shouldn't
e = ulong.max;
assert(e % 2 == 1);
r = b ^^ e;
assert(r == 1); // Passes, but shouldn't
assert(r == -1); // Fails, but shouldn't
}
I have marked this as a D runtime issue since it fails with all of DMD, GDC, and LDC. However, I couldn't actually find where the ^^ is defined in my quick search, so I don't know what the fix should look like. I suppose it might even be a hardware errata in my FPU (Xeon E3-1225 v3)...
Comment #1 by ibuclaw — 2015-07-08T07:48:38Z
Moving from druntime to phobos (std.math.pow is the component, unless the compiler was able to fold the operation).
Comment #2 by ibuclaw — 2015-07-08T08:04:03Z
This seems to be the correct extract from std.math to look at:
---
5675│ // Result is real only if y is an integer
5676│ // Check for a non-zero fractional part
5677├> if (y > -1.0 / real.epsilon && y < 1.0 / real.epsilon)
5678│ {
5679│ long w = cast(long)y;
5680│ if (w != y)
5681│ return sqrt(x); // Complex result -- create a NaN
5682│ if (w & 1) sign = -1.0;
5683│ }
5684│ x = -x;
---
If 'e = long.max', it is representable, and the (w & 1) check sets sign to '-1.0'
As instead 'e = long.max + 2', it is not representable, and so the signed-ness is not adjusted.
Checked result against libm routines, and yes it looks like a bug on our side.
Comment #3 by ibuclaw — 2015-07-08T17:50:39Z
OK, it looks like std.math.pow is only implemented to work with numbers only within the following boundaries.
-1.0/real.epsilon < X < 1/real.epsilon
I *could* relax this, but we'll need to test raising some random numbers to the power of an integer between long.max and ulong.max.
Comment #4 by thomas.bockman — 2015-07-08T19:42:14Z
I haven't tested it extensively yet, but here is a fixed version of that section:
// Result is real only if y is an integer
// Check for a non-zero fractional part
enum real maxPrecise = ulong.max;
enum real minPrecise = -maxPrecise;
if (y >= minPrecise && y <= maxPrecise)
{
real absY = abs(y);
ulong w = cast(ulong)absY;
if (w != absY)
return sqrt(x); // Complex result -- create a NaN
if (w & 1) sign = -1.0;
}
x = -x;
However, looking at the whole function, I think it would benefit from a major refactoring as well. So much duplicated logic...
Is that something I should do, or would it mostly likely be rejected on the basis of, "If it ain't broke, don't fix it"?
Comment #5 by ibuclaw — 2015-07-08T21:27:31Z
(In reply to thomas.bockman from comment #4)
> I haven't tested it extensively yet, but here is a fixed version of that
> section:
>
> // Result is real only if y is an integer
> // Check for a non-zero fractional part
> enum real maxPrecise = ulong.max;
> enum real minPrecise = -maxPrecise;
> if (y >= minPrecise && y <= maxPrecise)
> {
> real absY = abs(y);
> ulong w = cast(ulong)absY;
> if (w != absY)
> return sqrt(x); // Complex result -- create a NaN
> if (w & 1) sign = -1.0;
> }
> x = -x;
>
> However, looking at the whole function, I think it would benefit from a
> major refactoring as well. So much duplicated logic...
>
There is no duplicated logic. ;)
Comment #6 by thomas.bockman — 2015-08-31T19:01:38Z