Bug of the sqrt() compiled by DMD

[Johan]

On Sunday, 19 June 2022 at 08:36:33 UTC, Salih Dincer wrote:
>
> ```d
> import std.math, std.stdio;
> alias integer = ulong;
>
> void main()
> {
>   integer n = 94_906_260;
>   for( ; n <= 94_906_270; ++n)
>   {
>     integer root;
>     double root2, n2 = n * n;
> ```

For some of the values of `n` that you are using, a `double` 
cannot represent `n*n`.

94,906,260 * 94,906,260 = 9,007,198,187,187,600 = 0x1F FFFF C05E 
6D90

94,906,267 * 94,906,267 = 9,007,199,515,875,289 = 0x20 0000 0F90 
97D9

Both results are 56 bit integers.
An IEEE 64-bit double (D's `double`) has 52 bits of mantissa. It 
_can_ represent the first number without loss of precision 
because of the 4 zero bits at the end (56 - 4 = 52). It can _not_ 
represent the second number: instead of 0x20 0000 0F90 97D9, `n2` 
will store the value 0x20 0000 0F90 97D0, and you calculate the 
sqrt of that.

LDC creates instructions that do the sqrt in `double` domain, but 
DMD converts the `double` to a `real` (80bit floating point on 
x86) and does the sqrt in the `real` domain. Probably that is why 
it seems to work correctly for LDC and not for DMD 
(representation differences of the sqrt result for double and 
real, plus rounding when going back to `ulong`).

-Johan