www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Unexpected results with doubles

reply Joseph Malle <malle umich.edu> writes:
I am learning D.  I was working on Project Euler #199 (possible 
spoilers) and got some unexpected results.  It's probably a 
stupid mistake but I can't see it.

Here is an edited function from my program:

auto radius(const double r1, const double r2, const double r3) {
   auto const k1 = 1/r1;
   auto const k2 = 1/r2;
   auto const k3 = 1/r3;
   writeln();
   writeln("1   ", [k1, k2, k3]);
   writeln("2   ", [k1 * k2, k2 * k3, k3 * k1]);
   writeln("3   ", [k1 * k2 + k2 * k3 + k3 * k1]);
   assert(!isNaN(k1 * k2 + k2 * k3 + k3 * k1));
   writeln("4   ", [sqrt(k1 * k2 + k2 * k3 + k3 * k1)]);
   assert(!isNaN(sqrt(k1 * k2 + k2 * k3 + k3 * k1)));
   auto rv = 1 / (k1 + k2 + k3 + 2.0 * sqrt(k1 * k2 + k2 * k3 + k3 
* k1));
   assert(!isNaN(rv));
   writeln("radius ", [r1, r2, r3], " => ", rv);
   writeln();
   return rv;
}

Here is some output:

1   [41.7846, 6.4641, 6.4641]
2   [270.1, 41.7846, 270.1]
3   [581.985]
4   [24.1244]
radius [0.0239323, 0.154701, 0.154701] => 0.00971237


1   [41.7846, 6.4641, 6.4641]
2   [270.1, 41.7846, 270.1]
3   [581.985]
4   [nan]
radius [0.0239323, 0.154701, 0.154701] => 0.00971237


1   [41.7846, 6.4641, 6.4641]
2   [270.1, 41.7846, 270.1]
3   [581.985]
4   [nan]
radius [0.0239323, 0.154701, 0.154701] => 0.00971237

The "4   [nan]" is unexpected.  Each time they have the same 
input/same output.  But sometimes the 4th line is nan and 
sometimes it's not.  The asserts never fail.  I've seen this 
unexpected nan a few times with other inputs for this function.


If I change it to:

auto x = sqrt(k1 * k2 + k2 * k3 + k3 * k1);
writeln("4   ", [x]);
assert(!isNaN(x));

Then the assert fails.  I checked if the assert fails before the 
writeln too (as a sanity check) and yes, x is always NaN it seems.

I am doing $dmd -run on the command line.  Working on a 
reasonably up to date Mac.

$ dmd --version
DMD64 D Compiler v2.083.1
Copyright (C) 1999-2018 by The D Language Foundation, All Rights 
Reserved written by Walter Bright
Jan 07
parent reply "H. S. Teoh" <hsteoh quickfur.ath.cx> writes:
On Mon, Jan 07, 2019 at 07:57:14PM +0000, Joseph Malle via Digitalmars-d-learn
wrote:
[...]
 auto radius(const double r1, const double r2, const double r3) {
   auto const k1 = 1/r1;
   auto const k2 = 1/r2;
   auto const k3 = 1/r3;
   writeln();
   writeln("1   ", [k1, k2, k3]);
   writeln("2   ", [k1 * k2, k2 * k3, k3 * k1]);
   writeln("3   ", [k1 * k2 + k2 * k3 + k3 * k1]);
   assert(!isNaN(k1 * k2 + k2 * k3 + k3 * k1));
   writeln("4   ", [sqrt(k1 * k2 + k2 * k3 + k3 * k1)]);
   assert(!isNaN(sqrt(k1 * k2 + k2 * k3 + k3 * k1)));
   auto rv = 1 / (k1 + k2 + k3 + 2.0 * sqrt(k1 * k2 + k2 * k3 + k3 * k1));
   assert(!isNaN(rv));
   writeln("radius ", [r1, r2, r3], " => ", rv);
   writeln();
   return rv;
 }
 
 Here is some output:
 
 1   [41.7846, 6.4641, 6.4641]
 2   [270.1, 41.7846, 270.1]
 3   [581.985]
 4   [24.1244]
 radius [0.0239323, 0.154701, 0.154701] => 0.00971237
 
 
 1   [41.7846, 6.4641, 6.4641]
 2   [270.1, 41.7846, 270.1]
 3   [581.985]
 4   [nan]
 radius [0.0239323, 0.154701, 0.154701] => 0.00971237
 
 
 1   [41.7846, 6.4641, 6.4641]
 2   [270.1, 41.7846, 270.1]
 3   [581.985]
 4   [nan]
 radius [0.0239323, 0.154701, 0.154701] => 0.00971237
 
 The "4   [nan]" is unexpected.  Each time they have the same input/same
 output.  But sometimes the 4th line is nan and sometimes it's not.  The
 asserts never fail.  I've seen this unexpected nan a few times with other
 inputs for this function.
Either there's memory corruption somewhere, or there's a codegen bug in the compiler. Or the compiler somehow is malfunctioning with -run. Did you try compiling the program separately and running it? Does that make a difference?
 If I change it to:
 
 auto x = sqrt(k1 * k2 + k2 * k3 + k3 * k1);
 writeln("4   ", [x]);
 assert(!isNaN(x));
 
 Then the assert fails.  I checked if the assert fails before the
 writeln too (as a sanity check) and yes, x is always NaN it seems.
[...] The way to dig into the cause is to disassemble the radius() function and post the disassembly here. Then we can take a look to find out what's going on. What are the values of k1, k2, k3? T -- I think the conspiracy theorists are out to get us...
Jan 07
parent reply Joseph Malle <malle umich.edu> writes:
On Monday, 7 January 2019 at 20:18:28 UTC, H. S. Teoh wrote:
 Either there's memory corruption somewhere, or there's a 
 codegen bug in the compiler.
I think that must be it. I ran it with ldc instead of dmd and it worked fine (solved the original problem! woohoo).
 Or the compiler somehow is  malfunctioning with -run.  Did
 you try compiling the program separately and running it?
 Does that make a difference?
I did try doing it without -run and it had the same issue.
 The way to dig into the cause is to disassemble the radius() 
 function and post the disassembly here.  Then we can take a 
 look to find out what's going on.
What's the best way to get readable disassembly from dmd? I tried godbolt but there was lots more stuff than I expected...
Jan 07
parent "H. S. Teoh" <hsteoh quickfur.ath.cx> writes:
On Mon, Jan 07, 2019 at 09:42:42PM +0000, Joseph Malle via Digitalmars-d-learn
wrote:
 On Monday, 7 January 2019 at 20:18:28 UTC, H. S. Teoh wrote:
 Either there's memory corruption somewhere, or there's a codegen bug
 in the compiler.
I think that must be it. I ran it with ldc instead of dmd and it worked fine (solved the original problem! woohoo).
 Or the compiler somehow is  malfunctioning with -run.  Did you try
 compiling the program separately and running it?  Does that make a
 difference?
I did try doing it without -run and it had the same issue.
 The way to dig into the cause is to disassemble the radius()
 function and post the disassembly here.  Then we can take a look to
 find out what's going on.
What's the best way to get readable disassembly from dmd? I tried godbolt but there was lots more stuff than I expected...
Hmm. On Linux I'd use objdump: objdump -D myprogram | sed -ne'/radius/,/^$/p' > radius.disas This extracts the assembly dump of the radius() function into the file 'radius.disas'. Not sure how to do this on a Mac. (It should be possible to install objdump on Mac, but I don't own a Mac so I've no idea.) T -- Why do conspiracy theories always come from the same people??
Jan 07