From c117dd5fb183afb1a4790a6f6110d88704be6bf8 Mon Sep 17 00:00:00 2001
From: Eric Andersen <andersen@codepoet.org>
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 <math.h>
+
+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