1 files changed, 62 insertions, 0 deletions
diff --git a/src/base/math/sinh.c b/src/base/math/sinh.c
new file mode 100644
index 0000000..ce036ae
--- /dev/null
+++ b/src/base/math/sinh.c
@@ -0,0 +1,62 @@
+#include <u.h>
+#include <base.h>
+
+/*
+ * sinh(arg) returns the hyperbolic sine of its floating-
+ * point argument.
+ *
+ * The exponential function is called for arguments
+ * greater in magnitude than 0.5.
+ *
+ * A series is used for arguments smaller in magnitude than 0.5.
+ * The coefficients are #2029 from Hart & Cheney. (20.36D)
+ *
+ * cosh(arg) is computed from the exponential function for
+ * all arguments.
+ */
+
+static	double	p0  = -0.6307673640497716991184787251e+6;
+static	double	p1  = -0.8991272022039509355398013511e+5;
+static	double	p2  = -0.2894211355989563807284660366e+4;
+static	double	p3  = -0.2630563213397497062819489e+2;
+static	double	q0  = -0.6307673640497716991212077277e+6;
+static	double	q1   = 0.1521517378790019070696485176e+5;
+static	double	q2  = -0.173678953558233699533450911e+3;
+
+double
+math·sinh(double arg)
+{
+	double temp, argsq;
+	int sign;
+
+	sign = 0;
+	if(arg < 0) {
+		arg = -arg;
+		sign++;
+	}
+	if(arg > 21) {
+		temp = math·exp(arg)/2;
+		goto out;
+	}
+	if(arg > 0.5) {
+		temp = (math·exp(arg) - math·exp(-arg))/2;
+		goto out;
+	}
+	argsq = arg*arg;
+	temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg;
+	temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
+out:
+	if(sign)
+		temp = -temp;
+	return temp;
+}
+
+double
+math·cosh(double arg)
+{
+	if(arg < 0)
+		arg = - arg;
+	if(arg > 21)
+		return math·exp(arg)/2;
+	return (math·exp(arg) + math·exp(-arg))/2;
+}