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