www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.ldc - std.math unittests - accuracy of floating point results

reply "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.

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  [1]
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

[1] 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.
Mar 01 2015
next sibling parent reply "Kevin Brogan" <kevin brogan.ca> writes:
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  [1]
 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

 [1] 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.
Mar 01 2015
parent reply "Kai Nacke" <kai redstar.de> writes:
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  [1]
 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
Mar 01 2015
parent reply "Johan Engelen" <j j.nl> writes:
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
Mar 08 2015
parent Dan Olson <zans.is.for.cans yahoo.com> writes:
"Johan Engelen" <j j.nl> writes:


 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


Thanks back.  I'll have to try it out on ARM.
Mar 09 2015
prev sibling parent Dan Olson <zans.is.for.cans yahoo.com> writes:
"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
Mar 01 2015