Fix Issue 11652: support numerical ^^ complex operations in std.complex

[WebDrake]

This patch adds a Complex.opBinaryRight method for op == "^^" and numeric LHS.

The current design takes advantage of the fact that we know we're dealing with a number whose argument is either 0 or PI. An alternative design would be to just create a complex number whose real part is equal to the input, and then raise it to the power of this:

return typeof(return)(lhs, 0) ^^ this;

This would generate probably less precisely calculated values but would be 100% consistent with complex ^^ complex calculations. I leave the choice here up to review. :-)

[WebDrake]

Interesting, FreeBSD and Darwin seem to be failing for FP operations on the following unittests:

auto rer = a ^^ complex(2.0, 0.0);
assert(rer.re == a ^^ 2.0);
assert(rer.im == 0.0);

I know FP equality comparisons are risky, but I was rather hoping this one might be reliable ... oh well, back to approxEqual. :-P

[WebDrake]

To understand what I'm getting at with respect to precision of calculation, consider the case of -1 raised to the power of 0.5 + 2i. This should result in a complex number with absolute value of exp(-2 * PI) and argument of PI/2. That is, theoretically, a number whose real part is 0 (argument PI/2) and whose imaginary part is exp(-2 * PI).

Using the current opBinaryRight for numeric types to calculate (-1.0) ^^ complex(0.5, 2.0) one gets a value of

 -0.00000000000000000000005061713666775166342967400662379056270419250452749659588259
+ 0.00186744273170798887745425176376556919422000646591186523437500000000000000000000i

while if one calculates complex(-1.0, 0) ^^ complex(0.5, 2.0) one gets:

  0.00000000000000000011434411173245104255279602635219036750149656173942161420040153
+ 0.00186744273170798931113512075796734279720112681388854980468750000000000000000000i

So, the real part is further away from the "true" value of 0 by several orders of magnitude with the complex ^^ complex calculation. Further, if we look at absolute value and argument, we get:

Calculating exp(-2 * PI) and using the constant PI_2:

exp(-2 * PI) = 0.00186744273170798881403265983809983130647935922752367332577705383300781250000000
      PI / 2 = 1.57079632679489661925640447970309310221637133508920669555664062500000000000000000

numeric ^^ complex:

abs = 0.00186744273170798887745425176376556919422000646591186523437500000000000000000000
        (diff 6.34216e-20)
arg = 1.57079632679489655799898173427209258079528808593750000000000000000000000000000000
        (diff -6.12574e-17)

complex ^^ complex

abs = 0.00186744273170798931113512075796734279720112681388854980468750000000000000000000
         (diff 4.97102e-19)
arg = 1.57079632679489655799898173427209258079528808593750000000000000000000000000000000
         (diff -6.12574e-17)

[WebDrake]

I've pushed an updated patch to avoid the failing unittest. I've preserved the == 0.0 comparison just to see if it works; if not, I'll fix that up with approxEqual too.

[monarchdodra]

Do we really need to get the abs of lhs, if the first thing we do is check it's sign? At the very least, I'd only calculate it in the second branch (if at all).

Unless my floating point-fu is too low, shouldn't this work just fine?

        if (lhs >= 0)
        {
            return typeof(return)(lhs ^^ this.re, log(lhs) * this.im);
        }
        else
        {
            // theta = PI
            FPTemporary!R r = -lhs;
            FPTemporary!(CommonType!(T, R)) ab = r ^^ this.re * exp(-PI * this.im);
            FPTemporary!(CommonType!(T, R)) ar = PI * this.re + log(r) * this.im;
            return typeof(return)(ab * std.math.cos(ar), ab * std.math.sin(ar));
        }

[WebDrake]

Agree, I'll tweak.

[WebDrake]

This update should fix the errors identified for the internals of numeric ^^ complex and also address the absolute-value bits. I've reinstated some of the floating-point == floating-point unittests just to see how they roll now the internals are right; will force-push a corrected commit if they still fail.

It's a bit irritating that the cos/sin functions are so slightly inexact that complex ^^ complex with imaginary part zero results in a complex number whose imaginary part is almost but not quite zero. Floating point, we love you. :-)

[WebDrake]

Here's your extra unittest. I'll squash all of this when it's done, how is the rest of it looking?

[WebDrake]

.... and the exact-equality unittest fails on 32-bit systems still. Oh well. Fixed, should pass everything now.

[monarchdodra]

Seems good to me overall. Give it a good squashing, it it'll be good to go.

It's LGTM for me. But I wouldn't mind another reviewer giving it a thumbs up first before.

[WebDrake]

Will do :-) It would be nice to get Lars or Don's input as the original module authors.

[monarchdodra]

Will do :-) It would be nice to get Lars or Don's input as the original module authors.

Ping them: @kyllingstad , @donc .

[kyllingstad]

Well, the unittest which is currently failing (line 574) should pass in its current form, I think. That is, it should be safe to assume exact equality in that case.

[WebDrake]

I agree, but it's not obvious why. It's failing on 32-bit systems (or at least, 32-bit Darwin and FreeBSD). What I could do, if you guys don't mind the extra testing weight, is to push a version that writefln's the values to screen so that we can check them.

[WebDrake]

Let's see how these go. For the cases using (-a) it makes some sense as maybe cos and sin of PI are being calculated inexactly, but it is bizarre where we're taking powers of (positive) a.

[WebDrake]

Well, the results here look rather bizarre. The real parts really do look like they should be identical to the results of the numeric ^^ numeric operation. The imaginary parts for (-a) ^^ complex are slightly non-zero which goes with expectations that sin(n * PI) might deviate very slightly from its expected value of 0. The deviation for cos(n * PI) looks like it must be small enough to be unobservable.

Seeing this, I can't explain the equality failures that have been observed, so I'm going to reinstate equality comparisons where appropriate and see what happens. Meanwhile, does anyone else have any idea why I should see the results I have here but get equality failures?

[monarchdodra]

https://d.puremagic.com/test-results/pull.ghtml?projectid=1&runid=811862&logid=6

Are these not replayed on win 64?

Meanwhile, does anyone else have any idea why I should see the results I have here but get equality failures?

Honestly, floating point equality in D is... strange. In particular, at the end of the day, if you are comparing something smaller than real, there is basically 0 guarantee that a comparison will succeed, even on equal objects, as one may be stored in a double, and the other, (during the comparison) temporarily in a more precise real.

http://d.puremagic.com/issues/show_bug.cgi?id=8475
http://d.puremagic.com/issues/show_bug.cgi?id=8476

Relevant (I believe):
Notice that even with pure, the == will fail...

double getFloat() pure
{
   return sqrt(1.1);
}

void main()
{
    writeln(isIdentical(getFloat(), getFloat())); //true
    writeln(getFloat() is getFloat()); //true
    writeln(getFloat() == getFloat()); //False, wtf?
}

and

Maxim Fomin 2012-10-16 09:18:02 PDT

Regardless of whether (when comparing floats with reals) compiler should
compare with full or truncated to smaller type precision, current behavior
seems to be inconsistent with spec

The conclusion I came to is that unless a floating point represents an exact value, == is simply not reliable, even for exactly equal floating point types

[WebDrake]

Are these not replayed on win 64?

Well, they should be. If they're not, that probably means something about how Win64 handles output from the unittests.

Honestly, floating point equality in D is... strange.

So I see. Well, in that case I don't see what else I can do, except perhaps to switch to using isIdentical. I would have assumed that isIdentical and == ought to be the same for floating-point types, so I can only surmise it's some fundamental quirk of optimization in how == is implemented.

[WebDrake]

As per Ali Çehreli's advice on D.learn, let's try outputting those results with %a formatting to be sure of equality.

[WebDrake]

OK, so we're definitely dealing with a bug with == rather than a problem with the code. By the way, I think the other bug report you meant to link to was this: http://d.puremagic.com/issues/show_bug.cgi?id=8745

[WebDrake]

A note on why I just pushed another commit with == comparisons: I realized that isIdentical fails even on my 64-bit system if you do something like,

assert(isIdentical(rer.re, a ^^ 2.0));

which I guess could be put down to the fact that whereas rer.re is a double, a ^^ 2.0 is a value on the FPU, and ironically, whereas an equality comparison works out fine, an exact bitwise comparison may fail.

So, I'm guessing temporary variables are essential anyway, so I might as well give == one last chance.

[WebDrake]

Typical that the exact test run I most want to confirm the success or failure of, instead falls over on a thread priority error :-P

[WebDrake]

Darwin32 seems to be OK with the unittests as they have become (if not with thread priority...), so squashing. Assuming the tests go all right, is this OK with everyone? :-)

[WebDrake]

OK, I give up. Equality == FP comparisons are obviously not something that can be reliably used here, at least for 32-bit systems. Further, isIdentical comparisons fail even on 64-bit systems: e.g. on my own system, the unittest fails if we replace line 571 with assert(isIdentical(rer.re, rcheck));.

So, unless anyone can suggest an alternative, I think we have to concede the case to approxEqual. I take it everyone agrees that this is a fault of FP comparison rather than any issue with the code?

[WebDrake]

Last try before I revert to approxEqual. (And lo and behold, thread priority error strikes again :-P )

[WebDrake]

OK, looks like finally we have working unittests. Sorry for the noise and truckload of patches over what ought in principle to have been very simple.

Anyway, @monarchdodra, @kyllingstad, how does this look?

[monarchdodra]

Nitpick: The if seems a little verbose, it could be reduced to simply one of:

//Result ≡ PI (mod 2PI)
assert(approxEqual(std.math.abs(arg(rec1b)), PI, EPS)));
assert(approxEqual(std.math.cos(arg(rec1b)), -1, EPS)));
assert(approxEqual(std.math.sin(arg(rec1b)), 0, EPS)));

Furthermore, if you want to be very anal, I think the if could be wrong, in the sense that depending on where the rounding occurs in the call chain, the first result could be actually be -PI, whereas the second result could be PI

Doing this should be the simplest and most correct form, I think?

auto arg1b = arg(rec1b);
assert(approxEqual(arg1b, PI, EPS) || approxEqual(arg1b, -PI, EPS));

[WebDrake]

Furthermore, if you want to be very anal, I think the if could be wrong

Ahh, I see the concern. I don't mind trialling the abs of the argument; the reason I formulated it like this was so that code coverage analysis would let you see if there was any change in whether the imaginary part was turning out slightly negative or not.

Doesn't the simplest and most correct form you propose potentially fall victim to the same issue you identify with the use of if?

[monarchdodra]

Doesn't the simplest and most correct form you propose potentially fall victim to the same issue you identify with the use of if?

AFAIK, no. Since once stored in a variable, the rounding error will be bound to that variable's neighborhood, and in no way can it jump from PI to -PI. My fear with calling a function twice in a row would be a -1e19 being rounded (or not) to 0, prior to the call making the result of the arg call jump from -PI to PI. With a variable, there is only 1 call, so no problem.

I think.

[WebDrake]

Ahh, sorry, I missed the separate variable. :-) Patch coming up!

[WebDrake]

Here you go :-) Let's see whether the thread priority demon strikes this time ...

[monarchdodra]

Auto-merge toggled on

[WebDrake]

We have auto-merge now? Superb! When did that happen?

[monarchdodra]

We have auto-merge now? Superb! When did that happen?

Roughly about a month ago? It's pretty nice yeah. Unfortunalty, darwin's ridiculously high failure rate does get into its way :/