From c117dd5fb183afb1a4790a6f6110d88704be6bf8 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 00:15:59 +0000 Subject: Seems we were lacking an acos() implementation --- libm/double/acos.c | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 libm/double/acos.c (limited to 'libm') diff --git a/libm/double/acos.c b/libm/double/acos.c new file mode 100644 index 000000000..60f61dc98 --- /dev/null +++ b/libm/double/acos.c @@ -0,0 +1,58 @@ +/* acos() + * + * Inverse circular cosine + * + * + * + * SYNOPSIS: + * + * double x, y, acos(); + * + * y = acos( x ); + * + * + * + * DESCRIPTION: + * + * Returns radian angle between 0 and pi whose cosine + * is x. + * + * Analytically, acos(x) = pi/2 - asin(x). However if |x| is + * near 1, there is cancellation error in subtracting asin(x) + * from pi/2. Hence if x < -0.5, + * + * acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) ); + * + * or if x > +0.5, + * + * acos(x) = 2.0 * asin( sqrt((1-x)/2) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -1, 1 50000 3.3e-17 8.2e-18 + * IEEE -1, 1 10^6 2.2e-16 6.5e-17 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * asin domain |x| > 1 NAN + */ + +#define __USE_BSD +#include + +double acos(double x) +{ + if (x < -0.5) { + return (M_PI - 2.0 * asin( sqrt((1+x)/2) )); + } + if (x > 0.5) { + return (2.0 * asin( sqrt((1-x)/2) )); + } + + return(M_PI_2 - asin(x)); +} -- cgit v1.2.3