Floating Point Consistency

[Kantor__Baruch]

Hi

I'm upgrading a CVF application that relies heavily on single precision floating point calculations.

I'm using the "Enable Floating Point Consistency" and "Extend Precision of Single Precision Constants".

Currently I create IA-32 Fortran DLL on XE.

From comparing the results of the CVF and XE versions I see that the CVF application is more precise.

It looks like in CVF there are more significant digits (or half digits).

I've read previous posts regarding Floating Point Consistency, including referred documents and presentations.

I've tried all possible combinations of the Fortran Floating Point Settings and still can not recreate precisely CVF results.

Will be grateful for any insight on this.

Baruch

[jimdempseyatthecove]

>>Successes so far have been in adding double-precision changes here and there in the code.

Make sure to pay attention to the literals that will contain precision errors. Example:

DPvariable = 0.1

DPresult = DPvariable * 0.1

0.1 cannot be expressed exactly in binary as it is a repeating fraction. The error in the SP value will get promoted. IOW an lsb error in the SP variable is in the 24'th bit of the mantissa, equating to DP results in an error in the 24th bit of a 52-bit DP mantissa.

This holds for other inexact fractions.

Use 0.1_8, or 0.1_DP (where you declare DP appropriately as a parameter).

This issue will show up in expressions, parameter declarations, and DATA initializations.

Jim Dempsey

[jimdempseyatthecove]

or

DPvariable = 0.1D0

Jim

[wtstephens]

Still stuck on the equations: 

x = EXP(y)

result = (1.0-x)/(1.0+x)

Old CVF "Update A" and gfortran get different answers using the very same REAL*8 values -- same out to 16 digits.

I am going to try to use REAL*16 for "x", and then try to round it to 10 bytes, 12 bytes, etc to see if I can mimic whatever CVF is doing.

Anyone know how to do that rounding???

[jimdempseyatthecove]

>>I am going to try to use REAL*16 for "x", and then try to round it to 10 bytes, 12 bytes, etc to see if I can mimic whatever CVF is doing.

You may get closer results, but do not expect exact results.

What you are trying to do is to match inexact results from CVF with inexact results from REAL*16 using ifx.

>> x = EXP(y)

program Console19
    implicit none
    real(8) :: X1, X2, Y
    real(16) :: temp
    real(8), parameter :: bump = 0.1_8
    integer :: i
    Y = 0.0_8
    do i = 1,1000
        Y = bump * i
        X1 = EXP(Y)
        temp = Y
        temp = EXP(temp)
        X2 = temp
        if(X1 /= X2) then
            print *, "different ", Y, X1-X2
        endif
    end do
end program Console19
------------
 different    6.60000000000000       1.136868377216160E-013
 different    11.8000000000000       2.910383045673370E-011
 different    12.6000000000000       5.820766091346741E-011
 different    84.7000000000000       1.180591620717411E+021

Jim Dempsey

[wtstephens]

The best match that I am seeing at this point comes from skipping EXP() -- and doing the exponentiation of "e" (defined out to 40 digits) myself!

  ! ---------------------------------
  REAL*16 FUNCTION EXP_MIMIC(Y) RESULT(EXPY)
    IMPLICIT NONE
    REAL*8, INTENT(IN) :: Y

    ! EXPY = EXP(Y)
    ! EXPY = TRUNC16(EXPY,60) ! REMOVE THE LOWEST N BITS

    REAL*16, PARAMETER :: E1 = 2.7182818284590452353602874713526624977572_16
    EXPY = E1**Y
    RETURN
  ENDFUNCTION

 Then performing the (1-x)/(1+x) calculation with x as REAL*16, and the final result stored in REAL*8 -- which is virtually always just 1 bit difference!

Interestingly, there is an accuracy bump I got by truncating EXP at 60 bits (out of the 112) -- but not as accurate as doing the REAL*16 exponentiation. The CVF calculation just seems not completely reproducible.

Every day I wake up with another idea.