Mon 24 Jul 2017 11:49:30 AM UTC, original submission:
According to FITS standard (v3), "The [BZERO] value field shall contain a floating-point number ...". Hence, in Gnuastro, we ask CFITSIO to read it as double precision floating point.
But a double precision floating point can only safely convert an integer to a float and back to an integer upto 15 significant decimal digits. Hence, when using BZERO to read the array as a differently signed integer of the same width, we will get accurate results for signed 8-bit, unsigned 16-bit and unsigned 32-bit integer datasets.
However, the integer `9223372036854775808' that must be used for BZERO to read the dataset as an unsigned 64-bit integer has 19 significant decimal digits. So it (or integers close to it) cannot be accurately resolved in `double' type (to compare with the standard value).
As one example, attached I am sending `bzero-u64.c' (a small C program to show the problem). This value is first read as `uint64_t', then as a `double' and then saved into a double from a string (with `strtod'). But in the latter two cases, it is subtracted by 1 and 2 respectively. So when this program is compiled and run, we see that in printing and comparing, the compiled program can't distinguish the subtractions. Hence `d' is equal to both `d_str' and the string constant.
I also tried this with CFITSIO. I made a 3x3 `uint64_t' array and saved it as a FITS file with CFITSIO. In the produced FITS file, I manually decreased the `BZERO' value by 1. The FITS file with the modified BZERO keyword is attached as `u64.fits'. When I use CFITSIO to read the value of BZERO as double precision floating point, this wrong value is interpretted as the correct value and like the example C program, the equality condition succeeds.
I got in touch with William Pence (author of CFITSIO) about this issue. The main point he raised was that we can use larger floating point containers (for example a 128-bit floating point), so the issue that this number does not fit into a `double' only an implementation detail. When the users don't have access to such types, the safest way would be to read the BZERO value as a string and compare it with the string '9223372036854775808' (not as a number).
I am busy implementing this now. I will read the BZERO keyword as a string by default and then parse it as a `double' with `strtod' in Gnuastro. This way, we can have both the string and the number with only one attempt at reading the keyword. When the type of th e input dataset is a 64-bit integer, we will compare the BZERO value with a string and in the other cases, we will compare with the number.
|