- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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]
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Perhaps the module "decimal_arithmetic" in my Flibs project - http://flibs.sf.net - can help you.
Regards,
Arjen
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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) ) ....
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
![](/skins/images/895D6060305DF45A57FACF854B5A8CD1/responsive_peak/images/icon_anonymous_message.png)
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page