1 files changed, 138 insertions, 0 deletions
diff --git a/libm/double/unity.c b/libm/double/unity.c
new file mode 100644
index 000000000..9223e0edf
--- /dev/null
+++ b/libm/double/unity.c
@@ -0,0 +1,138 @@
+/*							unity.c
+ *
+ * Relative error approximations for function arguments near
+ * unity.
+ *
+ *    log1p(x) = log(1+x)
+ *    expm1(x) = exp(x) - 1
+ *    cosm1(x) = cos(x) - 1
+ *
+ */
+
+#include <math.h>
+
+#ifdef ANSIPROT
+extern int isnan (double);
+extern int isfinite (double);
+extern double log ( double );
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double exp ( double );
+extern double cos ( double );
+#else
+double log(), polevl(), p1evl(), exp(), cos();
+int isnan(), isfinite();
+#endif
+extern double INFINITY;
+
+/* log1p(x) = log(1 + x)  */
+
+/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 2.32e-20
+ */
+static double LP[] = {
+ 4.5270000862445199635215E-5,
+ 4.9854102823193375972212E-1,
+ 6.5787325942061044846969E0,
+ 2.9911919328553073277375E1,
+ 6.0949667980987787057556E1,
+ 5.7112963590585538103336E1,
+ 2.0039553499201281259648E1,
+};
+static double LQ[] = {
+/* 1.0000000000000000000000E0,*/
+ 1.5062909083469192043167E1,
+ 8.3047565967967209469434E1,
+ 2.2176239823732856465394E2,
+ 3.0909872225312059774938E2,
+ 2.1642788614495947685003E2,
+ 6.0118660497603843919306E1,
+};
+
+#define SQRTH 0.70710678118654752440
+#define SQRT2 1.41421356237309504880
+
+double log1p(x)
+double x;
+{
+double z;
+
+z = 1.0 + x;
+if( (z < SQRTH) || (z > SQRT2) )
+	return( log(z) );
+z = x*x;
+z = -0.5 * z + x * ( z * polevl( x, LP, 6 ) / p1evl( x, LQ, 6 ) );
+return (x + z);
+}
+
+
+
+/* expm1(x) = exp(x) - 1  */
+
+/*  e^x =  1 + 2x P(x^2)/( Q(x^2) - P(x^2) )
+ * -0.5 <= x <= 0.5
+ */
+
+static double EP[3] = {
+ 1.2617719307481059087798E-4,
+ 3.0299440770744196129956E-2,
+ 9.9999999999999999991025E-1,
+};
+static double EQ[4] = {
+ 3.0019850513866445504159E-6,
+ 2.5244834034968410419224E-3,
+ 2.2726554820815502876593E-1,
+ 2.0000000000000000000897E0,
+};
+
+double expm1(x)
+double x;
+{
+double r, xx;
+
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+#endif
+#ifdef INFINITIES
+if( x == INFINITY )
+	return(INFINITY);
+if( x == -INFINITY )
+	return(-1.0);
+#endif
+if( (x < -0.5) || (x > 0.5) )
+	return( exp(x) - 1.0 );
+xx = x * x;
+r = x * polevl( xx, EP, 2 );
+r = r/( polevl( xx, EQ, 3 ) - r );
+return (r + r);
+}
+
+
+
+/* cosm1(x) = cos(x) - 1  */
+
+static double coscof[7] = {
+ 4.7377507964246204691685E-14,
+-1.1470284843425359765671E-11,
+ 2.0876754287081521758361E-9,
+-2.7557319214999787979814E-7,
+ 2.4801587301570552304991E-5,
+-1.3888888888888872993737E-3,
+ 4.1666666666666666609054E-2,
+};
+
+extern double PIO4;
+
+double cosm1(x)
+double x;
+{
+double xx;
+
+if( (x < -PIO4) || (x > PIO4) )
+	return( cos(x) - 1.0 );
+xx = x * x;
+xx = -0.5*xx + xx * xx * polevl( xx, coscof, 6 );
+return xx;
+}