FORTRAN floating point accuracy problems

03312016, 07:30 PM
(This post was last modified: 04012016 06:19 AM by Dieter.)
Post: #21




RE: FORTRAN floating point accuracy problems
(03312016 02:45 PM)HP67 Wrote: 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. OK, so ddd actually is hhh. ;) What is the largest value you have to handle? Also, possibly two more digits could be obtained if the integer and fractional part of dms would be handled separately. But if 0.001 s is fine, the posted code is an easy solution. (03312016 02:45 PM)HP67 Wrote: FORTRAN IV doesn't have any integer division except for integer types and doesn't have 8 byte integers. That's exactly how the VBA code works. "Long" is a signed fourbyte integer, and "\" is what Fortran does with a simple "/" division of two integers (12/5 = 2). (03312016 02:45 PM)HP67 Wrote: In understanding your declarations, Excel VBA long means 8 byte int? No, that's the same four bytes, i.e. integers within ±2^31. You should be able to translate this more or less 1:1 into Fortran. "Long" just means "Long integer", compared to a standard integer within ±32768. Dieter 

04012016, 06:15 AM
(This post was last modified: 04012016 06:32 AM by Dieter.)
Post: #22




RE: FORTRAN floating point accuracy problems
(03312016 07:30 PM)Dieter Wrote: Also, possibly two more digits could be obtained if the integer and fractional part of dms would be handled separately. This could be a possible solution then: Code: Function dms2dd(dms As Double) Now nine decimals of dms are used for the conversion, i.e. the resolution is 0.00001 s. There is (almost) no limit for dms (max. ~2 billion hours). The result should have eight valid digits after the decimal point. Please note that VBA automatically rounds the result of "mmss = (dms  ddd) * 1000000000" to the nearest integer. You may have to code this as MMSS = IDINT((DMS–DDD)*1.0D9 + 0.5D0). Dieter 

04012016, 03:12 PM
Post: #23




RE: FORTRAN floating point accuracy problems
Thanks, can't respond now, will get back to you next week. Have a great weekend!
It ain't OVER 'till it's 2 PICK 

04032016, 03:49 PM
(This post was last modified: 04032016 05:08 PM by HP67.)
Post: #24




RE: FORTRAN floating point accuracy problems
Still having problems with specific input values. I tried my own version before I looked at yours and I tried yours also. Unsurprisingly yours worked a little better. But we are still having problems with specific values. BTW the output is shown to 7 decimal places but the input and expected values are only to 6 values. I showed extra digits since last time it seemed important.
Aside from specific problematic values (.15, .3) this works better than required. So, in that sense much progress. Thanks. Code:
This floating point hardware implementation uses a base 16 internal representation. Do you have any idea how we might change the algorithm to use powers of 16 instead of powers of 10 or if that has any chance of working better? I want to figure out how to code a double precision (8 byte) truncation function as I mentioned earlier. Your rounding routine is really great. If my implementation wasn't hampered by the existing 4 byte truncation function maybe all this would work a lot better? I have figured out how to convert the internal representation to decimal but there is some loss of precision on my HP calcs since they only carry 12 significant digits. I have not figured out enough about the math behind this to understand what I need to do to truncate the decimal portion while it's still in the internal base 16 format. I'm sure anybody mathematically inclined could see it right away but I'm still banging my head. Update: I see an instruction that could help with the the 8 byte truncation function. Will update on that in the next few days hopefully. It ain't OVER 'till it's 2 PICK 

04032016, 05:13 PM
(This post was last modified: 04032016 05:17 PM by Dieter.)
Post: #25




RE: FORTRAN floating point accuracy problems
(04032016 03:49 PM)HP67 Wrote: Still having problems with specific input values. I tried my own version before I looked at yours and I tried yours also. Unsurprisingly yours worked a little better. But we are still having problems with specific values. The problem seems to be the integer values mmss, mm and ss which are truncated so that a true result of mmss in case 3 may actually be 29.999999... instead of 30, which is then truncated to 29 minutes and 99.999... seconds. In VBA this problem does not occur as an assignment to an integer variable automatically rounds a real up or down to the next integer. Adding some explicit rounding should fix that: You could replace line 0006 MMSS = (DMS  DDD) * 1000000000 with MMSS = (DMS  DDD) * 1000000000 + 0.5D0 or even better MMSS = IDINT((DMS  DDD) * 1000000000 + 0.5D0) And once again the "expected" values are wrong. For instance Code: TEST DMS EXPECT RESULT DIFF ...should actually read... Code: TEST DMS EXPECT RESULT DIFF All these results are dead on. Actually only case 2, 3, 9 and 14 are off, and I think this can be fixed by adding the suggested integer rounding. The second program version has essentially the same issues: Line 0005 IDMS = DMS * 1.0D7 should better be rounded to IDMS = DMS * 1.0D7 + 0.5D0 or even better IDMS = IDINT(DMS * 1.0D7 + 0.5D0) (04032016 03:49 PM)HP67 Wrote: I want to figure out how to code a double precision (8 byte) truncation function as I mentioned earlier. Your rounding routine is really great. If my implementation wasn't hampered by the existing 4 byte truncation function maybe all this would work a lot better? Most probably, yes. Dieter 

04032016, 05:17 PM
Post: #26




RE: FORTRAN floating point accuracy problems
Again, I am only carrying 6 decimal places of expected results although I am displaying 7.
Thanks, I'll try your suggestions. But first I'm going to try to write an assembler subroutine for DINT. Then I can use your rounding routine which would really be great. It ain't OVER 'till it's 2 PICK 

04032016, 05:21 PM
(This post was last modified: 04032016 05:22 PM by Dieter.)
Post: #27




RE: FORTRAN floating point accuracy problems
(04032016 05:17 PM)HP67 Wrote: Again, I am only carrying 6 decimal places of expected results although I am displaying 7. This means that "differences" in the 7th place actually are not present. ;) Maybe you can update the "expect" values to seven decimals. (04032016 05:17 PM)HP67 Wrote: Thanks, I'll try your suggestions. But first I'm going to try to write an assembler subroutine for DINT. Then I can use your rounding routine which would really be great. Hmmm... I would love to know whether the suggested rounding fixes the issue. Maybe you can add these little changes and report the results... ?) Dieter 

04042016, 09:02 AM
Post: #28




RE: FORTRAN floating point accuracy problems
Can I chime in?
The following formulas, courtesy of John H. Meyers on comp.sys.hp48 do not require any fiddling at all: HR2HMS T := HR  FP(HR)*0.4; HMS := T  FP(T*100)*0.004; HMS2HR T := HMS + FP(HMS*100)/150; HR := T + FP(T)/1.5; Tried them on Free42S Binary. Both work like a charm. Code: 00 { 55Byte Prgm } Cheers, Werner 

04042016, 03:14 PM
Post: #29




RE: FORTRAN floating point accuracy problems
(04032016 05:21 PM)Dieter Wrote:(04032016 05:17 PM)HP67 Wrote: Again, I am only carrying 6 decimal places of expected results although I am displaying 7. Yes I know. Therefore I was not troubled by the diff values ;) (04032016 05:21 PM)Dieter Wrote: Maybe you can update the "expect" values to seven decimals. I'll have to look at it. It could be with 7 decimal places I can't fit 4 values on one line of FORTRAN code... (04032016 05:17 PM)HP67 Wrote: Thanks, I'll try your suggestions. But first I'm going to try to write an assembler subroutine for DINT. Then I can use your rounding routine which would really be great. (04032016 05:21 PM)Dieter Wrote: Hmmm... I would love to know whether the suggested rounding fixes the issue. Maybe you can add these little changes and report the results... ?) The good news is your rounding routine works perfectly up to 15 decimal places (possibly more, I will check) with my brand new shiny double precision (8 byte) DINT routine in assembler. I didn't have time to go back to our old versions of the driver code to see if it works better. I'll get back to that as soon as I can. Thanks for the help! I am really excited about having the rounding routine in my FORTRAN IV toolbag. I have wanted one for a long time and just couldn't figure out how to do it. I do have some rounding code I use for assembler control blocks but I just couldn't extrapolate it to real numbers. It ain't OVER 'till it's 2 PICK 

04042016, 03:19 PM
Post: #30




RE: FORTRAN floating point accuracy problems
(04042016 09:02 AM)Werner Wrote: Can I chime in? Hi. I know how to do the conversions by hand and with a calculator, etc. I need this for some code I am writing in circa 1974 FORTRAN IV. The math isn't the issue. The problem is dealing with loss of precision caused by the floating point hardware implementation which as far as I know isn't a problem on HP calculators or calculators generally since they use decimal arithmetic. It ain't OVER 'till it's 2 PICK 

04042016, 04:38 PM
Post: #31




RE: FORTRAN floating point accuracy problems
That's why I specifically tested them on Free42 Binary, which uses the binary floating point of the computer it's running on... all your examples ran perfectly.
Another issue is that since HH.MMSSssss cannot be correctly represented in binary floatingpoint, you shouldn't be using it? Use HHMMSS.ssss instead.. Cheers, Werner 

04042016, 04:49 PM
Post: #32




RE: FORTRAN floating point accuracy problems
The hardware implementation I'm using is not binary but I see where you are coming from and thanks for the effort. We have tried a few things and so far have been hampered by a few missing functions. We have got that mostly resolved. I just haven't had time to test again with the old versions of code that should work now.
It ain't OVER 'till it's 2 PICK 

04042016, 06:42 PM
Post: #33




RE: FORTRAN floating point accuracy problems
(04042016 09:02 AM)Werner Wrote: Can I chime in? Sure. 8) (04042016 09:02 AM)Werner Wrote: The following formulas, courtesy of John H. Meyers on comp.sys.hp48 do not require any fiddling at all: Ah, great. Yes, now I remember these formulas, too. Although I cannot definitely say (yet) if this solves all problems due to roundoff by nondecimal number representations, I think this should work very well. So a Fortran version of this sounds like a good idea. I would give it a try. Maybe (or even "probably") this will solve the whole issue very elegantly. Dieter 

04062016, 07:15 PM
(This post was last modified: 04062016 07:25 PM by Dieter.)
Post: #34




RE: FORTRAN floating point accuracy problems
(04042016 06:42 PM)Dieter Wrote:(04042016 09:02 AM)Werner Wrote: The following formulas, courtesy of John H. Meyers on comp.sys.hp48 do not require any fiddling at all: I now took a closer look at this, at here is what I found: 1. In most cases the code works surprisingly well. The key seems to be the division by 150 that shifts roundoff errors in the number representation by two digits to the right. Here they usually will not seriously affect the final result. 2. However, there are cases where errors may occur. I did some tests in Excel VBA. Here is an example: Consider the case dms = 1/5,000 00000 00001 = 0,19999 99999 99999 6000... which is displayed as 0,2 since Excel does not show more than 15 significant digits. Try 10*this – 1 and you see the correct value 0,99999 99999 99996. In the above code, the term 100*dms rounds to 20 even with "all" 15 digits displayed. Now INT of this yields 20 exactly (not 19!) and FP(100*dms) = 100*dms – INT(dms) returns –0,00000 00000 00039, i.e. a negative value! This messes up the whole calculation and the final result is 0,34444... hours instead of 0,33333... The key problem in the above example is the INT function that obviously returns the integer portion of what you see with 15 digits. This might avoid confusions for many users ("why is int(20) = 19?"), but numerically everything is going the wrong way from here. Well... "wrong" is a hard word. What is the conversion code supposed to do? Should it take 0,19999 99999 99999 6 as 19 minutes and 99,9999... seconds? Then the "wrong" result is correct. Or should it assume that this is supposed to be 0,2 exactly since the 16th digit can be off anytime? Then the mentioned problem occurs. So it seems the only reliable way is rounding the input value to an integer with as many digits as the working precision reliably allows. For instance 15 digits in IEEE754 double precision. Or at least round the fractional part of the input value to the maximum number of decimals that is reliably supported. So an input value like 30,2 would be rounded to 13 decimals (15 minus 2 for the integer part). Code for this approach has already been shown in this thread. Dieter 

04072016, 12:21 PM
Post: #35




RE: FORTRAN floating point accuracy problems
I have not had time to get to try the piece of code Werner posted from John Meyers and will probably not have time to get back to it for several weeks now.
I did want to post an update for Dieter on his rounding routine suggestion. It works beautifully with my new double precision truncation function. Code:
Thanks so much for the help! This is really a great tool to have, especially for FORTRAN IV where there is a very limited selection of builtins and no obvious way to do this for double precision values. I tried using this routine with our old code from earlier in the post but it did not help. We were still off badly on the same values we missed before. For some reason several of the results seem to have been rounded in the wrong direction. I tested the 1.25D0 value specifically because this was one of the values that seems to have gotten misadjusted upwards instead of rounded. It clearly works fine here. This doesn't make sense right now so I have to track it down later when I have time. I'm sorry I haven't been able to get back to the Meyer's code yet. I'll be very busy for another several weeks before I can get to it. Thanks guys. It ain't OVER 'till it's 2 PICK 

04072016, 01:36 PM
Post: #36




RE: FORTRAN floating point accuracy problems
I was waiting till you came up with a problematic case ;)
Actually, however, depending on how you look at it:  this is a case of wrong input (99 for seconds)  the result is correct (try it on Free42 Decimal, for instance) But I agree that you can't see it as it is rounded in the display.. BTW with x = 1e13/(5e13 + 1), on Free42 Binary, INT(100*x) delivers 19 and not 20 as you stated, and there are no negatives during the computation? Cheers, Werner 

04072016, 02:28 PM
Post: #37




RE: FORTRAN floating point accuracy problems
(04072016 01:36 PM)Werner Wrote:  this is a case of wrong input (99 for seconds) In decimal arithmetics everything is easy. The point here is that already the binary representation of the argument may be off (or "wrong", as you put it). Consider the case dms=30,2 or 151/5. The closest possible binary doubleprecision representation of this translates to 30,199 99999 99999 993 in decimal notation. Even if you manually enter 30,2 the program will work with 30,1999... so that all the mentioned problems occur. Rounding (up) the input to 15 significant digits solves this. (04072016 01:36 PM)Werner Wrote: BTW with x = 1e13/(5e13 + 1), on Free42 Binary, INT(100*x) delivers 19 and not 20 as you stated, That's an easy case that is handled correctly even by Excel. Replace the "E13" part by E14 or E15, depending on Free42's working precision, so that the "+1" occurs in the least significant digit. Dieter 

04112016, 06:19 AM
Post: #38




RE: FORTRAN floating point accuracy problems
Hello,
Perhaps you should try these formulas: For x>0: DEG>DMS(x)=(90*x+INT(60*x)+100*INT(x))/250 It's mathematically exact, but special care is needed in evaluating INT(60*x) on a calculator. => It seems better to round 60*x to last internal digit before applying INT. DMS>DEG(x)=(250*xINT(100*x)60*INT(x))/90 See you. 

04112016, 11:41 AM
Post: #39




RE: FORTRAN floating point accuracy problems
(04072016 02:28 PM)Dieter Wrote:(04072016 01:36 PM)Werner Wrote: BTW with x = 1e13/(5e13 + 1), on Free42 Binary, INT(100*x) delivers 19 and not 20 as you stated, But it's your example a few posts up : Quote:Consider the case dms = 1/5,000 00000 00001 of which you claimed Excel's INT function of 100*dms yielded 20: Quote:Now INT of this yields 20 exactly (not 19!) Well, never mind. If you do want to use d.ms format, which can't be represented exactly in binary floating point, you'll need rounding. Cheers, Werner 

04112016, 12:03 PM
(This post was last modified: 04112016 12:05 PM by Dieter.)
Post: #40




RE: FORTRAN floating point accuracy problems
(04112016 11:41 AM)Werner Wrote: But it's your example a few posts up : There was one zero missing. #) The example is supposed to have the 1 in the last significant (i.e. 15th) digit. So that should be 1/5,0000 00000 00001. Or 1E14/(5E14+1). That's why I said you should replace the "E13" parts by "E14". OTOH the decimal value 0,19999 99999 99999 6 was correct. In this case Excel indeed returns 20 for INT(100*dms). (04112016 11:41 AM)Werner Wrote: Well, never mind. If you do want to use d.ms format, which can't be represented exactly in binary floating point, you'll need rounding. That's what I said. And that's why I definitely prefer BCD arithmetics. BCD is the standard on pocket calculators where this d.ms notation is common, and here it works well. On the other hand most programming languages use binary arithmetics, so here a different approach should be preferred. For instance with a record, a struct or whatever it's called in your preferred language. Dieter 

« Next Oldest  Next Newest »

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