Welcome to our forums
I have problems with the accuracy of trigonometric functions:
REAL*8 PI,R1,R2
R1=-1.D0 PI=DACOS(R1) R2=DSIN(PI)
The result for R2 is not 0.D0 (as I expected), it is 1.22460635822377E-16 (as shown in the debugger).
This result leads to further problems when continuing my calculations.
Detlef Pannhorst
This amount of round-off error does not seem unreasonable to me. These are 64 bit real values and round-off error will arise at some point. With Win32 FTN95 you can use 80 bit reals for greater precision but this may not make any difference in this context.
Quoted from dpannhorst I have problems with the accuracy of trigonometric functions:
REAL*8 PI,R1,R2
R1=-1.D0 PI=DACOS(R1) R2=DSIN(PI)
The result for R2 is not 0.D0 (as I expected), it is 1.22460635822377E-16 (as shown in the debugger).
This result leads to further problems when continuing my calculations.
Detlef Pannhorst
The error is comparable to machine epsilon for REAL*8. If this error is not acceptable for subsequent calculations, you need to resort to one of several changes.
You can write your own wrapper around DACOS, in which you test for special argument values such as -1d0, 0d0, etc., and when one of them is matched, return the exact value immediately instead of calling the intrinsic routine.
Otherwise, you can use higher precision, but the issues may remain, albeit with a smaller but non-zero error.
Thirdly, you can examine your calculation procedure to see why it is so sensitive to such a small error, and make the calculation more robust.
I agree with the replies above. I have always found that the most accurate way of getting PI is to use 4.0d0*ATAN(1.0d0) rather than ACOS(-1.0d0). In any case, your code should not rely on SIN(PI) returning exactly zero.
Incidently, you don't need to use the specific functions DACOS, DSIN etc, just use the generic functions ACOS and SIN.
I also use the 4xatan(1) route (appropriately adorned).
Do the generic functions produce a real*10 result too?
As well as being careful about precision, also beware of very large arguments. These come about if you are not careful in Fourier series, when computing things like:
Sin { (2m-1) (2n-1) PI } etc
The best approach is to reduce the argument back into the 0 ... 2PI space first.
Eddie
The following code example demonstrates that the resulting value of PI gave the closest value of R2 to zero. As mecej4 points out, you should check why your further calculations require a more accurate value of zero than the numerical spacing that real*8 provides.
John
REAL*8 PI,R1,R2, one, two, three, four, dx
real*10 pi_10, one_10, two_10, three_10, four_10, r_10, dx_10
integer*4 i
!
one = 1
two = 2
three = 3
four = 4
write (*,*) 'R*8'
r2 = four * atan(one)
write (*,2001) r2, spacing(r2) ! alternative version of PI
R1 = -1
PI = DACOS (R1)
write (*,2001) pi, (pi-r2) ! pi and error
R2 = DSIN (PI)
write (*,2001) pi, r2, nearest(pi,r2)-pi
!
do i = 1,5
pi = nearest (pi, r2)
r2 = dsin(pi)
write (*,2001) pi, r2
end do
2001 format (f25.20,3es12.3)
!
write (*,*) 'R*10'
one_10 = 1
two_10 = 2
three_10 = 3
four_10 = 4
r_10 = -1
r_10 = ACOS (R_10)
write (*,2001) r_10, spacing(r_10)
pi_10 = four_10 * atan(one_10)
write (*,2001) pi_10, (pi_10-r_10)
r_10 = sin(pi_10)
write (*,2001) pi_10, r_10, nearest(pi_10,r_10)-pi_10
!
r_10 = sin(pi_10)
do i = 1,5
pi_10 = nearest (pi_10, r_10)
r_10 = sin(pi_10)
write (*,2001) pi_10, r_10
end do
end
Shall we take bets that the real reason the E-16 value doesn't suit, is a test for equality to zero?
My rules (for myself) are (a) never test the result of a calculation = zero, which is where warning 179 is excellent, but (b) always test an input value* that I am going on to use as a denominator before I use it, which is where warning 179 irritates the heck out of me!
Eddie
*or expression involving input values. You can get exactly zero if you input it.
Note to John: I assume that the real*8 value for 1, 2 and 4 is exact, but for 3 contains a tiny 'error' due to the characteristic-mantissa form being held in binary. Hence, in your code, does one+two = three?
Quoted from LitusSaxonicum ... Note to John: I assume that the real*8 value for 1, 2 and 4 is exact, but for 3 contains a tiny 'error' due to the characteristic-mantissa form being held in binary. Hence, in your code, does one+two = three?
Correct in the spirit of the note, but not literally so.
The integer 3 is exactly representable in IEEE floating format in 32, 64 and 80 bits. Numbers of the form 1.xxxx..._10 which when multiplied by a power of 2 yield an integer are also exactly represented if the product of the multiplication is representable in the number of bits available in the mantissa field.
Eddie,
1,2,3,4 are all exact binary representations, as are most integers that are converted. Also 0.5 and 0.25 should be, but most notably 0.1 and 0.2 are not. So, writing :- one = 1 pi = 4 * atan (one) should provide the maximum accuracy defined by the precision of 'one'.
When writing the above example, it took me a while to remember 'spacing', whose value is the key for identifying the accuracy available. The accuracy to test is 'spacing (pi)', and not 'spacing(sin(pi))' I was impressed that the value of pi retrieved, still gave the best answer for sin(pi), to the 'nearest' value of pi. I probably should have included the run trace, which identifies that the error is smaller than spacing (pi).
I also use a real value (parameter) to identify zero or some other unset value, such as
real*8, parameter :: dummy = -0.98765 array = dummy ! then test: if (array(i,j)==dummy) ... ! safely identified as not being updated
this always gives a warning that can be ignored
John
I suppose that I'm astonished that 3 and larger odd numbers, especially primes, are ever held exactly in a base 2 characteristic-mantissa form, and I program things with a healthy degree of scepticism about rounoff. I'm not surprised that 2 to the power n, where n is a positive or negative integer, is held exactly.
Eddie
Quoted from LitusSaxonicum I suppose that I'm astonished that 3 and larger odd numbers, especially primes, are ever held exactly in a base 2 characteristic-mantissa form, and I program things with a healthy degree of scepticism about rounoff. I'm not surprised that 2 to the power n, where n is a positive or negative integer, is held exactly.
Eddie
Hmmm, 3 = 21 + 20 😉
Eddie,
For pi = 4 * atan (one) I still use a parameter four, rather than 4, even though the standard implies the precision will be the same. chasing round-off errors can be difficult to identify (write), although NEAREST and SPACING are useful functions to quantify what to expect.
Regarding if (a==b) for reals, all stored reals have an excact value. It's just they might not be what we expect if they are derived from a calculation. 1.0 + 2.0 will always be exactly 3.0, but 0.1 + 0.2 will probably not be 0.3 ( I should try this for R*4, 8 & 10)
John
If you run this...
program test
real*8 x,y
real*10 r,s
write(6,'(A)')'3.1415926535897932384626' ! correct value
write(6,'(F24.22)')4.D0*datan(1.D0)
x = 4.
y = 1.
write(6,'(F24.22)')x*atan(y)
r = 4.
s = 1.
write(6,'(F24.22)')r*atan(s)
end
... you can see the improvement when going from real8 to real10. On the other hand, there is no difference using 'atan' or 'datan', neither writing '4.' or '4.D0'.
Regards - Wilfried
Wilfred,
I added to your example to show other alternatives for implied precision in a calculation. I also calculated the error, both using the atan calculation and then using a temporary variable of defined precision. If you note the reporting of the error term 'pi-a*atan(b)' it appears that atan(b) is calculated and retained as 80 bit precision, even though 32 bit is all that is implied. { I only note this and do not suggest that the compiler should be changed to produce 32bit precision for atan(b). } The code is: program test
real*4 a,b,c
real*8 x,y,z
real*10 r,s,t, pi
!
pi = 3.1415926535897932384626_3
write (6,'(A)') '3.1415926535897932384626' ! correct value
write (6,22) pi, 'pi_r*10'
22 format (F24.22,2x,a,2es12.2)
write (6,22) 4.D0*datan(1.D0), '4.D0*datan(1.D0)', pi-4.D0*datan(1.D0)
a = 4.
b = 1.
x = 4.
y = 1.
r = 4.
s = 1.
c = a*atan(b) ; write (6,22) a*atan(b), 'a*atan(b)', pi-a*atan(b), pi-c
z = x*atan(y) ; write (6,22) x*atan(y), 'x*atan(y)', pi-x*atan(y), pi-z
t = r*atan(s) ; write (6,22) r*atan(s), 'r*atan(s)', pi-r*atan(s), pi-t
!
c = 4*atan(b) ; write (6,22) 4*atan(b), '4*atan(b)', pi-4*atan(b), pi-c
z = 4*atan(y) ; write (6,22) 4*atan(y), '4*atan(y)', pi-4*atan(y), pi-z
t = 4*atan(s) ; write (6,22) 4*atan(s), '4*atan(s)', pi-4*atan(s), pi-t
!
end
Compilation was ftn95 pi.f95 /lgo The results are: [FTN95/Win32 Ver. 6.30.0 Copyright (c) Silverfrost Ltd 1993-2012] NO ERRORS [<TEST> FTN95/Win32 v6.30.0] Creating executable: c:\temp\lgotemp@.exe
Program entered
3.1415926535897932384626
3.1415926535897932386000 pi_r*10
3.1415926535897931160000 4.D0*datan(1.D0) 1.23E-16
3.1415927410125732423000 a*atan(b) 0.00E+00 -8.74E-08
3.1415926535897931160000 x*atan(y) 0.00E+00 1.23E-16
3.1415926535897932386000 r*atan(s) 0.00E+00 0.00E+00
3.1415927410125732423000 4*atan(b) 0.00E+00 -8.74E-08
3.1415926535897931160000 4*atan(y) 0.00E+00 1.23E-16
3.1415926535897932386000 4*atan(s) 0.00E+00 0.00E+00
John