Thread overview
Unexpected results with doubles
Jan 07, 2019
Joseph Malle
Jan 07, 2019
H. S. Teoh
Jan 07, 2019
Joseph Malle
Jan 07, 2019
H. S. Teoh
January 07, 2019
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



January 07, 2019
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...
January 07, 2019
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...
January 07, 2019
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??