Trace shows Linux has worse error than Windows in 1/(-1.)^1.2

Question

So I questioned why there are differences between Wolfram Cloud and Windows 11 Evaluate $\int_{0}^{\pi/2} \cos^a (x)\sin(ax) dx$ using Mathematica

E^(-((4 I \[Pi])/5)) Beta[-1, -(6/5), 11/5] // N

I FullSimplified as that gets better results (even though it should not, this is school math, there difference just happens somehow when -0.8090169943749475/0.587785252292473 changes to -0.8090169943749473/0.5877852522924732 with no reason whatsoever, but that is not what this question is about):

E^(-((4 I \[Pi])/5)) Beta[-1, -(6/5), 11/5] //FullSimplify//N

So, the result on windows is -6.6974 - 8.88178*10^-16 I (actually true value is -6.697403197294095 - 8.881784197001252-16 I) vs result on linux cloud that is -6.6974 - 1.33227*10^-15 I.

Why do you think that happens?

Well, I did trace as here: https://mathematica.stackexchange.com/a/271955/82985

and got the following result after 20 minutes of debug. Apparently it gets to

1/(-1.)^1.2 and not only lowers accuracy (why?) and thus loses one 1 in the complex part (on windows and linux), but it also makes a mistake in the real part (only on Linux).

That, that is why. You are aware this exists, right? https://core-math.gitlabpages.inria.fr/

Libm comparisions also exist: https://members.loria.fr/PZimmermann/papers/accuracy.pdf

Trace[-(-1)^(1/5) Beta[-1, -(6/5), 11/5] // N, TraceInternal -> True, 
 TraceDepth -> 10]

will give this on my Windows 11: https://pastebin.com/h9espi8G (wolfram cloud: https://pastebin.com/DGHQzaQX).

P.S. So there are 3 bugs here: on Linux rounding for floating precision (that cannot be controlled by N as in for 1/(-1.)^1.2) is horrible, even on Windows the error is just off-by-one in last digit, maybe because of AMD/Microsoft libm. Still an error. 2nd bug is that the general solver has somehow variable precision, not equal to this floating default and yet the precision in complex part is lower on one digit (or it just gets an error and sees 0 not 1). 3rd bug is that on windows it can sometimes just copy values through trace wrong with no reason (maybe it is some precision autotracking, dunno).

P.S.2 Opened a thread in community: https://community.wolfram.com/groups/-/m/t/2851685

Answer 1

This made me curious.

I work under macos on Apple Silicon. I get these results:

a = E^(-((4 I \[Pi])/5)) Beta[-1, -(6/5), 11/5] // N
b = E^(-((4 I \[Pi])/5)) Beta[-1, -(6/5), 11/5] // FullSimplify // N
Abs[Re[a] - Re[b]]/Abs[Re[b]]
Abs[Im[a] - Im[b]]/Abs[Im[b]]
Abs[a - b]/Abs[b]

-6.697403197294095 - 1.7763568394002505*^-15 I

-6.697403197294095 - 8.881784197001252*^-16 I

0.

1.

1.3261534262398443*^-16

So it looks as if we had 100% relative error in the imaginary part. But as a complex number, the relative error is almost as low as it can be. I agree, in a perfect world this should not happen, as double complex numbers have separate matissas and exponents for real and imaginary part.

Let's see where this error might come from; let's check the "complicated" functions first:

a1 = E^(-((4 I \[Pi])/5)) // N
b1 = E^(-((4 I \[Pi])/5)) // FullSimplify // N
Abs[Re[a1] - Re[b1]]/Abs[Re[b1]]
Abs[Im[a1] - Im[b1]]/Abs[Im[b1]]
Abs[a1 - b1]/Abs[b1]

-0.8090169943749473 - 0.5877852522924732 I

-0.8090169943749475 - 0.5877852522924731 I

1.3723111286221163*^-16

1.8888242266968722*^-16

1.5700924586837752*^-16

So both results are basically as good as it can get with double precision.

a2 = Beta[-1, -(6/5), 11/5] // N
b2 = Beta[-1, -(6/5), 11/5] // FullSimplify // N
Abs[Re[a2] - Re[b2]]/Abs[Re[b2]]
Abs[Im[a2] - Im[b2]]/Abs[Im[b2]]
Abs[a2 - b2]/Abs[b2]

5.418313004792032 - 3.936634828025925 I

5.418313004792032 - 3.936634828025925 I

Of course, the FullSimplify cannot do anything here.

Okay, but we have

a - a1 a2
b - b1 b2

0.+0.I

0.+0.I

So where the heck does the error in a - b come from? Just from the multiplication? Indeed! The simple "naive" way to compute the products would be

Re[a] == Re[a1] Re[a2] - Im[a1] Im[a2]
Im[a] == Re[a1] Im[a2] + Im[a1] Re[a2]

And apprently, that's what Mathematica does:

Re[a1] Re[a2] - Im[a1] Im[a2] - Re[a]
Re[a1] Im[a2] + Im[a1] Re[a2] - Im[a]
Re[b1] Re[b2] - Im[b1] Im[b2] - Re[b]
Re[b1] Im[b2] + Im[b1] Re[b2] - Im[b]

The problem with this approach is catastrophic cancellation in the imaginary parts:

Re[a1] Im[a2]
Im[a1] Re[a2]

3.1848044765212715

-3.1848044765212733

Re[b1] Im[b2]
Im[b1] Re[b2]

3.184804476521272

-3.184804476521273

As you can see the two summands are of almost equal magnitudes, but with opposite signs.

So it seems that this can be blamed solely on Mathematica's product routine for complex function. But I am pretty sure, this would also happen in C or C++ with the -Ofast or -ffast-math flags...

Answer 2

I'm not sure I understand the problem, despite attempts get clarification.

So far I'm not sure there is much more to this than to the following, this-is-school-math, behavior:

Sqrt[3. (4./3. - 1.) - 1.]
(*  0. + 1.49012*10^-8 I  *)

But in this instance, this behavior is predicted. In the OP's case, there is the added variability of different systems using different libraries on different CPUs.

About the variability, my Mac M1 Pro and the online Linux compute different values for $1/z$ for some complex values of $z$, and neither produces a correctly rounded result all the time. So far, all I've witnessed is incorrect rounding (<1ulp error).

If we consider complicated expressions (that is, involving more than one arithmetic operation), such as the OP's 1/(-1.)^1.2 that arises in Trace[N@Beta[-1, -(6/5), 11/5], TraceInternal -> True], we have to consider error propagation and additional rounding errors at each step. Even if operations are correctly rounded, the errors get propagated, and subtractive cancellation may lead to large relative error. One may also have correlated errors that cancel each other out. Different formulas have different numerical properties, even if they are school-math equivalent. For instance (-1.)^(-1.2) is a better way to compute 1/(-1.)^1.2, if you have a good Power function. (Possibly the Beta[] algorithm should be criticized for this.)

Now, interesting things happen in evaluating 1/(-1.)^1.2 on the Mac M1 and on Linux. The Mac computes z1 = (-1.)^1.2 correctly rounded and Linux does not; and neither the Mac nor Linux compute 1/z1 correctly rounded; however, on Linux, the two rounding errors cancel each other out, so that the result of 1/(-1.)^1.2 ends up correctly rounded. One might criticize the complex arithmetic of both systems for division not being correctly rounded; however, if complex division is done in real double-precision, then it is a complicated operation (in the above sense) and entails multiple roundings and error propagation.

That said, neither system computes 1/(-1)^(6/5) correctly rounded. That should be no surprise, since 6/5 is not the same as 1.2, since 1.2 equals 5404319552844595/4503599627370496 exactly, in binary double precision. (This applies to 11/5 and 2.2 in the Beta[] calculation; indeed, the error is 4 times as great, due to the way the bits line up in the binary expansion of the fractions. This rounding error is then propagated through the computation of Beta[].)

It is not made clear in the OP what is so horrible, but I infer it is that the imaginary part of a real result is nonzero in the first expression in the OP, whether simplified or not. I have not looked into how expensive it is to make complex arithmetic be correctly rounded, and that seems to be a significant source of error. (The error in the Beta[] cannot be completely traced.) The sources cited seem to consider real arithmetic and functions only. IEEE 754 does not treat complex functions (e.g. "$pow (x, y)$ signals the invalid operation exception for finite $x < 0$ and finite non-integer $y$"). I think most applications can tolerate errors of a couple ulp. This seems inevitable to me in computations that involve several operations such 1/(-1.)^1.2 or the Beta[] one (unless it is possible to write a special routine for your problem). Mathematica provides arbitrary precision for applications in which such small errors are too large.

Finally, I will point out that it is easy to write a numerically better routine in Mathematica for the OP's initial problem:

myBeta[x_, a_, b_] = 
  MeijerGReduce[Beta[x, a, b], x] // Activate // FunctionExpand;
Block[{Beta = myBeta},
 E^(-((4 I π)/5)) Beta[-1, -(6/5), 11/5] // FullSimplify // N
 ]
(*  -6.6974  *)

Appendix

What $1.2$ and $2.2$ equal in binary floating-point

Binary expansions of 1.2, 2.2. The color squares indicate a 1, the white squares indicate a 0, and the gray squares indicate a nonexistant or insignificant bit.

Answer 3

Okay, so the answer is that on Windows it is possible to use Intel C/C++ math.h libm library, which is much better than glibc on windows, indeed accuracy.pdf paper above shows that pow is just 0.817 for glibc (and gcc has its own pow, will not use glibc one) compared to 0.515 ulp for Intel C++ compiler library version 2023.0. And again, here we have complex numbers, that paper does not actually check those. That is worse than core-math, which is ideal 0.5 ulp.

I would love to clarify this and have https://gitlab.inria.fr/core-math/core-math/-/blob/master/src/binary32/pow/powf.c to be used on Mathematica on Linux.