I want to check the largest integer that can be held in various floating-point types in C without loss of precision. Here is a test program:
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#define FLOATTYPE long double
#define ONE ((FLOATTYPE)1.0)
#define TWO ((FLOATTYPE)2.0)
int
main(int argc,char*argv[]){
int i;
FLOATTYPE x;
x = ONE;
for(i=0;;++i){
printf("1.0<<%3d: x=%.0Lf",i,(long double)x);
if((x+ONE)!=x &&
(x+ONE)- x == ONE){
printf(" ... can increment without loss of precision\n");
}else{
printf(" ... cannot increment without loss of precision\n");
break;
}
x *= TWO;
}
printf("FLT_RADIX = %d\n",FLT_RADIX);
printf("FLT_MANT_DIG = %d\n",FLT_MANT_DIG);
printf("DBL_MANT_DIG = %d\n",DBL_MANT_DIG);
printf("LDBL_MANT_DIG = %d\n",LDBL_MANT_DIG);
printf("\nsizeof(FLOATTYPE) = %lu\n",sizeof(x));
}
Some results (using gcc-9 (Ubuntu 9.4.0-1ubuntu1~16.04) 9.4.0):
When FLOATTYPE is float: sizeof is 4, and the loop exits at i==24, which equals FLT_MANT_DIG.
When FLOATTYPE is double: sizeof is 8, and the loop exits at i==53, which equals DBL_MANT_DIG.
When FLOATTYPE is __float128: sizeof is 16, and the loop exits at i==113.
They all make sense. However:
FLOATTYPE is long double: sizeof is 16, and the loop exits at i==53, which does not equal LDBL_MANT_DIG (which is 64).It seems like long double is taking more memory than double, but not giving increased precision. How come?
Edit: more detail on the compiler etc:
This is on a Windows 10 Pro machine, hosting Ubuntu 16.04 in Window Subsystem for Linux 1. The compiler reports this from gcc-9 -v:
Using built-in specs.
COLLECT_GCC=gcc-9
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/9/lto-wrapper
OFFLOAD_TARGET_NAMES=nvptx-none:hsa
OFFLOAD_TARGET_DEFAULT=1
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu 9.4.0-1ubuntu1~16.04' --with-bugurl=file:///usr/share/doc/gcc-9/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,gm2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-9 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --with-system-zlib --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-9-SATzbE/gcc-9-9.4.0/debian/tmp-nvptx/usr,hsa --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 9.4.0 (Ubuntu 9.4.0-1ubuntu1~16.04)
The command to compile was simple: I commented out the definition of FLOATTYPE so that I could compile different versions, and ran:
gcc-9 test_precision0100.c -o test_precision0100_longdouble.exe -DFLOATTYPE="long double"
Then ran ./test_precision0100_longdouble.exe.
The compiler does not give any warning messages with -Wall -Wextra -pedantic -std=c99, other than unused-parameter for argc and argv.
I get the same results with FLOATTYPE defined in the code as supplied above. I also get the same anomalous results with the built-in gcc v5.4.0, but not on another machine hosting Ubuntu 18.04 on WSL2. The output looks how you'd expect from my description, ending:
1.0<< 50: x=1125899906842624 ... can increment without loss of precision
1.0<< 51: x=2251799813685248 ... can increment without loss of precision
1.0<< 52: x=4503599627370496 ... can increment without loss of precision
1.0<< 53: x=9007199254740992 ... cannot increment without loss of precision
FLT_RADIX = 2
FLT_MANT_DIG = 24
DBL_MANT_DIG = 53
LDBL_MANT_DIG = 64
sizeof(FLOATTYPE) = 16
Here is "test_precision0100.s" from "gcc -S test_precision0100.c" as above:
.file "test_precision0100.c"
.text
.section .rodata
.LC1:
.string "1.0<<%3d: x=%.0Lf"
.align 8
.LC2:
.string " ... can increment without loss of precision"
.align 8
.LC3:
.string " ... cannot increment without loss of precision"
.LC4:
.string "FLT_RADIX = %d\n"
.LC5:
.string "FLT_MANT_DIG = %d\n"
.LC6:
.string "DBL_MANT_DIG = %d\n"
.LC7:
.string "LDBL_MANT_DIG = %d\n"
.LC8:
.string "\nsizeof(FLOATTYPE) = %lu\n"
.text
.globl main
.type main, @function
main:
.LFB2:
.cfi_startproc
pushq %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp
.cfi_def_cfa_register 6
subq $48, %rsp
movl %edi, -36(%rbp)
movq %rsi, -48(%rbp)
fld1
fstpt -16(%rbp)
movl $0, -20(%rbp)
.L5:
movl -20(%rbp), %eax
pushq -8(%rbp)
pushq -16(%rbp)
movl %eax, %esi
movl $.LC1, %edi
movl $0, %eax
call printf
addq $16, %rsp
fldt -16(%rbp)
fld1
faddp %st, %st(1)
fldt -16(%rbp)
fucomip %st(1), %st
jp .L9
fldt -16(%rbp)
fucomip %st(1), %st
fstp %st(0)
je .L2
jmp .L7
.L9:
fstp %st(0)
.L7:
fldt -16(%rbp)
fld1
faddp %st, %st(1)
fldt -16(%rbp)
fsubrp %st, %st(1)
fld1
fucomip %st(1), %st
jp .L10
fld1
fucomip %st(1), %st
fstp %st(0)
jne .L2
movl $.LC2, %edi
call puts
fldt -16(%rbp)
fadd %st(0), %st
fstpt -16(%rbp)
addl $1, -20(%rbp)
jmp .L5
.L10:
fstp %st(0)
.L2:
movl $.LC3, %edi
call puts
nop
movl $2, %esi
movl $.LC4, %edi
movl $0, %eax
call printf
movl $24, %esi
movl $.LC5, %edi
movl $0, %eax
call printf
movl $53, %esi
movl $.LC6, %edi
movl $0, %eax
call printf
movl $64, %esi
movl $.LC7, %edi
movl $0, %eax
call printf
movl $16, %esi
movl $.LC8, %edi
movl $0, %eax
call printf
movl $0, %eax
leave
.cfi_def_cfa 7, 8
ret
.cfi_endproc
.LFE2:
.size main, .-main
.ident "GCC: (Ubuntu 9.4.0-1ubuntu1~16.04) 9.4.0"
.section .note.GNU-stack,"",@progbits
That is a long-standing WSL1 bug - https://github.com/microsoft/WSL/issues/830
You can try to work around it by adding:
#include <fpu_control.h>
...
int main() {
unsigned short cw = 0x37f;
_FPU_SETCW(cw);
to your program.