/* * Two possible C type double division bugs when compiling with gcc -Ofast * * Using: MinGW GCC compiler: gcc (MinGW.org GCC-8.2.0-5) 8.2.0 * Microsoft Windows 7 Professional 64 bit operating system * Lenovo Thinkpad T420 laptop 64 bit computer * * Also tried on: GNU GCC compiler: gcc (Ubuntu 7.2.0-8ubuntu3.2) 7.2.0 * Lubuntu 32 bit operating system * Samsung N220 netbook 32 bit computer * * Compiling with gcc: -O0, -O1, -O2, O3, -Os, -Og; this code works as intended. * * But compiling with gcc -Ofast, it sometimes calculates the wrong "a" or "x". * This behaviour seems consistent and reproducible on both the above setups, * although the wrong behaviour is slightly different on each setup. * * I know essentially nothing about what the "-Ofast" option does, * and don't know whether using the "-Ofast" option might be expected to change * the way type double division sometimes works, or whether this is something * that merits further investigation. * * At the end of this file are about 60 lines of the actual code and * C preprocessor directives. The rest of the 300 lines are information * for the two possible bugs, with explanations of what the code * is doing, and examples of the program output. * * I hope the following is reasonably clear, but I am happy to answer questions. * * * CJDN is Chronological Julian Day Number for dates. * * * In "qx(int jd)": double x = floor((jd - 1867216.25) / 36524.25); * (jd = 1867217 is CJDN of Gregorian date 0400-03-01; "x" is 0.0;) * jd = 1830692 is CJDN of Gregorian date 0300-03-01; "x" is -1.0; * jd = 1976789 = 1830692 + 146097 is CJDN of 0700-03-01; "x" should be 3.0, * * but on the Samsung "x" is 2.0 when compiled using "-Ofast"; * this is the case whether compiling with "#define xx" or "#define xxz"; * * but on the Lenovo compiled using "-Ofast": * * when compiled with "#define xx": "x" is correctly 3.0; * has: printf("# c_jd_to_civil_x: jd %d; x %4.1f; z %d;\n", jd, x, 0); * output: * // qx:xx: jd 1830692; x -1.0; z ; * // qx:xx: jd 1976789; x 3.0; z ; * ^^^^^^ correct * * when compiled with "#define xxz": "x" is wrongly 2.0; * has: printf("# c_jd_to_civil_x: jd %d; x %4.1f; z ;\n", jd, x); * output: * // qx:xxz: jd 1830692; x -1.0; z 0; * // qx:xxz: jd 1976789; x 2.0; z 0; * ^^^^^^ wrong * I am surprised this small change in the program affects * how "gcc -Ofast" compiles the calculation of "x". * * * In "qax(int yy)", having "#define" one or more of: * "ax_4_11", "ax_4_2", "ax_3_1", "ax_4_3", "ax_3_2"; * * "qax(int yy)" is a cut-down version of two date calculation functions. * (These potential problems were first noticed when comparing versions * which use doubles with faster possible replacements * which use only integer variables and integer arithmetic.) * * "x": when yy >= 7 and yy % 4 == 3: * On the Samsung netbook and the Lenovo laptop: * double x = floor((jd - 1867216.25) / 36524.25); * "x" should be yy - 4 but is yy - 5; * this is when yy = 7 + 4 * k, for integer k >= 0, * that is the CJDN for Gregorian date y-03-01, with y = 0700 + 400 * k, * that is when jd = 1976789 + 146097 * k; * * "a": when yy <= -1: * int y = yy * 100; * double a = floor(y / 100.0); * "a" should be yy * but on the Samsung netbook is yy - 1. * The Lenovo laptop results are similar, but slightly different; * some of the wrong results for "a" depend on the range of "yy" being used * and/or whether there is one or more loops of "yy": * * Seven runs with "gcc -Ofast" on the Lenovo setup showing clear patterns. * For both "a" and "x": * (i) first value is the doubles calculation, which is sometimes wrong; * (ii) second value is the integers calculation, which is correct; * (iii) the value after the "?" is (i) - (ii), which should be zero, * but is non-zero when (i) is wrong. * * only "#define ax_4_11": yy = -4 to 11; negative yy: least is ok, others not; * // yy -4 y -400: a -4.0 -4 ? 0 # jd 1575023: x -8.0 -8 ? 0 * // yy -3 y -300: a -4.0 -3 ? -1 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -3.0 -2 ? -1 # jd 1648071: x -6.0 -6 ? 0 * // yy -1 y -100: a -2.0 -1 ? -1 # jd 1684595: x -5.0 -5 ? 0 * // yy 0 y 0: a 0.0 0 ? 0 # jd 1721120: x -4.0 -4 ? 0 * // yy 1 y 100: a 1.0 1 ? 0 # jd 1757644: x -3.0 -3 ? 0 * // yy 2 y 200: a 2.0 2 ? 0 # jd 1794168: x -2.0 -2 ? 0 * // yy 3 y 300: a 3.0 3 ? 0 # jd 1830692: x -1.0 -1 ? 0 * // yy 4 y 400: a 4.0 4 ? 0 # jd 1867217: x 0.0 0 ? 0 * // yy 5 y 500: a 5.0 5 ? 0 # jd 1903741: x 1.0 1 ? 0 * // yy 6 y 600: a 6.0 6 ? 0 # jd 1940265: x 2.0 2 ? 0 * // yy 7 y 700: a 7.0 7 ? 0 # jd 1976789: x 2.0 3 ? -1 * // yy 8 y 800: a 8.0 8 ? 0 # jd 2013314: x 4.0 4 ? 0 * // yy 9 y 900: a 9.0 9 ? 0 # jd 2049838: x 5.0 5 ? 0 * // yy 10 y 1000: a 10.0 10 ? 0 # jd 2086362: x 6.0 6 ? 0 * // yy 11 y 1100: a 11.0 11 ? 0 # jd 2122886: x 6.0 7 ? -1 * * only "#define ax_4_2": yy = -4 to -2; negative yy: least is ok, others not; * // yy -4 y -400: a -4.0 -4 ? 0 # jd 1575023: x -8.0 -8 ? 0 * // yy -3 y -300: a -4.0 -3 ? -1 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -3.0 -2 ? -1 # jd 1648071: x -6.0 -6 ? 0 * * only "#define ax_3_1": yy = -3 to -1; negative yy: least is ok, others not; * // yy -3 y -300: a -3.0 -3 ? 0 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -3.0 -2 ? -1 # jd 1648071: x -6.0 -6 ? 0 * // yy -1 y -100: a -2.0 -1 ? -1 # jd 1684595: x -5.0 -5 ? 0 * * only "#define ax_4_3": yy = -4 to -3; negative yy: both are ok; * // yy -4 y -400: a -4.0 -4 ? 0 # jd 1575023: x -8.0 -8 ? 0 * // yy -3 y -300: a -3.0 -3 ? 0 # jd 1611547: x -7.0 -7 ? 0 * * only "#define ax_3_2": yy = -3 to -2; negative yy: both are ok; * // yy -3 y -300: a -3.0 -3 ? 0 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -2.0 -2 ? 0 # jd 1648071: x -6.0 -6 ? 0 * * only "#define ax_4_2" and "#define ax_4_3": yy = -4 to -2 and yy = -4 to -4; * negative yy: all are wrong; * // * // yy -4 y -400: a -5.0 -4 ? -1 # jd 1575023: x -8.0 -8 ? 0 * // yy -3 y -300: a -4.0 -3 ? -1 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -3.0 -2 ? -1 # jd 1648071: x -6.0 -6 ? 0 * // * // yy -4 y -400: a -5.0 -4 ? -1 # jd 1575023: x -8.0 -8 ? 0 * // yy -3 y -300: a -4.0 -3 ? -1 # jd 1611547: x -7.0 -7 ? 0 * * only "#define ax_3_1" and "#define ax_3_2": yy = -3 to -1 and yy = -3 to -2; * negative yy: all are wrong; * // * // yy -3 y -300: a -4.0 -3 ? -1 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -3.0 -2 ? -1 # jd 1648071: x -6.0 -6 ? 0 * // yy -1 y -100: a -2.0 -1 ? -1 # jd 1684595: x -5.0 -5 ? 0 * // * // yy -3 y -300: a -4.0 -3 ? -1 # jd 1611547: x -7.0 -7 ? 0 * // yy -2 y -200: a -3.0 -2 ? -1 # jd 1648071: x -6.0 -6 ? 0 * * * * *** "Bugs" - Use this guidance: http://www.mingw.org/Reporting_Bugs * *** https://www.chiark.greenend.org.uk/~sgtatham/bugs.html * ***** Try it on anoher machine - good idea - done - see above * * The following version information is for the Lenovo setup. * * MinGW Version * The MinGW Version is the version of the mingw-runtime or w32api that is installed. * You can look at the mingw/include/_mingw.h file for this version number. * * ** In C:\MinGW\include\_mingw.h is the version number one of these? * $Id: _mingw.h.in,v 6f940ffde558 2019/01/22 11:36:30 keith $ * __MINGW32_VERSION = 1,000,000 * major + 1,000 * minor + patch * #define __MINGW32_VERSION 5002002L * #define __MINGW32_MAJOR_VERSION 5 * #define __MINGW32_MINOR_VERSION 2 * #define __MINGW32_PATCHLEVEL 2 * * ** In C:\MinGW\include\w32api.h is the version number one of the following? * $Id: w32api.h.in,v 0bd67cc9bc86 2017/03/20 20:01:38 keithmarshall $ * #define __W32API_VERSION 5002002L * #define __W32API_MAJOR_VERSION 5 * #define __W32API_MINOR_VERSION 2 * #define __W32API_PATCHLEVEL 2 * * Build Environment * Your environment will probably not be the same as our defaults. * Please tell us whether you are using a Microsoft shell, (cmd.exe or command.com), * an MSYS shell, a Cygwin shell, or if you are running a cross-compiler environment. * If you are using an MSYS shell, please also include the version number, * as determined by running the command: uname -a * * ** On Lenovo it's a Microsoft shell, not sure whether cmd.exe or command.com. * ** On Samsung it's a Terminal * * Output on Lenovo: * * C:\MinGW\bin>gcc -v * Using built-in specs. * COLLECT_GCC=gcc * COLLECT_LTO_WRAPPER=c:/mingw/bin/../libexec/gcc/mingw32/8.2.0/lto-wrapper.exe * Target: mingw32 * Configured with: ../src/gcc-8.2.0/configure --build=x86_64-pc-linux-gnu * --host=mingw32 --target=mingw32 --prefix=/mingw --disable-win32-registry * --with-arch=i586 --with-tune=generic * --enable-languages=c,c++,objc,obj-c++,fortran,ada * --with-pkgversion='MinGW.org GCC-8.2.0-5' --with-gmp=/mingw * --with-mpfr=/mingw --with-mpc=/mingw --enable-static --enable-shared * --enable-threads --with-dwarf2 --disable-sjlj-exceptions * --enable-version-specific-runtime-libs --with-libiconv-prefix=/mingw * --with-libintl-prefix=/mingw --enable-libstdcxx-debug --with-isl=/mingw * --enable-libgomp --disable-libvtv --enable-nls * --disable-build-format-warnings * Thread model: win32 * gcc version 8.2.0 (MinGW.org GCC-8.2.0-5) * * C:\MinGW\bin>ld -v * GNU ld (GNU Binutils) 2.32 * */ #include #include /* For use with d > 0. */ #define DIV(n, d) ((n) >= 0 ? (n) / (d) : (-((-((n) + 1)) / (d)) - 1)) // Make it easy to experiment with different possibilities. // To exclude a possibility append a "q" to the "#define name". // "#define xx" has precedence over "#define xxx", // which has precedence over the "#define ax_*_*". // Show how "x" error can depend on whether a "printf" statement has a variable; // compile and run with only "#define xx", then with only "#define xxz". #define xx #define xxzq // Show how "a" and "x" errors can depend on the for loop or loops used; // compile and run with each of these "#define" only one at a time; // then compile and run with only "#define ax_4_2" and "#define ax_4_3"; // then compile and run with only "#define ax_3_1" and "#define ax_3_2". #define ax_4_11q #define ax_4_2q #define ax_3_1q #define ax_4_3q #define ax_3_2q #if defined(xx) || defined(xxz) static void qx(int jd) { double x = floor((jd - 1867216.25) / 36524.25); #ifdef xx printf("// qx:xx: jd %d; x %4.1f; z ;\n", jd, x); #else printf("// qx:xxz: jd %d; x %4.1f; z %d;\n", jd, x, 0); #endif } int main() { qx(1976789 - 146097); // Gregorian 0300-03-01 qx(1976789); // Gregorian 0700-03-01 return 0; } #else static void qax(int yy) { int y = yy * 100; double a = floor(y / 100.0); int ia = DIV(y, 100); int ib = 2 - ia + DIV(ia, 4); // Get Chronological Julian Day Number (CJDN) of Gregorian y-03-01. // 1721118 = CJDN of Julian Calendar 0000-03-01 // (1) calculate CJDN of Julian Calendar y-03-01 // (2) use "ib" to convert CJDN of Julian Calendar y-03-01 // to CJDN of Gregorian Calendar y-03-01 // "a" is used to calculate "b" // "x" is used to convert a CJDN to the calendar year, month, day int jd = 1721118 + y * 365 + DIV(y, 4) + ib; double x = floor((jd - 1867216.25) / 36524.25); int ix = DIV(jd * 4 - 7468865, 146097); printf("// yy %4d y %6d: a %5.1f %3d ? %2d", yy, y, a, ia, ((int) a) - ia); printf(" # jd %7d: x %5.1f %3d ? %2d", jd, x, ix, ((int) x) - ix); printf("\n"); } int main() { int yy; #ifdef ax_4_11 printf("//\n"); for (yy = -4; yy <= 11; yy++) qax(yy); #endif #ifdef ax_4_2 printf("//\n"); for (yy = -4; yy <= -2; yy++) qax(yy); #endif #ifdef ax_3_1 printf("//\n"); for (yy = -3; yy <= -1; yy++) qax(yy); #endif #ifdef ax_4_3 printf("//\n"); for (yy = -4; yy <= -3; yy++) qax(yy); #endif #ifdef ax_3_2 printf("//\n"); for (yy = -3; yy <= -2; yy++) qax(yy); #endif return 0; } #endif