Fix Issue 14786 - std.math.pow sometimes gets the sign of the result wrong.

[tsbockman]

This is a simple fix for Phobos Issue #14786.

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 cases that fails without this PR **/
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
}

Thanks to @ibuclaw for pointing me to the relevant code.

[burner]

these three lines look strange. aren't you comparing y with y in the end. isn't the only logical reason for that been != floating point imprecision

[tsbockman]

@burner The entire current implementation is definitely strange. It's sub-optimal performance-wise, and whoever worked on it last left a comment on the most important part which begins TODO: This is not accurate in practice...

Really I think the whole thing should probably be re-written, but for this pull request I decided to make the minimum change necessary to fix the bug I found. Anyway, after myself staring at the part you highlighted for a while, I eventually remembered why it's there:

It's not checking for floating-point imprecision, but rather whether y is an exact integer value - conceptually, if y == floor(y). It's expressed a little less clearly than that in the code because I need to compute cast(ulong) fabs(y), anyway, for the even/odd check immediately following. So, this way is (infinitesimally) faster.

The roots of negative numbers almost always have an imaginary component. So, a complex number would generally be needed to represent the result of pow(x, y) when (x < 0) && (y != floor(y)). When this case is detected, sqrt(x) is used to generate an imaginary-style NaN value to return.

[tsbockman]

Anyone want to review this? @9il or @ibuclaw perhaps?

[9il]

Any condition must be removed. We need to implement funcitons

bool floatingPointIsIntegral(F)(const F x) if(isFloatingPoint!F)
bool floatingPointIsOdd     (F)(const F x) if(isFloatingPoint!F)

first.

[tsbockman]

@9il I can see that there are some things in std.math - particularly in this function - that would benefit from refactoring.

But, do I really have to start that refactoring just to apply this tiny fix?

[9il]

@tsbockman Half of std.math requires to be rewritten because different kinds of reasons. So I don't know what exactly you what to refactor. FloatTraits can be used to write floatingPointIsIntegral and floatingPointIsOdd.

[9il]

Looks like current fix with ulong will fail anyway for integer pow greater then ulong.max

[tsbockman]

@9il

FloatTraits can be used to write floatingPointIsIntegral and floatingPointIsOdd.

That looks helpful. I will try re-writing this patch to use floatTraits.

Looks like current fix with ulong will fail anyway for integer pow greater then ulong.max

It works as intended for everything up to 80-bit reals. However, it would fail for 128-bit quads or double-doubles.

According to the language specification, 80-bit is largest floating point type currently supported by D.

Is this incorrect? Is real larger than 80 bits on some D platforms?

[9il]

According to the language specification, 80-bit is largest floating point type currently supported by D.

Is this incorrect?

Yes.

Is real larger than 80 bits on some D platforms?

No. However std.math should contain generic code that can be easily ported to 128 bit real.

[tsbockman]

@9il

No. However std.math should contain generic code that can be easily ported to 128 bit real.

I notice that there are references to ibmExtended format in std.math. What is that referring to?

The only thing my web search turned up that fits the name, is this weird format: Extended-precision 128-bit. However, the floatTraits numbers don't seem to match Wikipedia's description.

[dlang-bot]

Fix	Bugzilla	Description
✓	14786	The built-in exponentiation operator ^^ sometimes returns a value with the wrong sign.

[ibuclaw]

The only thing my web search turned up that fits the name, is this weird format

IBM 128-bit extended precision format is a pair of IEEE double precision numbers whose sum is equal to the extended precision value. The number with greater magnitude is first. This format has the same magnitude range as an IEEE double precision value, but effectively 106 bits of significand precision.

[ibuclaw]

IBM 128-bit extended precision format is a pair of IEEE double precision numbers

Or if you prefer an answer in code format, what you have is the following supported natively.

struct real
{
    double high;
    double low;
}

See https://www-01.ibm.com/support/knowledgecenter/ssw_aix_61/com.ibm.aix.genprogc/128bit_long_double_floating-point_datatype.htm

[tsbockman]

@ibuclaw @9il

See https://www-01.ibm.com/support/knowledgecenter/ssw_aix_61/com.ibm.aix.genprogc/128bit_long_double_floating-point_datatype.htm

What an ugly hack! There is no way that any attempt to support that format in std.math will succeed without some way to actually test it.

This format has ... effectively 106 bits of significand precision.

Per the docs you linked, this is incorrect:

The actual number of bits of precision can vary. If the low-order part is much less than 1 ULP of the high-order part, significant bits (either all 0's or all 1's) are implied between the significant of the high-order and low-order numbers. Certain algorithms that rely on having a fixed number of bits in the significand can fail when using 128-bit long double numbers.
...
The number of bits in the significand is not fixed, but for a correctly formatted number (except in the denormal range) the minimum number available is 106.

I have updated this PR to document its limitation to 64-bit precision with a static assert. Can we please just pull this? It's a strict improvement over what's currently shipping, and can always be replaced with a more flexible approach later.

Making everyone keep using the current broken implementation of pow until I properly support 128-bit types makes no sense, seeing as nothing else in std.math works with them, anyway.

floatTraits does not, and cannot with its current design, properly implement ibmExtended. I have some ideas about how to properly support uncommon floating point formats in the future, but that is a much larger project to do properly.

[9il]

What an ugly hack!

No need to support ibmExtended, just IEEE floating point standard.
https://en.wikipedia.org/wiki/IEEE_floating_point

Can we please just pull this?

Yes, but unittest for the largest and smallest odd pows, largest and smallest even pows should be added, e.g. for maxPrecise, maxPrecise - 1, minPrecise, minPrecise + 1

[tsbockman]

@9il

Yes, but unittest for the largest and smallest odd pows, largest and smallest even pows should be added, e.g. for maxPrecise, maxPrecise - 1, minPrecise, minPrecise + 1

I already had some unittests, but I realized they were only valid with 80-bit Intel extended real, and would have failed with a lower precision type, such as double. This is now fixed. The tests also now use a named constant, maxOdd, to make the intent more explicit.

I trimmed a couple of redundant instructions from the runtime code, while I was at it.

[ibuclaw]

No need to support ibmExtended, just IEEE floating point standard.

That's irrelevant to this PR though. The implementation should work regardless of underlying float. :-)

[ibuclaw]

I already had some unittests, but I realized they were only valid with 80-bit Intel extended real

Use wolfram to get the results for the tests if you haven't already done so. I cannot stress how important is to use a reliable source for unittests vs. whatever number returns when you run it on your machine.

[tsbockman]

@ibuclaw

That's irrelevant to this PR though. The implementation should work regardless of underlying float. :-)

The issue is that this code:

enum ulong_dig = 8 * ulong.sizeof;
static assert(real.mant_dig <= ulong_dig);
enum real maxOdd = ulong.max >> (ulong_dig - real.mant_dig);

const real absY = fabs(y);
if (absY <= maxOdd)
{
    const ulong w = cast(ulong)absY;
    if (w != absY)
        return sqrt(x); // Complex result -- create a NaN
    if (w & 1)
        sign = -1.0;
}
x = -x;

Is considerably faster than its format-agnostic equivalent, because floating-point modulus is slow:

// Untested; this is just an example
const real ym2 = fabs(y) % 2.0L;
if(ym2 == 1.0L)
    sign = -1.0;
else if(ym2 != 0.0L)
    return sqrt(x);  // Complex result -- create a NaN
x = -x;

I could use a static if to select which implementation to use based on the real format, but then the generic version would just be dead - and therefore untested - code, because the ulong version is preferred on all current D platforms.

I can do this if you really want, but I don't think intentionally shipping dead code is a good idea.

[tsbockman]

@ibuclaw

Use wolfram to get the results for the tests if you haven't already done so. I cannot stress how important is to use a reliable source for unittests vs. whatever number returns when you run it on your machine.

This is not that kind of math - my unittest basically just says, "Negative one to an odd power should return negative one."

My error was in hard-coding the maximum odd value, without taking into account the possibility of other real formats. I corrected this by using real.mant_dig to calculate the maximum odd value at compile time.

[ibuclaw]

I'm fine with having a fast path for 80-bit reals.

I did some further testing (also because I'm a skeptic of cast(ulong)real), and this method I can vouch to work fine on a few targets (only testing 80bit and 64bit reals though). Because it leverages libc/builtin functions it should be fine with whatever real type you throw at it, and maybe even get a small speed boost (targets of GDC may replace C math function calls with intrinsics).

Of course, this should be considered slower than cast(ulong) because of the use of floor.

if (x < 0)
{
    real w = floor(y);
    if (w != y)
        return sqrt(x); // Complex result -- create a NaN

    // For negative x, find out if the integer exponent is odd or even.
    w = ldexp(y, -1);
    w = floor(w);
    w = ldexp(w, 1);
    if (w != y)
        sign = -1.0;

    x = fabs(x);
}

[tsbockman]

@ibuclaw

I'm fine with having a fast path for 80-bit reals.

The ulong method is the fast path for any floating-point format with mant_dig <= 64, not just for 80-bit. That means any alternative generic method will be dead code on all supported platforms, and at high risk of bit-rot.

I did some further testing (also because I'm a skeptic of cast(ulong)real), and this method I can vouch to work fine on a few targets (only testing 80bit and 64bit reals though). Because it leverages libc/builtin functions it should be fine with whatever real type you throw at it, and maybe even get a small speed boost (targets of GDC may replace C math function calls with intrinsics).

Of course, this should be considered slower than cast(ulong) because of the use of floor.

I was able to optimize your method further:

if(x < 0)
{
    if(floor(y) != y)
        return sqrt(x); // Complex result -- create a NaN

    // For negative x, find out if the integer exponent is odd or even.
    const hy = scalbn(y, -1);
    if(floor(hy) != hy)
        sign = -1.0;

    x = fabs(x);
}

However, in my micro-benchmark even this streamlined version is still over 4x slower than the ulong method (on 64-bit) - no faster than the naive abs(y) % 2.0L == 1.0L approach.

EDIT: In case anyone wants to know, I tested primarily with GDC. The gap is smaller with DMD, but the conclusion is the same. The LDC results are more ecclectic...

[tsbockman]

Since there seems to be strong desire for a generic solution, I should probably mention that I have a proof-of-concept which is by far the fastest in my testing, beating even the ulong approach:

bool floatingPointIsRound(F)(const F x, const int expOf2 = 0)
    if(isFloatingPoint!F)
{
    alias FP = FloatParts!F;
    FP xp; xp.num = x;
    const xe = xp.exp;
    if(xe == FP.traits.expNaN)
        return false;

    const fde2 = FP.traits.fracDig + expOf2;
    if(xe >= fde2)
        return true;
    if(xe < expOf2)
        return xp.mant == 0;

    const remMask = ~(~cast(xp.Mant)0 << (fde2 - xe));
    static if(FP.traits.isWholeExplicit)
        return (xp.mant & remMask) == 0;
    else
        return (xp.frac & remMask) == 0;
}

This code should work correctly and very quickly for any vaguely normal floating-point format (i.e., not ibmExtended).

However, in order to express the algorithm so directly and cleanly, I had to re-write floatTraits and create a FloatParts templated union:

module floatparts;

import std.traits : Unqual, isFloatingPoint;
import std.math : RealFormat;

template floatTraits(F)
    if(isFloatingPoint!F)
{
    alias Format = RealFormat;

    enum mantDig = F.mant_dig;
    enum fracDig = mantDig - 1;
    static if(mantDig == 11) {
        enum format = Format.ieeeHalf;
        enum usedBytes = 2;
        alias MantInt = ushort;
        enum mantOffset = 0;
        enum isWholeExplicit = false;
        enum expDig = 5;
        enum expOffset = fracDig;
    } else static if(mantDig == 24) {
        enum format = Format.ieeeSingle;
        enum usedBytes = 4;
        alias MantInt = uint;
        enum mantOffset = 0;
        enum isWholeExplicit = false;
        enum expDig = 8;
        enum expOffset = fracDig;
    } else static if(mantDig == 53 && F.sizeof == 8) {
        enum format = Format.ieeeDouble;
        enum usedBytes = 8;
        alias MantInt = ulong;
        enum mantOffset = 0;
        enum isWholeExplicit = false;
        enum expDig = 11;
        enum expOffset = fracDig;
    } else static if(mantDig == 53 && F.sizeof == 12) {
        enum format = Format.ieeeExtended53;
        enum usedBytes = 10;
        alias MantInt = ulong;
        enum mantOffset = 64 - mantDig;
        enum isWholeExplicit = true;
        enum expDig = 15;
        enum expOffset = 64;
    } else static if(mantDig == 64) {
        enum format = Format.ieeeExtended;
        enum usedBytes = 10;
        alias MantInt = ulong;
        enum mantOffset = 0;
        enum isWholeExplicit = true;
        enum expDig = 15;
        enum expOffset = mantDig;
    } else static if(mantDig == 113) {
        enum format = Format.ieeeQuadruple;
        enum usedBytes = 16;
        alias MantInt = ucent;
        enum mantOffset = 0;
        enum isWholeExplicit = false;
        enum expDig = 15;
        enum expOffset = fracDig;
    } else
        static assert(false, "Unknown floating-point format for " ~ F.stringof);

    enum mantMax = ~cast(MantInt)0 >> (MantInt.sizeof*8 - mantDig);
    enum fracMax = ~cast(MantInt)0 >> (MantInt.sizeof*8 - fracDig);
    enum expBias = ~0 >>> ((int.sizeof*8) - (expDig - 1));
    enum expMin = -expBias;
    enum expMax = expBias;
    enum expNaN = expMax + 1;
}

union FloatParts(F)
    if(isFloatingPoint!F)
{
public:
    alias traits = floatTraits!F;
    alias Mant = traits.MantInt;

    F num = F.nan;
    alias num this;
    /+this(F num) {
        this.num = num; }+/

    @property Mant frac() const {
        return (mantBits & fracMask) >> traits.mantOffset; }
    @property void frac(Mant newFrac) {
        mantBits = (mantBits & ~fracMask) | ((newFrac << traits.mantOffset) & fracMask); }

    static if(traits.isWholeExplicit) {
        @property Mant whole() const {
            return (mantBits & wholeMask) >> wholeOffset; }
        @property void whole(Mant newWhole) {
            mantBits = (mantBits & ~wholeMask) | ((newWhole << wholeOffset) & wholeMask); }

        @property Mant mant() const {
            return (mantBits & mantMask) >> traits.mantOffset; }
        @property void mant(Mant newMant) {
            mantBits = (mantBits & ~mantMask) | ((newMant << traits.mantOffset) & mantMask); }

        alias explicit = mant;
    } else {
        @property uint whole() const {
            return (lastBits & expMask) != 0; }

        @property Mant mant() const {
            return (cast(Mant)whole << wholeOffset) | frac; }

        alias explicit = frac;
    }

    @property int exp() const {
        return ((lastBits & expMask) >> expShift) - traits.expBias; }
    @property void exp(int newExp) {
        lastBits = (lastBits & ~expMask) | (((newExp + traits.expBias) << expShift) & expMask); }

    @property int sign() const {
        return cast(short)(lastBits & 0x8000) >> 15; }
    @property void sign(int newSign) {
        lastBits = (lastBits & 0x7FFF) | (newSign & 0x8000); }

    @property bool isFinite() const {
        return (lastBits & expMask) != expMask; }
    @property bool isNormal() const {
        switch(lastBits & expMask) {
            case 0:
                return mantBits == 0;
            case expMask:
                return false;
            default:
                return !(traits.isWholeExplicit && (mantBits & wholeMask) == 0);
        }
    }
    @property bool isZero() const {
        static if(traits.usedBytes >= (Mant.sizeof + ushort.sizeof))
            return (mantBits == 0) && ((lastBits & ~0x8000) == 0);
        else {
            ptrdiff_t u16x = allBits.length - 1;
            if(lastBits == 0x8000)
                --u16x;
            while(u16x > 0) {
                if(allBits[u16x] != 0)
                    return false;
                --u16x;
            }
            return true;
        }
    }
    ref typeof(this) makeZero() {
        ushort keepSign = lastBits & 0x8000;
        mantBits = 0;
        lastBits = keepSign;
        return this;
    }
    ref typeof(this) makeZero(int newSign) {
        mantBits = 0;
        lastBits = (newSign & 0x8000);
        return this;
    }
    @property bool isNaN() const {
        return (lastBits & expMask) == expMask && frac != 0; }
    @property bool isSignaling() const {
        return (mantBits & quietMask) == 0; }
    ref typeof(this) makeNaN(bool quiet = true) {
        lastBits = lastBits | expMask;
        enum nanMant = wholeMask | (cast(Mant)1 << traits.mantOffset);
        if(quiet)
            mantBits = mantBits | quietMask | nanMant;
        else
            mantBits = (mantBits & ~quietMask) | nanMant;
        return this;
    }
    @property bool isInfinite() const {
        return (lastBits & expMask) == expMask && frac == 0; }
    ref typeof(this) makeInfinite() {
        lastBits = lastBits | expMask;
        mantBits = (mantBits & ~fracMask) | wholeMask;
        return this;
    }

private:
    enum Mant fracMask = traits.fracMax << traits.mantOffset;
    enum wholeOffset = traits.fracDig + traits.mantOffset;
    enum Mant wholeMask = traits.isWholeExplicit? cast(Mant)1 << wholeOffset : 0;
    enum Mant mantMask = wholeMask | fracMask;
    enum expShift = 15 - traits.expDig;
    enum expMask = 0x8000 - (1 << expShift);

    enum quietMask = cast(Mant)1 << (wholeOffset - 1);

    Mant mantBits;
    static assert(traits.usedBytes % 2 == 0);
    ushort[traits.usedBytes / 2] allBits;
    @property ushort lastBits() const {
        return allBits[$ - 1]; }
    @property void lastBits(ushort newBits) {
        allBits[$ - 1] = newBits; }
}

If there is interest, I could potentially finish this and refactor std.math to make use of it in place of the current ad-hoc bit twiddling.

[ibuclaw]

There are a few things I can see wrong with that. Feel free to raise a PR. Bearing in mind that this PR exists too: #3499

[tsbockman]

@ibuclaw

There are a few things I can see wrong with that.

No doubt - it's just a proof of concept at the moment.

Feel free to raise a PR. Bearing in mind that this PR exists too: #3499

I might do that, once my CheckedInt work is ready to merge/rot in the review queue.

[tsbockman]

@ibuclaw
What about this PR?

As I said, I am willing to add a slow, generic path (either % or ldexp/scalbn based). However, it will not actually be used on any current platform.

Is that what you want?

[tsbockman]

@ibuclaw @9il
I went ahead and added the generic ldexp method, in theory enabling support for 128-bit types.

Normally it does not even get compiled, of course, but I have tried to test it myself by temporarily commenting out the faster ulong method. It passed all Phobos unittests on my machine.

Good enough?

[9il]

Normally it does not even get compiled, of course, but I have tried to test it myself by temporarily commenting out the faster ulong method. It passed all Phobos unittests on my machine.

Looks good

[ibuclaw]

Wouldn't this static if condition always be false? The only way it could ever succeed is if CTFE started supporting ucent.

[9il]

ulong.max converted to real

[tsbockman]

@ibuclaw

Wouldn't this static if condition always be false?

maxOdd is of type real. If a more precise real type (such as IEEE quad) is ever implemented on some platform, the condition would be true.

I added that fall-back algorithm because you said,

The implementation should work regardless of underlying float.

[ibuclaw]

maxOdd is of type real.

Yep, I saw that after I posted. :-)

[tsbockman]

@ibuclaw Would you say this is ready to merge?

[ibuclaw]

I'm fine with this (sorry I've been busy with 2.067) (N.B. It looks like std.math is the only module that I have local patches for now - well done phobos devs!)

[tsbockman]

@9il You've already approved the latest version, yes? Anyone want to toggle auto-merge on?

[9il]

@9il You've already approved the latest version, yes?

Yes.

Anyone want to toggle auto-merge on?

ping @ibuclaw

[ibuclaw]

ping @ibuclaw

Typical. :-p

[ibuclaw]

Actually, sorry to be a bother, but curly braces on the next line. This format is used everywhere.

static if(...)
{
}
else
{
}

[9il]

... and rebase please.

[tsbockman]

@ibuclaw

Actually, sorry to be a bother, but curly braces on the next line. This format is used everywhere.

Oops. I knew about that rule; I just wasn't thinking about it last time.

@9il

... and rebase please.

OK.

[ibuclaw]

I'll re-test on my arm box just to be sure and will merge depending on when my last comment is addressed.

[tsbockman]

It should be fixed now. (I amended the last commit.)

[tsbockman]

@ibuclaw The corrected version has now passed all my tests, and the auto-tester.

[ibuclaw]

Auto-merge toggled on

[tsbockman]

Thanks.