FORTRAN floating point accuracy problems

03272016, 03:02 PM
Post: #1




FORTRAN floating point accuracy problems
I'm getting some unacceptable results from rounding in old FORTRAN. I know there are some very talented mathematicians on the board here so I wonder if you guys would be willing to help me. To say I'm bad at math would be an understatement! I have no idea where things are going wrong.
I'm trying to accomplish something very simple in a language which has very little in the way of helpful functions to do what I want. Given a floating point input value for degrees, minutes, seconds in the form dd.mmsst... I want to produce the value in decimal degrees. This is no problem of math even for me. The problem is extracting the digits with only a truncation function and nothing else helpful in the repetoire as far as I can see. Here is the piece of code: Code:
and here are the results of some tests: Code:
You can see in test case #3 above there is a big problem extracting the minutes. The code winds up with 29 instead of 30 because of the truncation function. I find this kind of problem is coming up a lot in various calculations I'm doing and I don't know how to code around it. Can anybody please explain how I can rearrange the computations or do something differently to avoid this problem? Thanks. It ain't OVER 'till it's 2 PICK 

03272016, 05:15 PM
(This post was last modified: 03272016 07:52 PM by Dieter.)
Post: #2




RE: FORTRAN floating point accuracy problems
(03272016 03:02 PM)HP67 Wrote: I'm getting some unacceptable results from rounding in old FORTRAN. Sure. This is caused by the inexact binary representation of decimal numbers. For instance, the value "30.3" here obviously has a binary value that translates to 30.29999999999... in decimal notation. This again is interpreted as 30 degrees, 29 minutes and 99.99999... seconds, to that the final result is 30.5111... degrees. I translated your Fortran code to Excel VBA, and indeed the same problems occured. There are several ways to overcome the problem. For a quick and dirty solution, simply add a tiny amount to the ss value in the first line of code, say, 1E12 if 15 digits are used and the angle is below 1000° . This will cause a slight error (1E–10 seconds or 3 E–15 degrees) in the final result. A better solution may check if the mm.ss value is very close to the next higher integer, i.e. mm.999.... Since the seconds cannot exceed 59.9999, even a test whether the ss portion of mm.ss exceeds 0.6 may be sufficient. At least if the user does not enter nonstandard values like 30°40'85". Here are two solutions in VBA. The first one checks if the ss part of mm.ss exceeds 0.99999... In this case the minutes are rounded up to the next higher integer, which also means that the seconds are zero. Code: Function dms2dd(dms) The second solution simply increases mm by a tiny amount which is close to one ULP of the input. For 52bit binary arithmetics you may try 2^–40 or 1E–12. Code: Function dms2dd(dms) A third solution would simply round the argument (i.e. dms) to the max. possible number of decimals. For instance, if dms=30.3000 and you're working with 15 digit precision, the 0.3000 part cannot have more than 13 digits following the decimal point (15 minus two for the integer part). So the fractional part cannot have 15 digits like 0.29999 99999 99997 (which causes the roundoff problem), but it has to be 0.30000 00000 000. In the final code this rounding procedure would replace the line with the "tiny" variable in the second solution: Code: Function dms2dd(dms) Maybe this is the preferred solution. BTW, your result table seems to be wrong. For instance, for dms(71.0030) the expected value is 71 + 30/3600 = 71.00833333... and not 71.00830. This is also true for several other results. Dieter 

03272016, 10:47 PM
Post: #3




RE: FORTRAN floating point accuracy problems
For this kind of problem, I find first dealing with the sign and then getting into an integer domain as fast as possible and then applying corrections there is usually the best approach.
You've got integer degrees already, multiply by 60 and get integer minutes, then again by sixty to get seconds (which will have to stay real). Then sanity check everything to deal with the corner cases. You'll end up with a pile of code along the lines of: Code: if seconds .ge. 60 then  Pauli 

03272016, 11:13 PM
(This post was last modified: 03272016 11:13 PM by Dieter.)
Post: #4




RE: FORTRAN floating point accuracy problems
(03272016 10:47 PM)Paul Dale Wrote: For this kind of problem, I find first dealing with the sign and then getting into an integer domain as fast as possible and then applying corrections there is usually the best approach. OK, but how do you prevent the roundoff problems due to binary arithmetics? Consider the example with 33,3000. Isolating the minutes leads to 29,99999... instead of 30, which then is interpreted as 29 minutes and 99,99999... seconds, i.e. 30 minutes and 39,9999... seconds. So some kind of rounding has to be applied. Dieter 

03282016, 03:48 AM
Post: #5




RE: FORTRAN floating point accuracy problems
The most recent time I did this was using decimal arithmetic and it took several goes.
Yes, rounding will be required. At least it is easy enough to detect the errant case  when the seconds value is greater than sixty there is possibly a problem. Figuring out where and how to round is another matter entirely  this is getting close to fudging the answer.  Pauli 

03282016, 09:29 AM
Post: #6




RE: FORTRAN floating point accuracy problems
Hi Dieter, thanks so much for your time and explanations. I'm not sure I got it all. I understand floating point representation is inexact. But I don't know how to fix it.
I did see in some cases adding a small correction would help but when I added it in all cases it made the results much worse overall. Below when you said to add the correction in the first line, did you mean the first line of my code (I think there it is already too late) or in the first line of your example code? (03272016 05:15 PM)Dieter Wrote: I translated your Fortran code to Excel VBA, and indeed the same problems occured. There are several ways to overcome the problem. For a quick and dirty solution, simply add a tiny amount to the ss value in the first line of code, say, 1E12 if 15 digits are used and the angle is below 1000° . This will cause a slight error (1E–10 seconds or 3 E–15 degrees) in the final result. (03272016 05:15 PM)Dieter Wrote: A better solution may check if the mm.ss value is very close to the next higher integer, i.e. mm.999.... Since the seconds cannot exceed 59.9999, even a test whether the ss portion of mm.ss exceeds 0.6 may be sufficient. At least if the user does not enter nonstandard values like 30°40'85". Yes, the numbers are generated internally so they should be wellformed. I think I understand your second suggestion but I'll have to look at your code and my code further. (03272016 05:15 PM)Dieter Wrote: The second solution simply increases mm by a tiny amount which is close to one ULP of the input. For 52bit binary arithmetics you may try 2^–40 or 1E–12. I will have to check. The computer I'm using does not use IEEE binary, it uses a proprietary hexadecimal format and I don't remember how many bits of precision are carried in the mantissa. I'll get back to you on that. How do you choose a correction value close to one ULP? You said for 52 bit binary we should use 2^40 or 1E12. What is the connection between these values and the bit precision? (03272016 05:15 PM)Dieter Wrote: A third solution would simply round the argument (i.e. dms) to the max. possible number of decimals. For instance, if dms=30.3000 and you're working with 15 digit precision, the 0.3000 part cannot have more than 13 digits following the decimal point (15 minus two for the integer part). So the fractional part cannot have 15 digits like 0.29999 99999 99997 (which causes the roundoff problem), but it has to be 0.30000 00000 000. In the final code this rounding procedure would replace the line with the "tiny" variable in the second solution: I don't have a rounding function and I don't understand how to write one. Is this solution still possible? Is it possible to write a general purpose rounding function? All I have is trig functions, integer (which does not simply take the integer portion and store in the 8 byte floating point value but rather it demotes the 8 byte floating point integer value to a 4 byte integer type. (03272016 05:15 PM)Dieter Wrote: BTW, your result table seems to be wrong. For instance, for dms(71.0030) the expected value is 71 + 30/3600 = 71.00833333... and not 71.00830. This is also true for several other results. Thanks, that will teach me to rely on test data from the web. I should have checked this with an HP calculator... Thank you very much. It ain't OVER 'till it's 2 PICK 

03282016, 09:33 AM
Post: #7




RE: FORTRAN floating point accuracy problems
Hi and thanks for helping with this.
(03272016 10:47 PM)Paul Dale Wrote: For this kind of problem, I find first dealing with the sign and then getting into an integer domain as fast as possible and then applying corrections there is usually the best approach. That would be ideal but here the truncation and loss of precision occurs during extracting the minutes and seconds. And once in integer values the results of division also produce truncated values. So, I can't see how to do this. (03272016 10:47 PM)Paul Dale Wrote: You've got integer degrees already, multiply by 60 and get integer minutes, then again by sixty to get seconds (which will have to stay real). Then sanity check everything to deal with the corner cases. If I multiply by 60 to move the minutes to the left of the decimal there is already some loss of precision for certain values like 30 minutes for example. I don't know how to extract the minutes or seconds. This is exactly the problem. It ain't OVER 'till it's 2 PICK 

03282016, 01:29 PM
(This post was last modified: 03282016 01:34 PM by Dieter.)
Post: #8




RE: FORTRAN floating point accuracy problems
(03282016 09:29 AM)HP67 Wrote: I did see in some cases adding a small correction would help but when I added it in all cases it made the results much worse overall. Below when you said to add the correction in the first line, did you mean the first line of my code (I think there it is already too late) or in the first line of your example code? The first line (after the declarations) of your code, the one that reads SS = (DMS  IDINT(DMS)) * 100.0D0. (03282016 09:29 AM)HP67 Wrote: Yes, the numbers are generated internally so they should be wellformed. That's good to know. However, the third and preferred solution should handle either case. (03282016 09:29 AM)HP67 Wrote: I will have to check. The computer I'm using does not use IEEE binary, it uses a proprietary hexadecimal format and I don't remember how many bits of precision are carried in the mantissa. I'll get back to you on that. How do you choose a correction value close to one ULP? You said for 52 bit binary we should use 2^40 or 1E12. What is the connection between these values and the bit precision? Let's assume the working precision is 15 decimal digits (like IEEE754 double precision). If 23 decimal digits are kept for the integer portion of the argument (whole degrees), the remaining 1213 hold the mmss part. So the last digit of mmss is (roughly) 1E12. Or (base 2) about 2^–39. The idea is to get this value large enough to do the rounding but small enough to minimize its mathematical impact, as it affects the last digit(s) of the final result. If you want to implement this idea I would recommend a value for "tiny" that is based on the argument dms, e.g. tiny=dms/1E14 for 15digit working precision. So this is how it is supposed to work: Code: entered argument dms: But the third approach should work even better: (03282016 09:29 AM)HP67 Wrote: I don't have a rounding function and I don't understand how to write one. Is this solution still possible? Is it possible to write a general purpose rounding function? All I have is trig functions, integer (which does not simply take the integer portion and store in the 8 byte floating point value but rather it demotes the 8 byte floating point integer value to a 4 byte integer type. Sure. Write your own rounding function. ;) CAVE: I DO NOT HAVE ANY EXPERIENCE WITH FORTRAN! But from what I understand there is a NINT resp. IDNINT function that rounds a value to its closest integer. For a function that rounds to n digits, simply multiply the argument by 10^n, round this to the next closest integer and divide by 10^n again. Here is the general idea for rounding a positive value X to D decimal places: Code: THIS IS NOT FORTRAN CODE, IT JUST SHOWS THE GENERAL IDEA! Please note that even R is rounded since we cannot be sure that calculating 10^13 will return exactly 10,000,000,000,000. So just to be sure... Once you know the number of digits to round (13 in the following example), rounding is easy. Simply define the factor R somewhere before (or use a constant  don't know how this is done in Fortran) and the rest is trivial: Code: R = 1.0D13 A closertoperfect solution would round to a variable number of digits, depending on the magnitude of the dms argument. But 13 digits works well in my Excel VBA implementation, at least for angles up to 999 degrees. Now try this and post your improved code here. Dieter 

03282016, 02:20 PM
(This post was last modified: 03282016 03:37 PM by HP67.)
Post: #9




RE: FORTRAN floating point accuracy problems
Thanks again, Dieter. Modern Fortran has many more functions including ceil, floor, and nearest integer but I'm using FORTRAN IV which has only truncation. I want to look over the hardware doc before I decide what to do. I don't understand the floating point instructions sine I never used them. I do remember there is some kind of hardware rounding function. It might be better to write whatever primitives I need in assembler and call them from FORTRAN. I'll post again after I have a better idea what hardware support there is.
What I should have mentioned: Modern Fortran has IEEE support and depending on the compiler that should include all IEEE rounding modes, etc. But, a long time ago (F90, possibly F77) Fortran added a feature called internal reads where you can "read" and "write" into variables without any I/O. Since FORTRAN always rounds on input/output this was a way to get a standard rounding mode for whatever precision you wanted. However it's not available as far back as I am working, so... It ain't OVER 'till it's 2 PICK 

03302016, 01:10 PM
Post: #10




RE: FORTRAN floating point accuracy problems
Dieter, something is wrong here.
The rounding function works to some extent when I specify a smallish value for decimal digits. If I round to 4 or 5 places it seems to work better than 13. I don't know why. In any case, the results when rounded are worse overall than without the rounding calculation. In many cases, it appears to be rounding down instead of up. I think the problem is caused by the lack of a nearest integer function but I am not sure. The IDINT function is a double precision (8 byte) that simply returns the integer part like IP in RPL. The MMSS value is for the following line (out of order because I added a WRITE statement to the subroutine to debug it. For example in the first test the input value is 89.11150 but MMSS is extracted and rounded to 13 digits according to your example above and comes out to 11.1499999649823) Code:
If I specify rounding to 4 places the results improve a lot but are still worse than when unrounded: Code:
Do you have any idea what's going on? Why should rounding to more places cause it to round down instead of up? Here is the code at this point: Code:
Code:
It ain't OVER 'till it's 2 PICK 

03302016, 03:41 PM
(This post was last modified: 03302016 03:55 PM by Dieter.)
Post: #11




RE: FORTRAN floating point accuracy problems
(03302016 01:10 PM)HP67 Wrote: Do you have any idea what's going on? Why should rounding to more places cause it to round down instead of up? Yes, I think I know what's going on here: Your Fortran seems to support a precision of about ten decimal digits. That's why rounding to 13 digits does not help. I think your DROUND function is fine, just replace the 4 digit to 8 digit rounding. Or, if you are sure that nothing beyond 1/100 second occurs, simply set it to 6 digit rounding. You may test this easily: simply write a short program that prints a periodic value like 1/3 or 1/9 with 15 digits. For instance, simply let the program print 1.0D0/3.0D0 with 15 decimal places and see what you get. If it's something like 0,33333 33334 8855... you know that 31 bits or 910 decimal places is all you got. ;) (03302016 01:10 PM)HP67 Wrote: If I specify rounding to 4 places the results improve a lot but are still worse than when unrounded: Sure. Rounding to 4 places will work fine if 0.mmss does not have more than 4 digits, i.e. no fractions of a second. And indeed for these cases your results are fine, only the "EXPECT" values are wrong again: dms2dd(71,0030) = 71,0083333.... and dms2dd(0,4949) = 0,83027777...., or dms2dd(65,11010) = 65,18361111... which is exactly what your program returns. The "expected" values 71,0080, 0,83020 and 65.18600 are WRONG! Actually all cases with not more than 4 decimals are perfectly correct. Even those with 5 places (i.e. fractions of a second) are exact if the argument is rounded to full seconds (4 digits). For instance, dms2dd(13,07244) actually evaluates dms2dd(13,0724), so mm.ss = 7,2400000 and the result for this is 13,12333. Adjust the roundoff function to 5 places (or 6, or 8) and these will return exact results as well. Now please try rounding to 6 or 8 digits and report what you got. And, more important, check your results with a reliable reference, e.g. an HP. ;) Dieter 

03302016, 04:21 PM
Post: #12




RE: FORTRAN floating point accuracy problems
(03302016 03:41 PM)Dieter Wrote: Yes, I think I know what's going on here: Your Fortran seems to support a precision of about ten decimal digits. That's why rounding to 13 digits does not help. I think your DROUND function is fine, just replace the 4 digit to 8 digit rounding. Or, if you are sure that nothing beyond 1/100 second occurs, simply set it to 6 digit rounding. According to the doc it's 14 hexadecimal digits which is supposed to be approximately equivalent to 16.8 decimal digits. I don't understand how it's encoded yet so I can't verify this. But it's IBM and the doc is never wrong. It's just not very helpful in this manual. I have another couple places I want to check. So, the 13 digits you mentioned should have been perfect. (03302016 03:41 PM)Dieter Wrote: You may test this easily: simply write a short program that prints a periodic value like 1/3 or 1/9 with 15 digits. For instance, simply let the program print 1.0D0/3.0D0 with 15 decimal places and see what you get. If it's something like 0,33333 33334 8855... you know that 31 bits or 910 decimal places is all you got. ;) Yes, and more to the point I'll run through a loop printing the output of the rounding function for 1/3 for 1..16 digits. [face palm] why didn't I think of this.. (03302016 01:10 PM)HP67 Wrote: If I specify rounding to 4 places the results improve a lot but are still worse than when unrounded: (03302016 03:41 PM)Dieter Wrote: Sure. Rounding to 4 places will work fine if 0.mmss does not have more than 4 digits, i.e. no fractions of a second. And indeed for these cases your results are fine, only the "EXPECT" values are wrong again: dms2dd(71,0030) = 71,0083333.... and dms2dd(0,4949) = 0,83027777...., or dms2dd(65,11010) = 65,18361111... which is exactly what your program returns. The "expected" values 71,0080, 0,83020 and 65.18600 are WRONG! They should be correct in the latest run. I checked all the calculations on an HP 48 and made a few changes before I started working on this today. I think the problem is lack of precision in the printout and the expected values are really correct but I'll read your post again and check it against my test. I'm only displaying 5 decimal places and the output routines round. I'll extend the size of the input values in the next printout. I was only using 5 decimal places because that will be enough for my end use of this routine. For debugging we should probably show more. (03302016 03:41 PM)Dieter Wrote: Even those with 5 places (i.e. fractions of a second) are exact if the argument is rounded to full seconds (4 digits). For instance, dms2dd(13,07244) actually evaluates dms2dd(13,0724), so mm.ss = 7,2400000 and the result for this is 13,12333. Adjust the roundoff function to 5 places (or 6, or 8) and these will return exact results as well. I don't think this is working. Although the input, expected, and actual results are shown they are rounded by the output routines. The difference value is valid at least to 5 decimal places rounded. In many cases the difference is nonzero within that precision so I think there is a still a problem. (03302016 03:41 PM)Dieter Wrote: Now please try rounding to 6 or 8 digits and report what you got. Yes, I'll run the test and post again tomorrow hopefully. (03302016 03:41 PM)Dieter Wrote: And, more important, check your results with a reliable reference, e.g. an HP. ;) If I had a fullsize keyboard on my HP calcs I would probably throw my PC out the window! Thank you very much for helping me with this. It's so frustrating not to get it when it comes to math. It ain't OVER 'till it's 2 PICK 

03302016, 06:38 PM
(This post was last modified: 03302016 06:48 PM by Dieter.)
Post: #13




RE: FORTRAN floating point accuracy problems
(03302016 04:21 PM)HP67 Wrote: According to the doc it's 14 hexadecimal digits which is supposed to be approximately equivalent to 16.8 decimal digits. I don't understand how it's encoded yet so I can't verify this. But it's IBM and the doc is never wrong. I am quite confident that there's not more than 31 or 32 bit. Take a look at the slightly inaccurate MMSS values you listed. These can be accurately represented in binary notation with 31 bits or less: Code: 30.299999999348074 = 65,068,754,533 / 2^31 32 bits mean 9,3 decimal digits, which matches the results. Since 2 digits are used for the integer portion, rounding the mmss part to 7 decimals shoud be OK. (03302016 04:21 PM)HP67 Wrote: Yes, and more to the point I'll run through a loop printing the output of the rounding function for 1/3 for 1..16 digits. [face palm] why didn't I think of this.. Sounds like a good idea. Ideally the number of rounded digits should depend on the integer part of the dms argument. 2 digits => round to 7 decimals. 3 digits => round to 6 decimals. 0 digits => round to 9 decimals. More generally: round to 8 – floor(log10(dms)) decimals. In other words: round to 9 significant (!) decimal digits. The dms argument obviously is not represented more accurately anyway. (03302016 04:21 PM)HP67 Wrote: They should be correct in the latest run. I checked all the calculations on an HP 48 and made a few changes before I started working on this today. I think the problem is lack of precision in the printout and the expected values are really correct but I'll read your post again and check it against my test. Well, 71°00'30" is 71 + 30/3600 = 71,008333333... degrees and not 71,0080°. There's not much room for different interpretations. ;) (03302016 04:21 PM)HP67 Wrote:(03302016 03:41 PM)Dieter Wrote: Even those with 5 places (i.e. fractions of a second) are exact if the argument is rounded to full seconds (4 digits). For instance, dms2dd(13,07244) actually evaluates dms2dd(13,0724), so mm.ss = 7,2400000 and the result for this is 13,12333. Adjust the roundoff function to 5 places (or 6, or 8) and these will return exact results as well. I do. :) Dieter 

03302016, 08:32 PM
Post: #14




RE: FORTRAN floating point accuracy problems
(03302016 06:38 PM)Dieter Wrote: 32 bits mean 9,3 decimal digits, which matches the results. Since 2 digits are used for the integer portion, rounding the mmss part to 7 decimals shoud be OK. Dieter  what is the meaning of "9,3 decimal digits" ? I get that due to internal representation inaccuracies you can only get somewhere between 9 and 10 reliable decimal digits, but what is the meaning or interpretation of the ",3" ? These threads are interesting to follow, I'm always learning something new, but I've not seen this usage before. Or maybe have simply forgotten it in the last 30+ years since doing this stuff on a daily basis. Bob Prosperi 

03302016, 09:09 PM
(This post was last modified: 03302016 09:17 PM by Dieter.)
Post: #15




RE: FORTRAN floating point accuracy problems
(03302016 08:32 PM)rprosperi Wrote: Dieter  what is the meaning of "9,3 decimal digits" ? I get that due to internal representation inaccuracies you can only get somewhere between 9 and 10 reliable decimal digits, but what is the meaning or interpretation of the ",3" ? In mathematics, results here and there may not exactly match real life. ;) The meaning is exactly as you suggested: "somewhere between 9 and 10 reliable decimal digits". Say, 9 digits with an error less than 1 ULP. I think it should be about 0,5 ULP in this case, as 2^–31 is approx. 4,6 E–10. Or 2^–32 = 2,3 E–10 for 32 bits. Which BTW should translate to 9,6 digits (sorry). (03302016 08:32 PM)rprosperi Wrote: These threads are interesting to follow, I'm always learning something new, but I've not seen this usage before. Or maybe have simply forgotten it in the last 30+ years since doing this stuff on a daily basis. Specifying a number of valid digits with fractional values is quite common. For instance in discussions on the accuracy of approximation methods. I remember such notations ("10,7 digits") in a thread on various methods for the Gamma function. In that case it was simply the base10 log of the relative error. ;) Dieter 

03312016, 09:00 AM
Post: #16




RE: FORTRAN floating point accuracy problems
(03302016 06:38 PM)Dieter Wrote: I am quite confident that there's not more than 31 or 32 bit. Take a look at the slightly inaccurate MMSS values you listed. These can be accurately represented in binary notation with 31 bits or less: I certainly cannot and am not arguing that your observation is not correct. Still, there is 56 bits of precision here for an 8 byte real. 56 bits for the fraction, 1 for the sign, and 7 for the exponent. I'm testing on two different compilers on two different hardware architectures (although the former being emulated on the latter) and both get similar results with what is supposed to be 8 byte reals but totally different data formats. I suspect an error in my code like an incorrect default declaration being used (resulting in a 4 byte real) but I can't find it. This should not be happening. It has been bothering me all along that the results should be this far off. Something is definitely wrong here. (03302016 04:21 PM)HP67 Wrote: Yes, and more to the point I'll run through a loop printing the output of the rounding function for 1/3 for 1..16 digits. [face palm] why didn't I think of this.. Results shown below. (03302016 04:21 PM)HP67 Wrote: They should be correct in the latest run. I checked all the calculations on an HP 48 and made a few changes before I started working on this today. I think the problem is lack of precision in the printout and the expected values are really correct but I'll read your post again and check it against my test. (03302016 06:38 PM)Dieter Wrote: Well, 71°00'30" is 71 + 30/3600 = 71,008333333... degrees and not 71,0080°. There's not much room for different interpretations. ;) What I was saying is that the output shown was rounded by the FORTRAN output routines. I increased the length of the output format and you can see this is not a real problem. There was one error remaining in the last test (dyslexia on my part) but that is now fixed. (03302016 04:21 PM)HP67 Wrote: I don't think this is working. (03302016 06:38 PM)Dieter Wrote: I do. :) Please see below. With rounding to 5 places some of the values are significantly off. I also tried using 9 places and 10 and it didn't help. Code:
It ain't OVER 'till it's 2 PICK 

03312016, 12:39 PM
Post: #17




RE: FORTRAN floating point accuracy problems
(03312016 09:00 AM)HP67 Wrote: [... a lot ] Thanks for all the input. Here and there I see clearer now.
This is the point where I fear I cannot be of much help any longer, as further debugging requires knowledge in Fortran. But maybe there are some experts here. Dieter 

03312016, 12:51 PM
(This post was last modified: 03312016 01:01 PM by HP67.)
Post: #18




RE: FORTRAN floating point accuracy problems
Hi and thank you very much. I was out for a walk and I realized the problem was staring me in the face. I was hoping I could update the thread before you saw it. Yes, it's the IDINT function and it happens in DROUND as you suspected.
FORTRAN IV's IDINT returns a 4 byte integer. There is no 8 byte INT() function in FORTRAN IV. I see three possible approaches: 1) If you can define ROUND in terms of MOD (there is an 8 byte mod function available) or some other way that doesn't require taking the integer portion with the IDINT function then we can fix the rounding routine. The magnitude of the dd.mmss or hh.mmss data itself never gets close to the 2^31 signed integer limit. 2) If you know an algorithm we could substitute for the existing IDINT function. If we can go through some process to find the integer portion we can do it in 8 byte arithmetic. 3) If I can figure out the instruction and data formats for this I should be able to write an 8 byte truncation function and I wouldn't need IDINT. Some further updates on specific points: (03312016 12:39 PM)Dieter Wrote: But from what I read, mixing several data types in a calculation can cause weird results. This is generally safe if you pay attention and specify the correct constants and don't use probably the one builtin that downgrades precision, IDINT! Certainly the value of the RHS is promoted to the type of the LHS. That doesn't always help, obviously. (03312016 12:39 PM)Dieter Wrote: R = IDINT(10.0D0 ** P + 0.5D0) The exponent is a 4 byte integer I checked earlier that the result of an 8 byte real raised to any legal variable is still an 8 byte real. That part is fine since we're using integral powers less than 2^31 here. What is the precision of R? R is a double precision (8 byte) floating point value that was promoted from a single precision 4 byte integer :( It ain't OVER 'till it's 2 PICK 

03312016, 01:39 PM
Post: #19




RE: FORTRAN floating point accuracy problems
(03312016 12:51 PM)HP67 Wrote: FORTRAN IV's IDINT returns a 4 byte integer. There is no 8 byte INT() function in FORTRAN IV. That's the problem. (03312016 12:51 PM)HP67 Wrote: I see three possible approaches: I see a fourth one: simply do it all in integer arithmetics. Right at the beginning the dms input is converted to an integer. If the largest 4bit integer is 2147483647, the largest dms input is 214,7483647. So we just multiply the input by 10.000.000 and use that as an integer. I tried this in Excel VBA and it works well: Code: Function dms2dd(dms as double) Note: the "\" operator means integer division. The "#" suffix marks a variable or value as double precision. I think you would use 60.0D0 and 3.6D6. The major limitation of this method is the domain of the dms argument. The largest angle is 214 degrees and the precision of dms is 0.001 seconds. In other words, everything beyond the 7th decimal is neglected. This way the result has 56 valid decimals (0,001/3600 =1,7E–6). You may check if Fortran truncates or rounds x when it is converted in the first line (x=dms*10000000). Maybe you will have to code this as X = IDINT(DMS*1.0D7+0.5D0). What do you think? Dieter 

03312016, 02:45 PM
(This post was last modified: 03312016 02:55 PM by HP67.)
Post: #20




RE: FORTRAN floating point accuracy problems
Thanks very much.
Mostly I'm using this to convert decimal time back and forth. .001 seconds should be plenty of resolution. I'll be happy with a reliable hh.mmsst for now. I think a rounding function and a double precision truncation function would be more generally useful but I'll study what you came up with. I can anticipate a few problems but fortunately I am almost always wrong about math so maybe this time it will work in my favor. FORTRAN IV doesn't have any integer division except for integer types and doesn't have 8 byte integers. When we get back to dividing again won't we go through the same or similar conversions/promotions as before? So I am expecting problems. Again, I'm probably wrong. I'll play around with your suggestion. Eventually I would like to come up with an 8 byte truncation function for reals. The rounding function is very handy also. In understanding your declarations, Excel VBA long means 8 byte int? No comparable type in FORTRAN IV. Will a 4 byte int work? I think you explained that it will, so I don't understand why you chose long in your code. It ain't OVER 'till it's 2 PICK 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)