diff --git a/BITGROOM/example/CMakeLists.txt b/BITGROOM/example/CMakeLists.txt index 2062fe2ff..ad4297f6d 100755 --- a/BITGROOM/example/CMakeLists.txt +++ b/BITGROOM/example/CMakeLists.txt @@ -43,6 +43,15 @@ if (H5PL_BUILD_TESTING) target_link_libraries (h5repack_floats_bitgroom PRIVATE dl) endif () + # Round-trip test for special floating-point values (NaN, +/- Inf, +/- zero). + # See the test definition lower in this file for the full description. + add_executable (h5filter_specials_bitgroom ${PROJECT_SOURCE_DIR}/h5filter_specials_bitgroom.c) + TARGET_C_PROPERTIES (h5filter_specials_bitgroom ${LIB_TYPE}) + target_link_libraries (h5filter_specials_bitgroom PRIVATE ${H5PL_HDF5_LINK_LIBS}) + if (NOT WIN32) + target_link_libraries (h5filter_specials_bitgroom PRIVATE dl) + endif () + macro (ADD_H5_TEST testname) add_test ( NAME ${testname}-clearall @@ -261,6 +270,44 @@ if (H5PL_BUILD_TESTING) if (NOT DISABLE_H5BITGROOM_ENCODER) #UD BITGROOM ADD_H5_UD_TEST (ud_convert 0 h5repack_floats.h5 --enable-error-stack -v -f UD=32022,0,5,3,4,0,0,0 -l CHUNK=4x8) + + # Special-value round-trip test: writes a small chunked dataset containing + # NaN, +/- Inf, +/- zero, plus normal values; applies the BitGroom filter + # via H5Pset_filter; reads back; and asserts in C that the special values + # are preserved bit-exact while at least one normal value was quantized. + # BitGroom alternates shave (AND) on even indices and set (OR) on odd + # indices; the special values straddle both. The verdict comes from the + # program's exit code; no h5dump comparison is involved, so cross-platform + # NaN/Inf printf formatting cannot affect the result. TESTLIBDIR is + # already set by ADD_H5_UD_TEST above. + add_test ( + NAME H5BITGROOM_UD-specials-roundtrip-clearall + COMMAND ${CMAKE_COMMAND} -E remove + h5filter_specials_bitgroom.h5 + h5filter_specials_bitgroom.out + h5filter_specials_bitgroom.out.err + ) + if (NOT "${last_test}" STREQUAL "") + set_tests_properties (H5BITGROOM_UD-specials-roundtrip-clearall PROPERTIES DEPENDS ${last_test}) + endif () + set (last_test "H5BITGROOM_UD-specials-roundtrip-clearall") + add_test ( + NAME H5BITGROOM_UD-specials-roundtrip + COMMAND "${CMAKE_COMMAND}" + -D "TEST_PROGRAM=$" + -D "TEST_FOLDER=${PROJECT_BINARY_DIR}" + -D "TEST_EXPECT=0" + -D "TEST_OUTPUT=h5filter_specials_bitgroom.out" + -D "TEST_SKIP_COMPARE=1" + -D "TEST_LIBRARY_DIRECTORY=${TESTLIBDIR}" + -D "TEST_ENV_VAR=HDF5_PLUGIN_PATH" + -D "TEST_ENV_VALUE=${CMAKE_BINARY_DIR}/plugins" + -P "${H5BITGROOM_RESOURCES_DIR}/runTest.cmake" + ) + set_tests_properties (H5BITGROOM_UD-specials-roundtrip PROPERTIES + WORKING_DIRECTORY "${PROJECT_BINARY_DIR}" + DEPENDS H5BITGROOM_UD-specials-roundtrip-clearall) + set (last_test "H5BITGROOM_UD-specials-roundtrip") endif () endif () diff --git a/BITGROOM/example/h5filter_specials_bitgroom.c b/BITGROOM/example/h5filter_specials_bitgroom.c new file mode 100644 index 000000000..ca87770cb --- /dev/null +++ b/BITGROOM/example/h5filter_specials_bitgroom.c @@ -0,0 +1,228 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + * Copyright by The HDF Group. * + * All rights reserved. * + * * + * This file is part of the HDF5 BitGroom filter plugin source. * + * The full copyright notice, including terms governing use, modification, * + * and redistribution, is contained in the file COPYING, which can be found * + * at the root of the source code distribution tree. If you do not have * + * access to this file, you may request a copy from help@hdfgroup.org. * + * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ + +/* + * Round-trip test for special floating-point values through the BitGroom + * filter. Verifies that NaN, +Inf, -Inf, -0.0, and +0.0 are passed through + * the filter unchanged (bit-identical) on both the float and double code + * paths, and that at least one normal value is quantized (proving the + * filter actually ran). The filter is loaded dynamically via HDF5's plugin + * mechanism (HDF5_PLUGIN_PATH); no h5dump comparison is involved so cross- + * platform NaN/Inf printf formatting cannot affect the result. + * + * BitGroom alternates shave (AND-mask) on even indices and set (OR-mask) + * on odd indices; the special values are placed at both even and odd + * positions to exercise both loops. + */ + +#include "hdf5.h" +#include +#include +#include + +#define FILENAME "h5filter_specials_bitgroom.h5" +#define DSET_F32 "f32" +#define DSET_F64 "f64" +#define N 16 +#define N_SPECIALS 5 +#define H5Z_FILTER_BITGROOM 32022 + +/* Build exact bit patterns via memcpy so the compiler cannot constant-fold + * away -0.0f or normalize a quiet NaN to a different payload. */ +static void +init_input_f32(float *buf, uint32_t *bits_out) +{ + const uint32_t specials[N_SPECIALS] = { + 0x7FC00000U, /* quiet NaN */ + 0x7F800000U, /* +Inf */ + 0xFF800000U, /* -Inf */ + 0x80000000U, /* -0.0 */ + 0x00000000U, /* +0.0 */ + }; + for (int i = 0; i < N_SPECIALS; i++) + memcpy(&buf[i], &specials[i], sizeof(float)); + for (int i = N_SPECIALS; i < N; i++) + buf[i] = ((float)i - 3.7f) * 0.137f; + memcpy(bits_out, buf, N * sizeof(uint32_t)); +} + +static void +init_input_f64(double *buf, uint64_t *bits_out) +{ + const uint64_t specials[N_SPECIALS] = { + 0x7FF8000000000000ULL, /* quiet NaN */ + 0x7FF0000000000000ULL, /* +Inf */ + 0xFFF0000000000000ULL, /* -Inf */ + 0x8000000000000000ULL, /* -0.0 */ + 0x0000000000000000ULL, /* +0.0 */ + }; + for (int i = 0; i < N_SPECIALS; i++) + memcpy(&buf[i], &specials[i], sizeof(double)); + for (int i = N_SPECIALS; i < N; i++) + buf[i] = ((double)i - 3.7) * 0.137; + memcpy(bits_out, buf, N * sizeof(uint64_t)); +} + +static int +roundtrip_f32(hid_t file_id) +{ + hid_t space_id = H5I_INVALID_HID, dcpl_id = H5I_INVALID_HID, dset_id = H5I_INVALID_HID; + int ret = 1; + hsize_t dims[1] = {N}, chunk[1] = {N}; + /* NSD=3, sizeof(data)=4, has_mss_val=0, mss_val_byt_1to4=0, mss_val_byt_5to8=0 */ + const unsigned int cd_values[5] = {3, 4, 0, 0, 0}; + float input[N], output[N]; + uint32_t input_bits[N], output_bits[N]; + int any_quantized = 0; + + init_input_f32(input, input_bits); + + if ((space_id = H5Screate_simple(1, dims, NULL)) < 0) + goto done; + if ((dcpl_id = H5Pcreate(H5P_DATASET_CREATE)) < 0) + goto done; + if (H5Pset_chunk(dcpl_id, 1, chunk) < 0) + goto done; + if (H5Pset_filter(dcpl_id, H5Z_FILTER_BITGROOM, H5Z_FLAG_MANDATORY, 5, cd_values) < 0) + goto done; + if ((dset_id = + H5Dcreate(file_id, DSET_F32, H5T_IEEE_F32LE, space_id, H5P_DEFAULT, dcpl_id, H5P_DEFAULT)) < 0) + goto done; + if (H5Dwrite(dset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, input) < 0) + goto done; + H5Dclose(dset_id); + dset_id = H5I_INVALID_HID; + + if ((dset_id = H5Dopen(file_id, DSET_F32, H5P_DEFAULT)) < 0) + goto done; + if (H5Dread(dset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, output) < 0) + goto done; + memcpy(output_bits, output, sizeof(output_bits)); + + for (int i = 0; i < N_SPECIALS; i++) { + if (input_bits[i] != output_bits[i]) { + fprintf(stderr, "FAIL [f32, idx %d]: special value mutated: in=0x%08x out=0x%08x\n", i, + input_bits[i], output_bits[i]); + goto done; + } + } + for (int i = N_SPECIALS; i < N; i++) { + if (input_bits[i] != output_bits[i]) { + any_quantized = 1; + break; + } + } + if (!any_quantized) { + fprintf(stderr, "FAIL [f32]: no normal value was quantized; filter did not run\n"); + goto done; + } + + ret = 0; + +done: + if (dset_id >= 0) + H5Dclose(dset_id); + if (dcpl_id >= 0) + H5Pclose(dcpl_id); + if (space_id >= 0) + H5Sclose(space_id); + return ret; +} + +static int +roundtrip_f64(hid_t file_id) +{ + hid_t space_id = H5I_INVALID_HID, dcpl_id = H5I_INVALID_HID, dset_id = H5I_INVALID_HID; + int ret = 1; + hsize_t dims[1] = {N}, chunk[1] = {N}; + /* NSD=3, sizeof(data)=8, has_mss_val=0, mss_val_byt_1to4=0, mss_val_byt_5to8=0 */ + const unsigned int cd_values[5] = {3, 8, 0, 0, 0}; + double input[N], output[N]; + uint64_t input_bits[N], output_bits[N]; + int any_quantized = 0; + + init_input_f64(input, input_bits); + + if ((space_id = H5Screate_simple(1, dims, NULL)) < 0) + goto done; + if ((dcpl_id = H5Pcreate(H5P_DATASET_CREATE)) < 0) + goto done; + if (H5Pset_chunk(dcpl_id, 1, chunk) < 0) + goto done; + if (H5Pset_filter(dcpl_id, H5Z_FILTER_BITGROOM, H5Z_FLAG_MANDATORY, 5, cd_values) < 0) + goto done; + if ((dset_id = + H5Dcreate(file_id, DSET_F64, H5T_IEEE_F64LE, space_id, H5P_DEFAULT, dcpl_id, H5P_DEFAULT)) < 0) + goto done; + if (H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, input) < 0) + goto done; + H5Dclose(dset_id); + dset_id = H5I_INVALID_HID; + + if ((dset_id = H5Dopen(file_id, DSET_F64, H5P_DEFAULT)) < 0) + goto done; + if (H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, output) < 0) + goto done; + memcpy(output_bits, output, sizeof(output_bits)); + + for (int i = 0; i < N_SPECIALS; i++) { + if (input_bits[i] != output_bits[i]) { + fprintf(stderr, "FAIL [f64, idx %d]: special value mutated: in=0x%016llx out=0x%016llx\n", i, + (unsigned long long)input_bits[i], (unsigned long long)output_bits[i]); + goto done; + } + } + for (int i = N_SPECIALS; i < N; i++) { + if (input_bits[i] != output_bits[i]) { + any_quantized = 1; + break; + } + } + if (!any_quantized) { + fprintf(stderr, "FAIL [f64]: no normal value was quantized; filter did not run\n"); + goto done; + } + + ret = 0; + +done: + if (dset_id >= 0) + H5Dclose(dset_id); + if (dcpl_id >= 0) + H5Pclose(dcpl_id); + if (space_id >= 0) + H5Sclose(space_id); + return ret; +} + +int +main(void) +{ + hid_t file_id; + int r_f32 = 1, r_f64 = 1; + + file_id = H5Fcreate(FILENAME, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + if (file_id < 0) { + fprintf(stderr, "FAIL: H5Fcreate(%s) failed\n", FILENAME); + return 1; + } + + r_f32 = roundtrip_f32(file_id); + r_f64 = roundtrip_f64(file_id); + + H5Fclose(file_id); + + if (r_f32 == 0 && r_f64 == 0) { + printf("PASS: BitGroom special-value round-trip (f32 + f64)\n"); + return 0; + } + return 1; +} diff --git a/BITGROOM/src/H5Zbitgroom.c b/BITGROOM/src/H5Zbitgroom.c index 59859e9ca..a4f773c0f 100644 --- a/BITGROOM/src/H5Zbitgroom.c +++ b/BITGROOM/src/H5Zbitgroom.c @@ -78,7 +78,7 @@ 2 /**< Ordinal position of missing value flag in parameter list (cd_params array) */ #define CCR_FLT_PRM_PSN_MSS_VAL \ 3 /**< Ordinal position of missing value in parameter list (cd_params array) \ - NB: Missing value (_FillValue) uses two cd_params slots so it can be single or double-precision. \ + NB: HDF5 dataset fill value uses two cd_params slots so it can be single or double-precision. \ Single-precision values are read as first 4-bytes starting at cd_params[4] \ (and cd_params[5] is ignored), while double-precision values are read as first 8-bytes starting \ at cd_params[4] and ending with cd_params[5]. */ @@ -91,12 +91,6 @@ #ifndef NC_DOUBLE #define NC_DOUBLE 6 #endif /* !NC_DOUBLE */ -#ifndef NC_FILL_FLOAT -#define NC_FILL_FLOAT (9.9692099683868690e+36f) /* near 15 * 2^119 */ -#endif /* !NC_FILL_FLOAT */ -#ifndef NC_FILL_DOUBLE -#define NC_FILL_DOUBLE (9.9692099683868690e+36) -#endif /* !NC_FILL_DOUBLE */ /* Minimum number of explicit significant bits to preserve when zeroing/bit-masking floating point values Codes will preserve at least two explicit bits, IEEE significant representation contains one implicit bit @@ -485,8 +479,10 @@ ccr_bgr(const int nsd, const int type, const size_t sz, const int has_mss_val, p double prc_bnr_xct; /* [nbr] Binary digits of precision, exact */ double mss_val_cmp_dbl; /* Missing value for comparison to double precision values */ + double val_dbl; /* [frc] Copy of input value to avoid indirection */ float mss_val_cmp_flt; /* Missing value for comparison to single precision values */ + float val_flt; /* [frc] Copy of input value to avoid indirection */ int bit_xpl_nbr_sgn = -1; /* [nbr] Number of explicit bits in significand */ int bit_xpl_nbr_zro; /* [nbr] Number of explicit bits to zero */ @@ -527,11 +523,13 @@ ccr_bgr(const int nsd, const int type, const size_t sz, const int has_mss_val, p switch (type) { case NC_FLOAT: - /* Missing value for comparison is _FillValue (if any) otherwise default NC_FILL_FLOAT/DOUBLE */ + /* Comparison value is the dataset's HDF5 fill value (when user-defined), + * otherwise NaN -- which compares unequal to every finite value, so the + * fill-value test is an inert no-op when none was set on the dataset. */ if (has_mss_val) mss_val_cmp_flt = *mss_val.fp; else - mss_val_cmp_flt = NC_FILL_FLOAT; + mss_val_cmp_flt = NAN; bit_xpl_nbr_sgn = bit_xpl_nbr_sgn_flt; bit_xpl_nbr_zro = bit_xpl_nbr_sgn - prc_bnr_xpl_rqr; if (bit_xpl_nbr_zro > bit_xpl_nbr_sgn - NCO_PPC_BIT_XPL_NBR_MIN) @@ -549,21 +547,27 @@ ccr_bgr(const int nsd, const int type, const size_t sz, const int has_mss_val, p msk_f32_u32_one = ~msk_f32_u32_zro; // msk_f32_u32_hshv=msk_f32_u32_one & (msk_f32_u32_zro >> 1); /* Set one bit: the MSB of LSBs */ - /* Bit-Groom: alternately shave and set LSBs */ + /* Bit-Groom: alternately shave and set LSBs. + * Do not quantize the dataset fill value, +/- zero, NaN, or +/- Inf. + * isfinite(val_flt) covers !isnan(val_flt) and additionally blocks + * +/- Inf: the shave (AND) loop preserves Inf bit patterns, but the + * set (OR) loop turns +/- Inf (exponent all 1, mantissa 0) into a + * NaN (exponent all 1, mantissa non-zero) by ORing low mantissa bits. */ for (idx = 0L; idx < sz; idx += 2L) - if (op1.fp[idx] != mss_val_cmp_flt) + if ((val_flt = op1.fp[idx]) != mss_val_cmp_flt && val_flt != 0.0f && isfinite(val_flt)) u32_ptr[idx] &= msk_f32_u32_zro; for (idx = 1L; idx < sz; idx += 2L) - if (op1.fp[idx] != mss_val_cmp_flt && - u32_ptr[idx] != 0U) /* Never quantize upwards floating point values of zero */ + if ((val_flt = op1.fp[idx]) != mss_val_cmp_flt && val_flt != 0.0f && isfinite(val_flt)) u32_ptr[idx] |= msk_f32_u32_one; break; case NC_DOUBLE: - /* Missing value for comparison is _FillValue (if any) otherwise default NC_FILL_FLOAT/DOUBLE */ + /* Comparison value is the dataset's HDF5 fill value (when user-defined), + * otherwise NaN -- which compares unequal to every finite value, so the + * fill-value test is an inert no-op when none was set on the dataset. */ if (has_mss_val) mss_val_cmp_dbl = *mss_val.dp; else - mss_val_cmp_dbl = NC_FILL_DOUBLE; + mss_val_cmp_dbl = NAN; bit_xpl_nbr_sgn = bit_xpl_nbr_sgn_dbl; bit_xpl_nbr_zro = bit_xpl_nbr_sgn - prc_bnr_xpl_rqr; if (bit_xpl_nbr_zro > bit_xpl_nbr_sgn - NCO_PPC_BIT_XPL_NBR_MIN) @@ -580,13 +584,12 @@ ccr_bgr(const int nsd, const int type, const size_t sz, const int has_mss_val, p /* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */ msk_f64_u64_one = ~msk_f64_u64_zro; // msk_f64_u64_hshv=msk_f64_u64_one & (msk_f64_u64_zro >> 1); /* Set one bit: the MSB of LSBs */ - /* Bit-Groom: alternately shave and set LSBs */ + /* See float branch above for why isfinite(val_dbl) is used instead of !isnan(val_dbl). */ for (idx = 0L; idx < sz; idx += 2L) - if (op1.dp[idx] != mss_val_cmp_dbl) + if ((val_dbl = op1.dp[idx]) != mss_val_cmp_dbl && val_dbl != 0.0 && isfinite(val_dbl)) u64_ptr[idx] &= msk_f64_u64_zro; for (idx = 1L; idx < sz; idx += 2L) - if (op1.dp[idx] != mss_val_cmp_dbl && - u64_ptr[idx] != 0ULL) /* Never quantize upwards floating point values of zero */ + if ((val_dbl = op1.dp[idx]) != mss_val_cmp_dbl && val_dbl != 0.0 && isfinite(val_dbl)) u64_ptr[idx] |= msk_f64_u64_one; break; default: diff --git a/BITROUND/example/CMakeLists.txt b/BITROUND/example/CMakeLists.txt index 9be01cc55..c4e19c07f 100644 --- a/BITROUND/example/CMakeLists.txt +++ b/BITROUND/example/CMakeLists.txt @@ -43,6 +43,15 @@ if (H5PL_BUILD_TESTING) target_link_libraries (h5repack_floats_bitround PRIVATE dl) endif () + # Round-trip test for special floating-point values (NaN, +/- Inf, +/- zero). + # See the test definition lower in this file for the full description. + add_executable (h5filter_specials_bitround ${PROJECT_SOURCE_DIR}/h5filter_specials_bitround.c) + TARGET_C_PROPERTIES (h5filter_specials_bitround ${LIB_TYPE}) + target_link_libraries (h5filter_specials_bitround PRIVATE ${H5PL_HDF5_LINK_LIBS}) + if (NOT WIN32) + target_link_libraries (h5filter_specials_bitround PRIVATE dl) + endif () + macro (ADD_H5_TEST testname) add_test ( NAME ${testname}-clearall @@ -266,6 +275,43 @@ if (H5PL_BUILD_TESTING) if (NOT DISABLE_H5BITROUND_ENCODER) #UD BITROUND ADD_H5_UD_TEST (ud_convert 0 h5repack_floats.h5 --enable-error-stack -v -f UD=32023,0,5,3,4,0,0,0 -l CHUNK=4x8) + + # Special-value round-trip test: writes a small chunked dataset containing + # NaN, +/- Inf, +/- zero, plus normal values; applies the Granular BitRound + # filter via H5Pset_filter; reads back; and asserts in C that the special + # values are preserved bit-exact while at least one normal value was + # quantized. The verdict comes from the program's exit code; no h5dump + # comparison is involved, so cross-platform NaN/Inf printf formatting + # cannot affect the result. TESTLIBDIR is already set by ADD_H5_UD_TEST + # above. + add_test ( + NAME H5BITROUND_UD-specials-roundtrip-clearall + COMMAND ${CMAKE_COMMAND} -E remove + h5filter_specials_bitround.h5 + h5filter_specials_bitround.out + h5filter_specials_bitround.out.err + ) + if (NOT "${last_test}" STREQUAL "") + set_tests_properties (H5BITROUND_UD-specials-roundtrip-clearall PROPERTIES DEPENDS ${last_test}) + endif () + set (last_test "H5BITROUND_UD-specials-roundtrip-clearall") + add_test ( + NAME H5BITROUND_UD-specials-roundtrip + COMMAND "${CMAKE_COMMAND}" + -D "TEST_PROGRAM=$" + -D "TEST_FOLDER=${PROJECT_BINARY_DIR}" + -D "TEST_EXPECT=0" + -D "TEST_OUTPUT=h5filter_specials_bitround.out" + -D "TEST_SKIP_COMPARE=1" + -D "TEST_LIBRARY_DIRECTORY=${TESTLIBDIR}" + -D "TEST_ENV_VAR=HDF5_PLUGIN_PATH" + -D "TEST_ENV_VALUE=${CMAKE_BINARY_DIR}/plugins" + -P "${H5BITROUND_RESOURCES_DIR}/runTest.cmake" + ) + set_tests_properties (H5BITROUND_UD-specials-roundtrip PROPERTIES + WORKING_DIRECTORY "${PROJECT_BINARY_DIR}" + DEPENDS H5BITROUND_UD-specials-roundtrip-clearall) + set (last_test "H5BITROUND_UD-specials-roundtrip") endif () endif () diff --git a/BITROUND/example/h5filter_specials_bitround.c b/BITROUND/example/h5filter_specials_bitround.c new file mode 100644 index 000000000..2e0d8384e --- /dev/null +++ b/BITROUND/example/h5filter_specials_bitround.c @@ -0,0 +1,224 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + * Copyright by The HDF Group. * + * All rights reserved. * + * * + * This file is part of the HDF5 Granular BitRound filter plugin source. * + * The full copyright notice, including terms governing use, modification, * + * and redistribution, is contained in the file COPYING, which can be found * + * at the root of the source code distribution tree. If you do not have * + * access to this file, you may request a copy from help@hdfgroup.org. * + * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ + +/* + * Round-trip test for special floating-point values through the Granular + * BitRound filter. Verifies that NaN, +Inf, -Inf, -0.0, and +0.0 are passed + * through the filter unchanged (bit-identical) on both the float and double + * code paths, and that at least one normal value is quantized (proving the + * filter actually ran). The filter is loaded dynamically via HDF5's plugin + * mechanism (HDF5_PLUGIN_PATH); no h5dump comparison is involved so cross- + * platform NaN/Inf printf formatting cannot affect the result. + */ + +#include "hdf5.h" +#include +#include +#include + +#define FILENAME "h5filter_specials_bitround.h5" +#define DSET_F32 "f32" +#define DSET_F64 "f64" +#define N 16 +#define N_SPECIALS 5 +#define H5Z_FILTER_GRANULARBR 32023 + +/* Build exact bit patterns via memcpy so the compiler cannot constant-fold + * away -0.0f or normalize a quiet NaN to a different payload. */ +static void +init_input_f32(float *buf, uint32_t *bits_out) +{ + const uint32_t specials[N_SPECIALS] = { + 0x7FC00000U, /* quiet NaN */ + 0x7F800000U, /* +Inf */ + 0xFF800000U, /* -Inf */ + 0x80000000U, /* -0.0 */ + 0x00000000U, /* +0.0 */ + }; + for (int i = 0; i < N_SPECIALS; i++) + memcpy(&buf[i], &specials[i], sizeof(float)); + for (int i = N_SPECIALS; i < N; i++) + buf[i] = ((float)i - 3.7f) * 0.137f; + memcpy(bits_out, buf, N * sizeof(uint32_t)); +} + +static void +init_input_f64(double *buf, uint64_t *bits_out) +{ + const uint64_t specials[N_SPECIALS] = { + 0x7FF8000000000000ULL, /* quiet NaN */ + 0x7FF0000000000000ULL, /* +Inf */ + 0xFFF0000000000000ULL, /* -Inf */ + 0x8000000000000000ULL, /* -0.0 */ + 0x0000000000000000ULL, /* +0.0 */ + }; + for (int i = 0; i < N_SPECIALS; i++) + memcpy(&buf[i], &specials[i], sizeof(double)); + for (int i = N_SPECIALS; i < N; i++) + buf[i] = ((double)i - 3.7) * 0.137; + memcpy(bits_out, buf, N * sizeof(uint64_t)); +} + +static int +roundtrip_f32(hid_t file_id) +{ + hid_t space_id = H5I_INVALID_HID, dcpl_id = H5I_INVALID_HID, dset_id = H5I_INVALID_HID; + int ret = 1; + hsize_t dims[1] = {N}, chunk[1] = {N}; + /* NSD=3, sizeof(data)=4, has_mss_val=0, mss_val_byt_1to4=0, mss_val_byt_5to8=0 */ + const unsigned int cd_values[5] = {3, 4, 0, 0, 0}; + float input[N], output[N]; + uint32_t input_bits[N], output_bits[N]; + int any_quantized = 0; + + init_input_f32(input, input_bits); + + if ((space_id = H5Screate_simple(1, dims, NULL)) < 0) + goto done; + if ((dcpl_id = H5Pcreate(H5P_DATASET_CREATE)) < 0) + goto done; + if (H5Pset_chunk(dcpl_id, 1, chunk) < 0) + goto done; + if (H5Pset_filter(dcpl_id, H5Z_FILTER_GRANULARBR, H5Z_FLAG_MANDATORY, 5, cd_values) < 0) + goto done; + if ((dset_id = + H5Dcreate(file_id, DSET_F32, H5T_IEEE_F32LE, space_id, H5P_DEFAULT, dcpl_id, H5P_DEFAULT)) < 0) + goto done; + if (H5Dwrite(dset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, input) < 0) + goto done; + H5Dclose(dset_id); + dset_id = H5I_INVALID_HID; + + if ((dset_id = H5Dopen(file_id, DSET_F32, H5P_DEFAULT)) < 0) + goto done; + if (H5Dread(dset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, output) < 0) + goto done; + memcpy(output_bits, output, sizeof(output_bits)); + + for (int i = 0; i < N_SPECIALS; i++) { + if (input_bits[i] != output_bits[i]) { + fprintf(stderr, "FAIL [f32, idx %d]: special value mutated: in=0x%08x out=0x%08x\n", i, + input_bits[i], output_bits[i]); + goto done; + } + } + for (int i = N_SPECIALS; i < N; i++) { + if (input_bits[i] != output_bits[i]) { + any_quantized = 1; + break; + } + } + if (!any_quantized) { + fprintf(stderr, "FAIL [f32]: no normal value was quantized; filter did not run\n"); + goto done; + } + + ret = 0; + +done: + if (dset_id >= 0) + H5Dclose(dset_id); + if (dcpl_id >= 0) + H5Pclose(dcpl_id); + if (space_id >= 0) + H5Sclose(space_id); + return ret; +} + +static int +roundtrip_f64(hid_t file_id) +{ + hid_t space_id = H5I_INVALID_HID, dcpl_id = H5I_INVALID_HID, dset_id = H5I_INVALID_HID; + int ret = 1; + hsize_t dims[1] = {N}, chunk[1] = {N}; + /* NSD=3, sizeof(data)=8, has_mss_val=0, mss_val_byt_1to4=0, mss_val_byt_5to8=0 */ + const unsigned int cd_values[5] = {3, 8, 0, 0, 0}; + double input[N], output[N]; + uint64_t input_bits[N], output_bits[N]; + int any_quantized = 0; + + init_input_f64(input, input_bits); + + if ((space_id = H5Screate_simple(1, dims, NULL)) < 0) + goto done; + if ((dcpl_id = H5Pcreate(H5P_DATASET_CREATE)) < 0) + goto done; + if (H5Pset_chunk(dcpl_id, 1, chunk) < 0) + goto done; + if (H5Pset_filter(dcpl_id, H5Z_FILTER_GRANULARBR, H5Z_FLAG_MANDATORY, 5, cd_values) < 0) + goto done; + if ((dset_id = + H5Dcreate(file_id, DSET_F64, H5T_IEEE_F64LE, space_id, H5P_DEFAULT, dcpl_id, H5P_DEFAULT)) < 0) + goto done; + if (H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, input) < 0) + goto done; + H5Dclose(dset_id); + dset_id = H5I_INVALID_HID; + + if ((dset_id = H5Dopen(file_id, DSET_F64, H5P_DEFAULT)) < 0) + goto done; + if (H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, output) < 0) + goto done; + memcpy(output_bits, output, sizeof(output_bits)); + + for (int i = 0; i < N_SPECIALS; i++) { + if (input_bits[i] != output_bits[i]) { + fprintf(stderr, "FAIL [f64, idx %d]: special value mutated: in=0x%016llx out=0x%016llx\n", i, + (unsigned long long)input_bits[i], (unsigned long long)output_bits[i]); + goto done; + } + } + for (int i = N_SPECIALS; i < N; i++) { + if (input_bits[i] != output_bits[i]) { + any_quantized = 1; + break; + } + } + if (!any_quantized) { + fprintf(stderr, "FAIL [f64]: no normal value was quantized; filter did not run\n"); + goto done; + } + + ret = 0; + +done: + if (dset_id >= 0) + H5Dclose(dset_id); + if (dcpl_id >= 0) + H5Pclose(dcpl_id); + if (space_id >= 0) + H5Sclose(space_id); + return ret; +} + +int +main(void) +{ + hid_t file_id; + int r_f32 = 1, r_f64 = 1; + + file_id = H5Fcreate(FILENAME, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + if (file_id < 0) { + fprintf(stderr, "FAIL: H5Fcreate(%s) failed\n", FILENAME); + return 1; + } + + r_f32 = roundtrip_f32(file_id); + r_f64 = roundtrip_f64(file_id); + + H5Fclose(file_id); + + if (r_f32 == 0 && r_f64 == 0) { + printf("PASS: BitRound special-value round-trip (f32 + f64)\n"); + return 0; + } + return 1; +} diff --git a/BITROUND/src/H5Zgranularbr.c b/BITROUND/src/H5Zgranularbr.c index 8e468da81..fa950760d 100644 --- a/BITROUND/src/H5Zgranularbr.c +++ b/BITROUND/src/H5Zgranularbr.c @@ -75,7 +75,7 @@ 2 /**< Ordinal position of missing value flag in parameter list (cd_params array) */ #define CCR_FLT_PRM_PSN_MSS_VAL \ 3 /**< Ordinal position of missing value in parameter list (cd_params array) \ - NB: Missing value (_FillValue) uses two cd_params slots so it can be single or double-precision. \ + NB: HDF5 dataset fill value uses two cd_params slots so it can be single or double-precision. \ Single-precision values are read as first 4-bytes starting at cd_params[4] \ (and cd_params[5] is ignored), while double-precision values are read as first 8-bytes starting \ at cd_params[4] and ending with cd_params[5]. */ @@ -88,12 +88,6 @@ #ifndef NC_DOUBLE #define NC_DOUBLE 6 #endif /* !NC_DOUBLE */ -#ifndef NC_FILL_FLOAT -#define NC_FILL_FLOAT (9.9692099683868690e+36f) /* near 15 * 2^119 */ -#endif /* !NC_FILL_FLOAT */ -#ifndef NC_FILL_DOUBLE -#define NC_FILL_DOUBLE (9.9692099683868690e+36) -#endif /* !NC_FILL_DOUBLE */ /* Minimum number of explicit significant bits to preserve when zeroing/bit-masking floating point values Codes will preserve at least two explicit bits, IEEE significant representation contains one implicit bit @@ -636,15 +630,23 @@ ccr_gbr(const int nsd, const int type, const size_t sz, const int has_mss_val, p switch (type) { case NC_FLOAT: - /* Missing value for comparison is _FillValue (if any) otherwise default NC_FILL_FLOAT/DOUBLE */ + /* Comparison value is the dataset's HDF5 fill value (when user-defined), + * otherwise NaN -- which compares unequal to every finite value, so the + * fill-value test is an inert no-op when none was set on the dataset. */ if (has_mss_val) mss_val_cmp_flt = *mss_val.fp; else - mss_val_cmp_flt = NC_FILL_FLOAT; + mss_val_cmp_flt = NAN; bit_xpl_nbr_sgn = bit_xpl_nbr_sgn_flt; u32_ptr = op1.ui32p; + /* Do not quantize the dataset fill value, +/- zero, NaN, or +/- Inf. + * isfinite(val) covers !isnan(val) and additionally blocks +/- Inf, + * which would otherwise reach the bit-shift below with an out-of-range + * shift count (undefined behavior). netcdf-c upstream does not gate Inf; + * we diverge here because the standalone filter's loop body computes + * per-value bit masks that Inf inputs corrupt. */ for (idx = 0L; idx < sz; idx++) { - if ((val = op1.fp[idx]) != mss_val_cmp_flt && u32_ptr[idx] != 0U) { + if ((val = op1.fp[idx]) != mss_val_cmp_flt && val != 0.0 && isfinite(val)) { mnt = frexp(val, &xpn_bs2); /* DGG19 p. 4102 (8) */ mnt_fabs = fabs(mnt); /* dgt_nbr = floor(log10(|val|)) + 1, computed via deterministic @@ -657,10 +659,8 @@ ccr_gbr(const int nsd, const int type, const size_t sz, const int has_mss_val, p * exactly 0.5 (val is an exact power of 2). This integer-only * formulation is bit-deterministic across platforms; the * original libm-based form is not. */ - prc_bnr_xpl_rqr = (mnt_fabs == 0.0) - ? 0 - : (unsigned short)abs(((mnt_fabs == 0.5) ? xpn_bs2 + 1 : xpn_bs2) - - qnt_pwr); /* Protect against mnt = -0.0 */ + prc_bnr_xpl_rqr = + (unsigned short)abs(((mnt_fabs == 0.5) ? xpn_bs2 + 1 : xpn_bs2) - qnt_pwr); prc_bnr_xpl_rqr--; /* 20211003 Reduce formula result by 1 bit: Passes all tests, improves CR by ~10% */ @@ -681,16 +681,19 @@ ccr_gbr(const int nsd, const int type, const size_t sz, const int has_mss_val, p } break; case NC_DOUBLE: - /* Missing value for comparison is _FillValue (if any) otherwise default NC_FILL_FLOAT/DOUBLE */ + /* Comparison value is the dataset's HDF5 fill value (when user-defined), + * otherwise NaN -- which compares unequal to every finite value, so the + * fill-value test is an inert no-op when none was set on the dataset. */ if (has_mss_val) mss_val_cmp_dbl = *mss_val.dp; else - mss_val_cmp_dbl = NC_FILL_DOUBLE; + mss_val_cmp_dbl = NAN; bit_xpl_nbr_sgn = bit_xpl_nbr_sgn_dbl; u64_ptr = op1.ui64p; + /* See float branch above for why isfinite(val) is used instead of !isnan(val). */ for (idx = 0L; idx < sz; idx++) { - if ((val = op1.dp[idx]) != mss_val_cmp_dbl && u64_ptr[idx] != 0U) { + if ((val = op1.dp[idx]) != mss_val_cmp_dbl && val != 0.0 && isfinite(val)) { mnt = frexp(val, &xpn_bs2); /* DGG19 p. 4102 (8) */ mnt_fabs = fabs(mnt); /* dgt_nbr = floor(log10(|val|)) + 1, computed via deterministic @@ -703,10 +706,8 @@ ccr_gbr(const int nsd, const int type, const size_t sz, const int has_mss_val, p * exactly 0.5 (val is an exact power of 2). This integer-only * formulation is bit-deterministic across platforms; the * original libm-based form is not. */ - prc_bnr_xpl_rqr = (mnt_fabs == 0.0) - ? 0 - : (unsigned short)abs(((mnt_fabs == 0.5) ? xpn_bs2 + 1 : xpn_bs2) - - qnt_pwr); /* Protect against mnt = -0.0 */ + prc_bnr_xpl_rqr = + (unsigned short)abs(((mnt_fabs == 0.5) ? xpn_bs2 + 1 : xpn_bs2) - qnt_pwr); prc_bnr_xpl_rqr--; /* 20211003 Reduce formula result by 1 bit: Passes all tests, improves CR by ~10% */