forums.silverfrost.com

dpannhorst · Joined: 29 Aug 2005 Posts: 165 Location: Berlin, Germany

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

PaulLaidler · Posted: Fri Nov 02, 2012 9:37 am Post subject:

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.

mecej4 · Joined: 31 Oct 2006 Posts: 1899

davidb · Joined: 17 Jul 2009 Posts: 560 Location: UK

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.
_________________
Programmer in: Fortran 77/95/2003/2008, C, C++ (& OpenMP), java, Python, Perl

LitusSaxonicum · Posted: Sat Nov 03, 2012 9:00 pm Post subject:

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

JohnCampbell · Joined: 16 Feb 2006 Posts: 2615 Location: Sydney

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

LitusSaxonicum · Posted: Sun Nov 04, 2012 1:40 pm Post subject:

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?

mecej4 · Joined: 31 Oct 2006 Posts: 1899

JohnCampbell · Joined: 16 Feb 2006 Posts: 2615 Location: Sydney

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

LitusSaxonicum · Posted: Mon Nov 05, 2012 1:15 pm Post subject:

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

davidb · Joined: 17 Jul 2009 Posts: 560 Location: UK

JohnCampbell · Joined: 16 Feb 2006 Posts: 2615 Location: Sydney

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

Wilfried Linder · Posted: Wed Nov 07, 2012 9:46 am Post subject:

If you run this...

JohnCampbell · Joined: 16 Feb 2006 Posts: 2615 Location: Sydney

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: