Approximation in Real(8)

[h_amini]

Subroutine D2RS(A, N, D2R) reads double precision A and rounds it to N decimal points to give D2R. In some cases there is a problem in the last digit of D2R. For instance when A = .04D0, the subroutine gives D2R = 4.000000000000001E-02. I think the reason is that the closest approximation in REAL(8) to .04D0 is 4.000000000000001E-02 as line 8 shows too, but it causes problem in my program. So could someone let me know how to solve this problem.

Hamid

[bash]      PROGRAM TEST
      IMPLICIT NONE
      DOUBLE PRECISION D2R
      INTEGER N 
      N = 7 ! The number of places after the decimal point
      
      CALL D2RS(.04D0, N, D2R); PRINT*, 'D2R =', D2R
8     PRINT*, '.4D0 * ... =', .4D0 * 10.D0 ** -1
  
      END PROGRAM TEST
      
      
      
      SUBROUTINE D2RS(A, N, D2R)
      IMPLICIT NONE
      REAL(8) D2R
      DOUBLE PRECISION A, B, C
      INTEGER E, J, N, I
      CHARACTER*23 TEMP
      
      TEMP = ''
      WRITE(TEMP, FMT = '(E23.16E2)') A
    ! Number without exponent
      READ(UNIT = TEMP, FMT = "(E.)") B 
    ! The (N+1)-th digit after decimal point
      READ(UNIT = TEMP, FMT = "(<3+N>X,I1)") J 
    ! The number of digits in the exponent
      READ(UNIT = TEMP, FMT = "(20X,I3)") E

      !PRINT*, '1B =', B
!_____________________
      IF (J .GE. 5) B = B + SIGN(10.D0 ** -N,B) ! Update last digit
      !PRINT*, '1J', J
      
      IF (J .EQ. 4) THEN
      LOOP_A : DO I = 1, 8 ! Check next digit
        READ(UNIT = TEMP, FMT = "(<3+N+I>X,I1)") J
        !PRINT*, 'J', J
        IF (J .LT. 4) EXIT
        IF (J .GT. 4) THEN
          B = B + SIGN(10.D0 ** -N,B) ! Update last digit
          EXIT LOOP_A
        END IF
      ! If J = 4 --> continue checking next numbers
      END DO LOOP_A
      END IF
!_____________________
      !PRINT*, 'TEMP =', TEMP, ', B =', B, ', E =', E
      
      D2R = B * 10.D0 ** E
      !print*,'B',B,'E',E,'D2R',D2R ! temp

      RETURN
      END SUBROUTINE D2RS[/bash]

[IDZ_A_Intel]

A binary system representing fractions (e.g. IEEE 754) cannot always represent decimal fractions as you would expect from a decimal system. There are some libraries that internally store numbers as (equivilent to) strings and perform decimal math on these.

With a decimal math system you can precisely represent 0.04

However, you cannot precisely represent 1/3 as it is an infinately repeating fraction.l

In a binary system that represents fractions the internal representation for the decimal expression 1/10 is also an infinately repeating fraction. Only fractions representable by sums of 1/powers of 2 where the number of bits difference between the largest fraction and smallest fraction lay within the number of bits in the mantissa, are exactly representable.

0.25 = 1/4 (requires 1 bit exact)

0.04 == 1/25 = 1/32 + 1/64... requires more bits than in mantissa for REAL(8) therefor approximate

If you require exact monitary values in cents, then use cents rather than dollars, (or mills, or billionths, ...)

Jim Dempsey

[h_amini]

Thank you very much for this. My numbers are dimensionless so unites cannot be converted. As I understand the only solution to represent 0.04 precisely is using some libraries. Could you let me know how to get them?

Another thing is in my program I compare the rounded 0.04 with its precise value, for example

CALLD2RS(.04D0,N,D2R) ! It gives: D2R = 4.000000000000001E-02

IF(D2R .EQ. .04D0) PRINT*, OK

Clearly the condition in the IF statement is not met, so is there any solution for comparing numbers up to a given decimal point instead of precisely?

Hamid

[Arjen_Markus]

Perhaps the module "decimal_arithmetic" in my Flibs project - http://flibs.sf.net - can help you.

Regards,

Arjen

[TimP]

Quoting h.amini

Another thing is in my program I compare the rounded 0.04 with its precise value, for example

CALLD2RS(.04D0,N,D2R) ! It gives: D2R = 4.000000000000001E-02

IF(D2R .EQ. .04D0) PRINT*, OK

Clearly the condition in the IF statement is not met, so is there any solution for comparing numbers up to a given decimal point instead of precisely?

It's not clear that the condition will fail. If both the constant and the variable take the nearest 64-bit binary representation to decimal .04, they should compare equal. As you probably know, a common solution is to compare with a relative tolerance, e.g.

if( abs(d2r - .04d0) < 2 * spacing(.04d0) ) ....

[IDZ_A_Intel]

If it is always two digits to the right of the decimal point, you can read values in using a -2P scale factor in a format, which will multiply the value by 100 and allow exact representation (assuming the value isn't too large). You can then write it out using 2P to divide by 100.

[h_amini]

Much appreciated all the comments.

@Arjen: The project seems very interesting and Ill have a look on the detail as soon as I am free.

@Tim: The constant does not take the nearest binary representation and is precisely .04D0 so the condition is not met. Finally I used comparing tolerance to solve my problem but I wondered if there might be a more straightforward solution such as a function for this purpose which seems not.

@Steve: The numbers are more general to read with the scale factor but as mentioned I checked the tolerance for this problem.