The digits in the MGRS representation are given by multiplying x and y
by 1000000 (to convert to microns) and taking the floor of the result.
The expectation is that the multiplication and the floor produce
correctly rounded results. E.g.,
y=63.811z=1000000*yfloor(z)
should give 63811000. Evidentally with your setup, you're getting
63810999. This might be consequence of z being held in a register with
extra precision; so it's possible that declaring z to be volatile will
fix this. Perhaps g++ 10 needs a flag to modify this behavior. I can't
really explore this until I get my hands on a system with g++ 10. So I
propose to table this issue for now.
If you want to probe deeper with the information I've provided that
would be great.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Anonymous
Anonymous
-
2021-02-04
Using gdb, the value of y, entered by interpreting the string "63.811", is displayed as follows, after being loaded into an extended double precision register in the x87 floating point unit:
The question remains, of course, why these tests ever succeed without the adjustment. To answer that, I am in the fortunate situation that under different circumstances, namely
$ uname-aLinuxtn41.thorkilnaur.dk5.6.13-100.fc30.i686#1SMPFriMay1500:02:53UTC2020i686i686i386GNU/Linux
$ g++--versiong++(GCC)9.3.120200408(RedHat9.3.1-2)Copyright(C)2019FreeSoftwareFoundation, Inc.
Thisisfreesoftware; see the source for copying conditions. There is NOwarranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$
the tests actually do succeed.
So looking into the details under these Linux circumstances, I find, using gdb, that the values of y and y*1000000 are displayed as follows
The difference comes about when the floor() of y*1000000 is to be calculated: Under Cygwin, the floor() is calculated with inline code (the x87 frndint instruction with the RC rounding control set to "round down toward minus infinity"), resulting in the undesired value 63810999.
Under Linux, in contrast, a floor() function is called to perform this calculation. And the parameter of this function is passed in memory which means that the value held in the 80-bit extended precision floating point register must be stored as an ordinary 64-bit double prior to the call. And storing implies rounding to the shorter fraction of the 64-bit double.
The 80-bit value of y*1000000 is interpreted in binary as follows, grouping the bits into 1-bit sign, 15-bit exponent, 1-bit integer part, and 63-bit fraction (see the Intel manual "SECTION 8.2 X87 FPU DATA TYPES" for details about these floating point formats):
we can see what has happened: Rounding to the nearest 52-bit fraction has actually resulted in an integer, rather than a value slightly less than an integer. And hence, the floor() calculates the desired result. A lucky case, one may suppose.
In any case, the suggested adjustment works for the Linux circumstances as well: The tests still succeed with the adjustment applied.
Regards
Thorkil Naur
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
round will make the last digit 1. Check the documentation for MGRS::Forward.
I will also note that if you compile GeographicLib to use long double instead of double (with cmake -D GEOGRAPHICLIB_PRECISION=3 ...), then floor(63.811L * 1000000) gives 63811000 consistent with the behavior I expect for double.
So the "right" way to fix this problem is with something like
My hesitation about making this change is that the problem doesn't occur with g++ 10.2.1 on Linux. So presumably the failure you're seeing requires a rather specific set of conditions. Could you tell me what preprocessor symbols I can use to check for Cygwin?
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Anonymous
Anonymous
-
2021-02-05
Indeed, using round() instead of floor() was a mistake on my side, sorry about that.
Nevertheless, I hesitate to attribute the problem to any specific Cygwin circumstances. I am presently more inclined to say that we must accept both results (61811000 and 61810999) of calculating floor(y*1000000) with y=61.811: When excess precision is used, "it is unpredictable when rounding to the types specified in the source code takes place" (quoting the man g++ description of the -fexcess-precision=style option). And whether and when such rounding takes place clearly affects the outcome of the floor() calculation.
In an attempt to shed some light on this, I have distilled what seems to be the essence of the matter into the attached test program mfTest.cpp: Interpretation of a string into a floating point value, followed by the calculation of floor() of that value, after multiplying by 1000000.
Running this program:
$ uname-aLinuxtn41.thorkilnaur.dk5.6.13-100.fc30.i686#1SMPFriMay1500:02:53UTC2020i686i686i386GNU/Linux
$ g++--versiong++(GCC)9.3.120200408(RedHat9.3.1-2)Copyright(C)2019FreeSoftwareFoundation, Inc.
Thisisfreesoftware; see the source for copying conditions. There is NOwarranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ g++mfTest.cpp-O3-omfTest
$ ./mfTest61.811mfTest: 2021-Feb-0507.17mainT: sizeof(double)=8, prec=20megaFloor(t="double", s="61.811"): floor(61.810999999999999943*1000000)=61810999mainT: sizeof(longdouble)=12, prec=25megaFloor(t="long double", s="61.811"): floor(61.81099999999999999866773*1000000)=61811000
$
The precisions used when presenting the floating point values are excessive, of course. Nevertheless, I trust that presenting the floating point value resulting from the interpretation of the string "61.811" as "61.8109999..." indicates that the floating point value is strictly less than exactly 61.811. Hence, the double calculation could in this case be judged "correct" and the long double calculation "incorrect". But the situation is, in my judgement, that both results are equally acceptable.
The unpredictability of the calculations is illustrated by repeating the run of the program, now compiled with -O0 instead of -O3:
In this case, both the double calculation and the long double calculation could be judged "incorrect".
To avoid depending on unpredictabilities, one way would be to simply adjust a test such as GeoConvert17 to use --input-string "38n 500000 63.8115", adding a guard digit "5" to the input y value.
Running the GeographicLib-1.51 tests on Windows 7 with
the following failures happen:
Not a significant problem for me, to be sure.
Best regards
Thorkil Naur
The digits in the MGRS representation are given by multiplying x and y
by 1000000 (to convert to microns) and taking the floor of the result.
The expectation is that the multiplication and the floor produce
correctly rounded results. E.g.,
should give 63811000. Evidentally with your setup, you're getting
63810999. This might be consequence of z being held in a register with
extra precision; so it's possible that declaring z to be volatile will
fix this. Perhaps g++ 10 needs a flag to modify this behavior. I can't
really explore this until I get my hands on a system with g++ 10. So I
propose to table this issue for now.
If you want to probe deeper with the information I've provided that
would be great.
Using gdb, the value of y, entered by interpreting the string "63.811", is displayed as follows, after being loaded into an extended double precision register in the x87 floating point unit:
The value of y*1000000 is displayed as follows:
So it is understandable that floor()ing this value results in the undesired result, 63810999.
It appears that round() should be used rather than floor() in MGRS::Forward(...) as follows:
With this adjustment, the failing tests succeed.
The question remains, of course, why these tests ever succeed without the adjustment. To answer that, I am in the fortunate situation that under different circumstances, namely
the tests actually do succeed.
So looking into the details under these Linux circumstances, I find, using gdb, that the values of y and y*1000000 are displayed as follows
so exactly the same as under Cygwin.
The difference comes about when the floor() of y*1000000 is to be calculated: Under Cygwin, the floor() is calculated with inline code (the x87 frndint instruction with the RC rounding control set to "round down toward minus infinity"), resulting in the undesired value 63810999.
Under Linux, in contrast, a floor() function is called to perform this calculation. And the parameter of this function is passed in memory which means that the value held in the 80-bit extended precision floating point register must be stored as an ordinary 64-bit double prior to the call. And storing implies rounding to the shorter fraction of the 64-bit double.
The 80-bit value of y*1000000 is interpreted in binary as follows, grouping the bits into 1-bit sign, 15-bit exponent, 1-bit integer part, and 63-bit fraction (see the Intel manual "SECTION 8.2 X87 FPU DATA TYPES" for details about these floating point formats):
The 64-bit value stored for this is displayed as follows as a byte sequence, the least significant byte first:
This sequence is interpreted as follows, grouping into 1-bit sign, 11-bit exponent, 0-bit (implied) integer part, and 52-bit fraction:
Aligning the fractions
we can see what has happened: Rounding to the nearest 52-bit fraction has actually resulted in an integer, rather than a value slightly less than an integer. And hence, the floor() calculates the desired result. A lucky case, one may suppose.
In any case, the suggested adjustment works for the Linux circumstances as well: The tests still succeed with the adjustment applied.
Regards
Thorkil Naur
Indeed the problem will occur if
flooris inlined and the argumenty * mult_is not rounded to adouble.However your proposed solution, replacing
floorbyroundwon't work. Then you'll get the wrong answer forroundwill make the last digit1. Check the documentation forMGRS::Forward.I will also note that if you compile GeographicLib to use
long doubleinstead ofdouble(withcmake -D GEOGRAPHICLIB_PRECISION=3 ...), thenfloor(63.811L * 1000000)gives63811000consistent with the behavior I expect fordouble.So the "right" way to fix this problem is with something like
My hesitation about making this change is that the problem doesn't occur with
g++ 10.2.1on Linux. So presumably the failure you're seeing requires a rather specific set of conditions. Could you tell me what preprocessor symbols I can use to check for Cygwin?Indeed, using round() instead of floor() was a mistake on my side, sorry about that.
Nevertheless, I hesitate to attribute the problem to any specific Cygwin circumstances. I am presently more inclined to say that we must accept both results (61811000 and 61810999) of calculating floor(y*1000000) with y=61.811: When excess precision is used, "it is unpredictable when rounding to the types specified in the source code takes place" (quoting the man g++ description of the -fexcess-precision=style option). And whether and when such rounding takes place clearly affects the outcome of the floor() calculation.
In an attempt to shed some light on this, I have distilled what seems to be the essence of the matter into the attached test program mfTest.cpp: Interpretation of a string into a floating point value, followed by the calculation of floor() of that value, after multiplying by 1000000.
Running this program:
The precisions used when presenting the floating point values are excessive, of course. Nevertheless, I trust that presenting the floating point value resulting from the interpretation of the string "61.811" as "61.8109999..." indicates that the floating point value is strictly less than exactly 61.811. Hence, the double calculation could in this case be judged "correct" and the long double calculation "incorrect". But the situation is, in my judgement, that both results are equally acceptable.
The unpredictability of the calculations is illustrated by repeating the run of the program, now compiled with -O0 instead of -O3:
In this case, both the double calculation and the long double calculation could be judged "incorrect".
To avoid depending on unpredictabilities, one way would be to simply adjust a test such as GeoConvert17 to use --input-string "38n 500000 63.8115", adding a guard digit "5" to the input y value.
Using this idea, I have tried:
With these changes, all tests succeed on both Cygwin and Linux.
Best regards
Thorkil Naur