allan@alice.att.com (Allan R. Wilks) (11/12/90)
Is it reasonable to expect the cosine of (exactly) 0 to be (exactly) 1?
$ cat > a.c <<!
#include <math.h>
main()
{
printf("%g\n", 1.0-cos(0.0));
}
!
$ cc a.c -lm
$ a.out
1.11022e-16
$ cc a.c -lm43
$ a.out
0
Allan R. Wilks
AT&T Bell Labs
allan@research.att.comlmo@lsr-vax.UUCP ("Lance M. Optican - LMO") (11/14/90)
----- Begin Included Message -----
From uunet!vmb.brl.mil!info-iris-request Tue Nov 13 03:06:31 1990
Return-Path: <uunet!vmb.brl.mil!info-iris-request>
Received: from uunet.UUCP by (4.1/SMI-4.0)
id AA21136; Tue, 13 Nov 90 03:06:20 EST
Received: from VMB.BRL.MIL by uunet.uu.net (5.61/1.14) with SMTP
id AA15087; Tue, 13 Nov 90 01:37:08 -0500
Received: from vmb.brl.mil by VMB.brl.MIL id aa10955; 13 Nov 90 0:36 EST
Received: from vgr.brl.mil by VMB.BRL.MIL id aa10926; 13 Nov 90 0:23 EST
Received: from ucbvax.Berkeley.EDU by VGR.BRL.MIL id aa06765; 13 Nov 90 0:10 EST
Received: by ucbvax.Berkeley.EDU (5.63/1.42)
id AA06726; Mon, 12 Nov 90 21:10:44 -0800
Received: from USENET by ucbvax.Berkeley.EDU with netnews
for info-iris@brl.mil (info-iris@brl.mil)
(contact usenet@ucbvax.Berkeley.EDU if you have questions)
Date: 12 Nov 90 15:48:57 GMT
From: "Allan R. Wilks" <uunet!ucbvax.berkeley.edu!alice!allan>
Organization: AT&T Bell Laboratories, Murray Hill NJ
Subject: trig functions
Message-Id: <11609@alice.att.com>
Sender: uunet!BRL.MIL!info-iris-request
To: info-iris@BRL.MIL
Status: R
Is it reasonable to expect the cosine of (exactly) 0 to be (exactly) 1?
$ cat > a.c <<!
#include <math.h>
main()
{
printf("%g\n", 1.0-cos(0.0));
}
!
$ cc a.c -lm
$ a.out
1.11022e-16
$ cc a.c -lm43
$ a.out
0
Allan R. Wilks
AT&T Bell Labs
allan@research.att.com
----- End Included Message -----
It is reasonable to expect the maht library to perform well! I
tried Allan Wilks program on both Sun/3 and Sun/4 machines, and
they both gave "0" as the difference between 1.0 and cos(0.0).
Who at SGI is responsible for the math library? What standards
of compatibility are enforced with other machine architectures?
---------------------------------------------------+-----------------------
Lance M. Optican | uunet!lsr-vax!lmo
National Eye Institute, NIH, Bethesda, MD | (301) 496-3549
---------------------------------------------------+-----------------------dik@cwi.nl (Dik T. Winter) (11/14/90)
(I did already respond to the original article, but follow-up now.) In article <9011131959.AA11389@> lmo@lsr-vax.UUCP ("Lance M. Optican - LMO") writes: > ----- Begin Included Message ----- ...[lots of included message omitted]... > Is it reasonable to expect the cosine of (exactly) 0 to be (exactly) 1? > ----- End Included Message ----- > It is reasonable to expect the maht library to perform well! Yes, but that does not imply that cos(0)=1. > What standards > of compatibility are enforced with other machine architectures? What do you mean? Which machine architecture? Honeywell? (where cos(0)>1). The problem is more complex than you think. When implementing cos and sin you have two algorithmic choices (if you use additional precision you've got a few more choices): 1. Promise that the functions are monotonic at pi/4, but do not promise that cos(0)=1. This is what you get if you implement Cody-Waite's functions. You might even end up with cos(0)>1. When you reduce the range of the argument to (-pi/2,pi/2) and do an approximation there you end up with this behaviour (not exactly true, but nearly). 2. Promise that cos(0)=1, but do not promise that the functions are monotonic at pi/4 (i.e. in the neighbourhood of pi/4 it can happen that x>y but sin(x)<sin(y)). This is the way I prefer it and implement them. You get this behaviour if you reduce the range to (-pi/4,pi/4) and choose whether you need a cosine or a sine approximation. This will also make the numerical derivative at x=0 correct; but that is a moot point as it will be very bad at x=pi/4! In case 1 you can still make a choice: 3. Special case x=0. This may create non-monotonic behaviour in the neighbourhood of x=0 (it might even be that some results in the neighbourhood are larger than 1). Still, even if cos(0)!=1 it may very well be within the accuracy constraints! You should be much more apprehensive of implementations where sin(pi)=0. Throw them away, they are wrong! -- dik t. winter, cwi, amsterdam, nederland dik@cwi.nl
shenkin@cunixf.cc.columbia.edu (Peter S. Shenkin) (11/14/90)
In article <9011131959.AA11389@> lmo@lsr-vax.UUCP ("Lance M. Optican - LMO") writes: >From: "Allan R. Wilks" <uunet!ucbvax.berkeley.edu!alice!allan> >> [[ cos(0) on sgi gives 1.11022e-16 ]] > >It is reasonable to expect the maht library to perform well! I >tried Allan Wilks program on both Sun/3 and Sun/4 machines, and >they both gave "0" as the difference between 1.0 and cos(0.0). >Who at SGI is responsible for the math library? What standards >of compatibility are enforced with other machine architectures? Compatibility is not necessarily the issue, because it's possible to be compatible with a lousy standard, and machine architectures probably have less to do with it than the approximation algorithm being used. Wilks and LMO are seeing something they don't like on the Iris, and something they do like on the Sun. Before they cast stones, they should test cos(x) for every possible double-precision of x ( :-) )and make sure they still like the Sun version better. There's nothing about x=0 that is magic, and that should make it give a more accurate value of cos(x) than any other number. The classic criterion for a numerical approximation method is that it give a maximum error less than some stated value over its entire domain. Clearly, if 1.11022e-16 is less than the certified error, then this value meets the criterion. There are other issues as well, such as monotonicity, but certainly just the fact that ( cos(x) != 0. ) on SGI is not troublesome to me. -P. ************************f*u*cn*rd*ths*u*cn*gt*a*gd*jb************************** Peter S. Shenkin, Department of Chemistry, Barnard College, New York, NY 10027 (212)854-1418 shenkin@cunixc.cc.columbia.edu(Internet) shenkin@cunixc(Bitnet) ***"In scenic New York... where the third world is only a subway ride away."***
mark@mips.COM (Mark G. Johnson) (11/15/90)
In article <9011131959.AA11389@> lmo@lsr-vax.UUCP ("Lance M. Optican - LMO") writes: >> printf("%g\n", 1.0-cos(0.0)); >>$ cc a.c -lm >>$ a.out >>1.11022e-16 >>$ cc a.c -lm43 >>$ a.out >>0 > >It is reasonable to expect the maht library to perform well! I >tried Allan Wilks program on both Sun/3 and Sun/4 machines, and >they both gave "0" as the difference between 1.0 and cos(0.0). >Who at SGI is responsible for the math library? What standards >of compatibility are enforced with other machine architectures? > How about the accuracy in the vicinity of PI/4 (45 degrees)?? Here's a little program that uses trig functions to take the sine of PI/4. The correct answer is, of course, sqrt(1/2). But none of the machines (DEC-3100, Sun SPARCstation, SGI Iris) get this. Conclusion: don't be irate about noise-in-the-LSB when the argument is 0*PI; there's inaccuracy elsewhere too. ---------------------------------------------------------------------------- #include <stdio.h> #include <math.h> main() { double w, x, y, z ; w = 0.5 * acos(0.0) ; /* generate PI / 4 */ x = sin(w); /* the sine of (PI/4) is sqrt(1/2) */ y = 0.5 - (x * x) ; /* how far is it off? */ z = sqrt(0.5) - x ; /* compute the delta by another method */ printf("approx PI/4: w = %le\n", w); printf(" sine of that: x = %le\n", x); printf(" delta between x and sqrt(1/2): y = %le\n", y); printf(" delta between x and sqrt(1/2): z = %le (method 2)\n", z); } ---------------------------------------------------------------------------- SunOS Release 4.0.3c (GENERIC) #1: Thu May 25 17:17:12 PDT 1989 $ a.out approx PI/4: w = 7.853982e-01 sine of that: x = 7.071068e-01 delta between x and sqrt(1/2): y = 1.110223e-16 delta between x and sqrt(1/2): z = 1.110223e-16 (method 2) -- -- Mark Johnson MIPS Computer Systems, 930 E. Arques M/S 2-02, Sunnyvale, CA 94086 (408) 524-8308 mark@mips.com {or ...!decwrl!mips!mark}
khb@chiba.Eng.Sun.COM (Keith Bierman fpgroup) (11/16/90)
In article <43209@mips.mips.COM> mark@mips.COM (Mark G. Johnson) writes:
...
Conclusion: don't be irate about noise-in-the-LSB when the argument
is 0*PI; there's inaccuracy elsewhere too.
$ a.out
approx PI/4: w = 7.853982e-01
sine of that: x = 7.071068e-01
delta between x and sqrt(1/2): y = 1.110223e-16
delta between x and sqrt(1/2): z = 1.110223e-16 (method 2)
....
And with C1.0
chiba:/net/chiba/home2/khb>cc -fast sgi2.c && a.out
approx PI/4: w = 7.853982e-01
sine of that: x = 7.071068e-01
delta between x and sqrt(1/2): y = 1.110223e-16
delta between x and sqrt(1/2): z = 1.110223e-16 (method 2)
IEEE 754 does not require correctly rounded results for trig nor for
data conversions. Some folks are striving very hard to get "nearly
correctly rounded results" sorry for the inconvience along the way.
--
----------------------------------------------------------------
Keith H. Bierman kbierman@Eng.Sun.COM | khb@chiba.Eng.Sun.COM
SMI 2550 Garcia 12-33 | (415 336 2648)
Mountain View, CA 94043drb@eecg.toronto.edu (David R. Blythe) (11/17/90)
In article <11609@alice.att.com> allan@alice.att.com (Allan R. Wilks) writes: > >Is it reasonable to expect the cosine of (exactly) 0 to be (exactly) 1? No it is not. Its reasonable to expect a result which meets the accuracy specification for that routine. In the case of libm43 cos() has an accuracy of one unit in the last place (1 ulp) which for ieee double precision is 1.1e-16. This is documented in the intro manual pages for the math library (i.e man math). The faster cody&waite implementation of cos() in the libm library has an accuracy of 2 ulps. -drb drb@clsc.utoronto.ca
drb@eecg.toronto.edu (David R. Blythe) (11/17/90)
In article <43209@mips.mips.COM> mark@mips.COM (Mark G. Johnson) writes: > >Conclusion: don't be irate about noise-in-the-LSB when the argument >is 0*PI; there's inaccuracy elsewhere too. > >---------------------------------------------------------------------------- >#include <stdio.h> >#include <math.h> >main() >{ > double w, x, y, z ; > w = 0.5 * acos(0.0) ; /* generate PI / 4 */ > x = sin(w); /* the sine of (PI/4) is sqrt(1/2) */ > y = 0.5 - (x * x) ; /* how far is it off? */ > z = sqrt(0.5) - x ; /* compute the delta by another method */ > printf("approx PI/4: w = %le\n", w); > printf(" sine of that: x = %le\n", x); > printf(" delta between x and sqrt(1/2): y = %le\n", y); > printf(" delta between x and sqrt(1/2): z = %le (method 2)\n", z); >} >---------------------------------------------------------------------------- > >$ a.out >approx PI/4: w = 7.853982e-01 > sine of that: x = 7.071068e-01 > delta between x and sqrt(1/2): y = 1.110223e-16 > delta between x and sqrt(1/2): z = 1.110223e-16 (method 2) Trying to judge the accuracy of sin() using equally inaccurate (or worse) acos() and sqrt() routines isn't very effective or convincing. -drb drb@clsc.utoronto.ca