Bug of the sqrt() compiled by DMD

[kdevel]

On Sunday, 19 June 2022 at 20:26:27 UTC, Salih Dincer wrote:
> On Sunday, 19 June 2022 at 10:41:32 UTC, Johan wrote:
>> 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.
>
> I think the problem is not with bits.

It is not debatable that 94906267² = 9007199515875289 has no 
exact representation in the double type:

```
import std.stdio;

union Double {
    double d;
    ulong u;
}

void main ()
{
    Double d = Double (94906267);
    d.d *= 94906267;
    for (int i = -4; i < 5; ++i) {
       Double e = d;
       e.d += i;
       writefln!"%4d %19.2f %14x" (i, e.d, e.u);
    }
}
```

prints

```
   -4 9007199515875284.00 4340000007c84bea
   -3 9007199515875284.00 4340000007c84bea
   -2 9007199515875286.00 4340000007c84beb
   -1 9007199515875288.00 4340000007c84bec
    0 9007199515875288.00 4340000007c84bec
    1 9007199515875288.00 4340000007c84bec
    2 9007199515875290.00 4340000007c84bed
    3 9007199515875292.00 4340000007c84bee
    4 9007199515875292.00 4340000007c84bee
```

> Because the size of the number is only twice uint.max.

The mantissa of double has only 53 bits which is less than twice 
the number of bits in an (u)int. 94906267 occupies 27 bits:

```
                                |                          |
            94906267 
0000000000000000000000000000000000000101101010000010011110011011
    9007199515875289 
0000000000100000000000000000000000001111100100001001011111011001
```

> When integer type is changed qua long, the problem is fixed.

I would say that the problem is revealed. IMHO the result should 
not depend on wether the long type is signed or not. dmd seems to 
generate different code for the signed/unsigned cases:

```
import std.math;
ulong r1 (double d)
{
    return cast (ulong) sqrt (d);
}

long r2 (double d)
{
    return cast (long) sqrt (d);
}
```

dmd -O -c / objdump -Cdarw shows what goes on:
```
Disassembly of section .text.ulong sub.r1(double):

0000000000000000 <ulong sub.r1(double)>:
    0:   55                      push   %rbp
    1:   48 8b ec                mov    %rsp,%rbp
    4:   48 83 ec 10             sub    $0x10,%rsp
    8:   f2 0f 11 45 f0          movsd  %xmm0,-0x10(%rbp)
    d:   dd 45 f0                fldl   -0x10(%rbp)
   10:   d9 fa                   fsqrt
   12:   b9 00 00 00 80          mov    $0x80000000,%ecx
   17:   c7 45 f0 00 00 00 00    movl   $0x0,-0x10(%rbp)
   1e:   89 4d f4                mov    %ecx,-0xc(%rbp)
   21:   c7 45 f8 3e 40 bf 0f    movl   $0xfbf403e,-0x8(%rbp)
   28:   db 6d f0                fldt   -0x10(%rbp)
   2b:   d8 d9                   fcomp  %st(1)
   2d:   df e0                   fnstsw %ax
   2f:   d9 7d fc                fnstcw -0x4(%rbp)
   32:   d9 6d fa                fldcw  -0x6(%rbp)
   35:   f6 c4 01                test   $0x1,%ah
   38:   74 18                   je     52 <ulong 
sub.r1(double)+0x52>
   3a:   db 6d f0                fldt   -0x10(%rbp)
   3d:   de e9                   fsubrp %st,%st(1)
   3f:   df 7d f0                fistpll -0x10(%rbp)
   42:   48 8b 45 f0             mov    -0x10(%rbp),%rax
   46:   48 c1 e1 20             shl    $0x20,%rcx
   4a:   48 03 c1                add    %rcx,%rax
   4d:   d9 6d fc                fldcw  -0x4(%rbp)
   50:   eb 0a                   jmp    5c <ulong 
sub.r1(double)+0x5c>
   52:   df 7d f0                fistpll -0x10(%rbp)
   55:   48 8b 45 f0             mov    -0x10(%rbp),%rax
   59:   d9 6d fc                fldcw  -0x4(%rbp)
   5c:   c9                      leaveq
   5d:   c3                      retq
         ...

Disassembly of section .text.long sub.r2(double):

0000000000000000 <long sub.r2(double)>:
    0:   55                      push   %rbp
    1:   48 8b ec                mov    %rsp,%rbp
    4:   48 83 ec 10             sub    $0x10,%rsp
    8:   f2 0f 11 45 f0          movsd  %xmm0,-0x10(%rbp)
    d:   dd 45 f0                fldl   -0x10(%rbp)
   10:   d9 fa                   fsqrt
   12:   dd 5d f0                fstpl  -0x10(%rbp)
   15:   f2 0f 10 45 f0          movsd  -0x10(%rbp),%xmm0
   1a:   f2 48 0f 2c c0          cvttsd2si %xmm0,%rax
   1f:   c9                      leaveq
   20:   c3                      retq
   21:   00 00                   add    %al,(%rax)
         ...

```

The unsigned part