FORTRAN floating point accuracy problems

[Werner]

(04-11-2016 12:03 PM)Dieter Wrote:  In this case Excel indeed returns 20 for INT(100*dms).

Indeed it does.. but it's an example of Excel 'fiddling' - it adjusted the results to match its decimal audience ;-)

Free42 Binary returns 20 only with 1e16 - because then 5e16 + 1 equals 5e16.
Excel really should use decimal arithmetic..

Cheers, Werner

[Dieter]

(04-11-2016 12:35 PM)Werner Wrote:  

(04-11-2016 12:03 PM)Dieter Wrote:  In this case Excel indeed returns 20 for INT(100*dms).

Indeed it does.. but it's an example of Excel 'fiddling' - it adjusted the results to match its decimal audience ;-)
(...)
Excel really should use decimal arithmetic..

Well, Excel and mathematics... that's two different worlds. Especially as far as more advanced functions are concerned. I remember how older versions handled the Normal quantile function – the results were inaccurate (about 5 digits) and below a certain limit the result simply was... ±500000. *facepalm*

(04-11-2016 12:35 PM)Werner Wrote:  Free42 Binary returns 20 only with 1e16 - because then 5e16 + 1 equals 5e16.

Hmm....
Free42 Binary uses standard IEEE754 double precision and the computer's FPU, that's 16 valid digits at best.
Now 1E15/(5E15+1) = 0,19999 99999 99999 960000..., which, even when rounded to 16 digits, equals 0,2.
If IP(100*this) is 19, this of course is mathematically correct. But does Free42 Binary display 100x the above fraction as 19,999... or 20? So does it display 20 and at the same time return 19 as the integer portion of this?

Things like these are the reason why I prefer Free42 Decimal. ;-)

Dieter

[Werner]

1E15/(5E15 +1) on Free42 Bin equals (hex digits

1.9999999999998 * 2^(-3)

In the display of course, this shows as 0.2, but after

1e8 x FP

you see

9.99999996275e-1

And the '6' in its decimal representation is really the 16th digit, not the 17th (the leading '1' is the implied bit, not a full digit), it doesn't get rounded.

When you enter 0.2, it's hex representation is

1.999999999999A * 2^(-3)

Cheers, Werner

[EdS2]

I don't think binary arithmetic is to blame here - rather, it's the lack of a true conversion to integer. Either truncation or rounding would be fine, so long as they are accurate to all digits. The problem here is a little extra "cosmetic" rounding which seems to be difficult or impossible to avoid.

It might be OK if you take the integer part and then check by subtraction that it's the right value - now you have a program instead of a formula.

[HP67]

Guys I have been preparing for a big holiday coming up so I haven't been able to get back to this thread. I am still interested and will be able to update in a few weeks. Sorry for the delay and thanks for all your comments and suggestions.

[HP67]

Ok guys, I am almost back. Dieter, do you have an idea how to proceed now that we have your rounding function available?

Anybody else, ideas I should try? I'll go over the thread again since I have forgotten where we left off.

Thanks.

[Dieter]

(05-03-2016 11:54 AM)HP67 Wrote:  Ok guys, I am almost back. Dieter, do you have an idea how to proceed now that we have your rounding function available?

Using the rounding function you could try something like this:

Code:

Function dms2dd(dms)

Const safedigits = 15  ' or in your case: 16

dd = Int(dms)

If dms > 1 Then
   k = Int(log10(dd))    ' I assume Fortran features an exact log10 function
Else
   k = 0
End If

n = safedigits - 3 - k
mmss = Round(100 * (dms - dd), n)
mm = Int(mmss)
ss = (mmss - mm) * 100

dms2dd = dd + mm / 60 + ss / 3600

End Function

Here "safedigits" is the number of significant decimals we can consider exact. That's 15 for IEEE754 DP and 16 in your case.

The algorithm essentially determines how many digits after the decimal point can be assumed exact (15 resp. 16 minus the number of digits of the integer part) and rounds the mm.ss part accordingly.

Dieter

[HP67]

Thanks, I'll try this and post back when I can.

It would be nice if there were some way of determining safedigits environmentally, without having to code a constant, so the code would be somewhat portable.

[Dieter]

(05-03-2016 01:02 PM)HP67 Wrote:  It would be nice if there were some way of determining safedigits environmentally, without having to code a constant, so the code would be somewhat portable.

Some languages feature a function that returns the smalles possible machine epsilon. Others could use something like this:

Code:

Function safedigits()

e = 1.0   ' e has the data type for which we want to know the number of valid decimal digits

Do While 1.0 + e > 1.0
   e = e * 0.5
Loop

safedigits = Int(-Log10(e))

End Function

However, things may be slightly different with base-16 encoded numbers like the ones used by your hardware: cf. here.

Dieter

[HP67]

Code:


T1DMS2DD

 TEST          DMS       EXPECT       RESULT         DIFF

    1   89.1115000   89.1875000   89.1875000    0.0000000
    2   12.1500000   12.2500000   12.2500000    0.0
    3   33.3000000   33.5000000   33.5000000    0.0
    4   71.0030000   71.0083333   71.0083333    0.0000000
    5   42.2453000   42.4147222   42.4147222    0.0000000
    6   38.4225000   38.7069444   38.7069444    0.0000000
    7   29.3030000   29.5083333   29.5083333    0.0000000
    8    0.4949000    0.8302777    0.8302778    0.0000001
    9   75.1500000   75.2500000   75.2500000    0.0
   10   43.2230000   43.3750000   43.3750000    0.0000000
   11    9.3345000    9.5625000    9.5625000    0.0000000
   12   33.5752200   33.9645000   33.9645000    0.0
   13   13.0724420   13.1234500   13.1234500    0.0000000
   14   21.3000000   21.5000000   21.5000000    0.0
   15   59.4721120   59.7892000   59.7892000    0.0
   16   65.1100960   65.1836000   65.1836000    0.0000000

Outstanding! Thanks, Dieter!

I'll get back to you on how I made out with your program that calculates the epsilon later.

[Dieter]

(05-04-2016 11:05 AM)HP67 Wrote:  Outstanding! Thanks, Dieter!

Fine. There is just one issue: if the input is very small, i.e. dms is a fraction of a second, or maybe even something below 1E–16, the input may be rounded down to zero. This can be avoided by replacing the respective lines by this:

Code:

If dms >= 1 Then
   k = int(log10(dms))
Else
   k = int(log10(dms)) - 1
End If

If there was a FLOOR function the test was not required at all: k=FLOOR(log10(dms).
Maybe you can write one in Assembler. ;-)

Now add a few lines that handle dms<0 (there is a DSIGN function in Fortran IV) or dms=0 (simply return zero), and you're done.

Edit: this could be a solution.

Code:

Function dms2dd(dms)
Const safedigits = 16
adms = Abs(dms)

If adms >= 1 Then
   k = Int(log10(adms))
Else
   If adms = 0 Then
      k = 0
   Else
      k = Int(log10(adms)) - 1
   End If
End If

n = safedigits - 3 - k

dd = Int(adms)
mmss = Round(100 * (adms - dd), n)
mm = Int(mmss)
ss = (mmss - mm) * 100

dms2dd = dsign(dd + mm / 60 + ss / 3600, dms)

End Function

If I get it right, dsign transfers the sign of the second argument to the first one, i.e. dsgn(x, y) = abs(x) * sign(y), which is quite handy here.

Dieter

[HP67]

(05-04-2016 01:34 PM)Dieter Wrote:  Fine. There is just one issue: if the input is very small, i.e. dms is a fraction of a second, or maybe even something below 1E–16, the input may be rounded down to zero. This can be avoided by replacing the respective lines by this:

Code:

If dms >= 1 Then
   k = int(log10(dms))
Else
   k = int(log10(dms)) - 1
End If

If there was a FLOOR function the test was not required at all: k=FLOOR(log10(dms).
Maybe you can write one in Assembler. ;-)

Actually I had a FLOOR function in FORTRAN that was based on the less than optimal IDINT function that I replaced with my handy-dandy assembler DINT. So this is available now. This is what I have been using, does it look ok or can it be improved:

Code:


      DOUBLE PRECISION FUNCTION DFLOOR(X)                               00040003
      REAL*8 X                                                          00050000
      REAL*8 DINT                                                       00070004
C                                                                       00080004
      DFLOOR = DINT(X)                                                  00090004
      IF (X .LT. 0.0D0) DFLOOR = DFLOOR - 1.0D0                         00100002
      RETURN                                                            00110000
      END                                                               00120000

(05-04-2016 01:34 PM)Dieter Wrote:  Now add a few lines that handle dms<0 (there is a DSIGN function in Fortran IV) or dms=0 (simply return zero), and you're done.

Edit: this could be a solution.

Code:

Function dms2dd(dms)
Const safedigits = 16
adms = Abs(dms)

If adms >= 1 Then
   k = Int(log10(adms))
Else
   If adms = 0 Then
      k = 0
   Else
      k = Int(log10(adms)) - 1
   End If
End If

n = safedigits - 3 - k

dd = Int(adms)
mmss = Round(100 * (adms - dd), n)
mm = Int(mmss)
ss = (mmss - mm) * 100

dms2dd = dsign(dd + mm / 60 + ss / 3600, dms)

End Function

If I get it right, dsign transfers the sign of the second argument to the first one, i.e. dsgn(x, y) = abs(x) * sign(y), which is quite handy here.

Dieter

Given our beautiful results above is this necessary? If you think it's a good idea I'll look into that later since I have shut the machine off already.

Thanks again.

[Dieter]

(05-04-2016 02:50 PM)HP67 Wrote:  Actually I had a FLOOR function in FORTRAN that was based on the less than optimal IDINT function that I replaced with my handy-dandy assembler DINT. So this is available now. This is what I have been using, does it look ok or can it be improved:

It can be improved so that it returns correct results for integer x<0. The given implementation will return FLOOR(–3) = –4, which of course is wrong. You might test if x<0 and x≠int(x).

(05-04-2016 02:50 PM)HP67 Wrote:  Given our beautiful results above is this necessary?

The suggested changes provide correct results for arguments that are extremely small, zero or negative. I think this makes sense.

Dieter

[HP67]

Thanks. I have no idea how I could have missed that. Something looked wrong to me. Fixed!

I'll get back to the other improvements asap but have been unable to get back to this for a long time since some big things have been going on.