std.math unittests - accuracy of floating point results
Mar 01, 2015
Johan Engelen
Mar 02, 2015
Kevin Brogan
Mar 02, 2015
Kai Nacke
Mar 08, 2015
Johan Engelen
Mar 10, 2015
Dan Olson
Mar 02, 2015
Dan Olson
```Hi all,
I am working on making more of the std.math unittests pass (I'm new to the project, and it is a nice niche thing to tinker on, learning the codebase, workflow, etc.).
I've hit on a problem that I do not know how to handle: floating point comparison.

There are some tests that check whether exp(x) works well, including overflow checks for different x. See phobos/std/math.d line 2083. The checks are defined for 80-bit reals, and I am converting them to 64-bit reals (Win64). The problem is that the checks are bit-precise (i.e. assert(x == y)), but the calculation results are sometimes 1 ulp off. For example, the results of two tests:

std.math.E = 0x4005bf0a8b145769 = 2.7182818284590450
exp(1.0L)  = 0x4005bf0a8b145769 = 2.7182818284590450  
Wolfram Alpha                   = 2.718281828459045235...

E*E*E     = 0x403415e5bf6fb105 = 20.085536923187664
exp(3.0L) = 0x403415e5bf6fb106 = 20.085536923187668
Wolfram Alpha                  = 20.08553692318766774...

I do not know how I can make the second test pass, without breaking the first one. I feel the tests are too strict and should allow an error of 1 ulp.
dmd 2.066.1 passes these unittests with values corresponding Wolfram Alpha.

(Incidentally, an inaccuracy of 1 ulp also haunts a std.csv unittest, but I do not yet know why exactly)

How should I go about fixing these unittests for us?

Thanks,
Johan

 The correct result for exp(1.0L) I was able to obtain by enabling the LLVM intrinsic for exp, although there is a comment saying that that actually causes unittest failure. Without the LLVM intrinsic, exp(1.0L) is 1 ulp off.
```
```On Sunday, 1 March 2015 at 23:08:54 UTC, Johan Engelen wrote:
> Hi all,
>  I am working on making more of the std.math unittests pass (I'm new to the project, and it is a nice niche thing to tinker on, learning the codebase, workflow, etc.).

Guess what I'm doing. :-)

> I've hit on a problem that I do not know how to handle: floating point comparison.
>
> There are some tests that check whether exp(x) works well, including overflow checks for different x. See phobos/std/math.d line 2083. The checks are defined for 80-bit reals, and I am converting them to 64-bit reals (Win64). The problem is that the checks are bit-precise (i.e. assert(x == y)), but the calculation results are sometimes 1 ulp off. For example, the results of two tests:
>
> std.math.E = 0x4005bf0a8b145769 = 2.7182818284590450
> exp(1.0L)  = 0x4005bf0a8b145769 = 2.7182818284590450  
> Wolfram Alpha                   = 2.718281828459045235...
>
> E*E*E     = 0x403415e5bf6fb105 = 20.085536923187664
> exp(3.0L) = 0x403415e5bf6fb106 = 20.085536923187668
> Wolfram Alpha                  = 20.08553692318766774...
>
> I do not know how I can make the second test pass, without breaking the first one. I feel the tests are too strict and should allow an error of 1 ulp.

They are too strict. floating point math is not exact between different architectures, or even compilation flags. You can get a different result just because the compiler reordered two operations.

> dmd 2.066.1 passes these unittests with values corresponding Wolfram Alpha.
>
> (Incidentally, an inaccuracy of 1 ulp also haunts a std.csv unittest, but I do not yet know why exactly)
>
> How should I go about fixing these unittests for us?
>
> Thanks,
>   Johan
>
>  The correct result for exp(1.0L) I was able to obtain by enabling the LLVM intrinsic for exp, although there is a comment saying that that actually causes unittest failure. Without the LLVM intrinsic, exp(1.0L) is 1 ulp off.

Take a look at bool approxEqual in std.math. A lot of unit tests use it already, but some of them don't. Every floating point comparison should be using them, in my opinion.

```
```"Johan Engelen" <j@j.nl> writes:

> Hi all,
>  I am working on making more of the std.math unittests pass (I'm new
> to the project, and it is a nice niche thing to tinker on, learning
> the codebase, workflow, etc.).
> I've hit on a problem that I do not know how to handle: floating point
> comparison.

Hi Johan.  I think you are running into similar problems I ran into with 64-bit reals on ARM.  Hopefully there is a common solution for all 64-bit real architectures.

>
> (Incidentally, an inaccuracy of 1 ulp also haunts a std.csv unittest, but I do not yet know why exactly)

I encountered this too.  I found that std.conv.parse!double with 64-bit real is often off by 1 ulp compared to strod().  For 80-bit reals, they match.  I think the fix for this problem is a change to std.conv.parse.

> How should I go about fixing these unittests for us?

You could look in https://github.com/smolt/phobos/blob/ios/std/math.d and see if any changes I made for iOS may help you.  I have all remaining problems marked with D versions that are prefixed with "WIP_", like WIP_FloatPrecIssue.

Here are the phobos files possibly of interest:
std/csv.d
std/internal/math/errorfunction.d
std/internal/math/gammafunction.d
std/math.d

Hope this helps,
Dan
```
```Hi Kevin! Hi Johan!

On Monday, 2 March 2015 at 00:30:28 UTC, Kevin Brogan wrote:
> On Sunday, 1 March 2015 at 23:08:54 UTC, Johan Engelen wrote:
>> Hi all,
>> I am working on making more of the std.math unittests pass (I'm new to the project, and it is a nice niche thing to tinker on, learning the codebase, workflow, etc.).
>
> Guess what I'm doing. :-)
>
>> I've hit on a problem that I do not know how to handle: floating point comparison.
>>
>> There are some tests that check whether exp(x) works well, including overflow checks for different x. See phobos/std/math.d line 2083. The checks are defined for 80-bit reals, and I am converting them to 64-bit reals (Win64). The problem is that the checks are bit-precise (i.e. assert(x == y)), but the calculation results are sometimes 1 ulp off. For example, the results of two tests:
>>
>> std.math.E = 0x4005bf0a8b145769 = 2.7182818284590450
>> exp(1.0L)  = 0x4005bf0a8b145769 = 2.7182818284590450  
>> Wolfram Alpha                   = 2.718281828459045235...
>>
>> E*E*E     = 0x403415e5bf6fb105 = 20.085536923187664
>> exp(3.0L) = 0x403415e5bf6fb106 = 20.085536923187668
>> Wolfram Alpha                  = 20.08553692318766774...
>>
>> I do not know how I can make the second test pass, without breaking the first one. I feel the tests are too strict and should allow an error of 1 ulp.
>
> They are too strict. floating point math is not exact between different architectures, or even compilation flags. You can get a different result just because the compiler reordered two operations.

I agree, too. Bitwise comparison is worse if you are working with floating points.

BTW: These fixes should go upstream. Please create Phobos PRs for them. I am happy to cherry-pick these changes as soon as they committed.

Regards,
Kai
```
```Thanks for the link to your branch Dan.

A few 64-bit fixes for gammafunction.d :
https://github.com/D-Programming-Language/phobos/pull/3045

```
```"Johan Engelen" <j@j.nl> writes: