ENH: Import the csinh/csin implementation from FreeBSD

The code from FreeBSD was lightly adapted to fit with the numpy style.

With this commit, npy_csinh(f) and npy_csin(f) pass all of the tests in
test_c99complex.c.

numpy/core/src/npymath/npy_math_complex.c.src

@@ -482,10 +482,9 @@ static @ctype@ _npy_scaled_cexp@c@(@type@ x, @type@ y, npy_int expt)
 #ifndef HAVE_CSIN@C@
 @ctype@ npy_csin@c@(@ctype@ z)
 {
-    @type@ x, y;
-    x = npy_creal@c@(z);
-    y = npy_cimag@c@(z);
-    return npy_cpack@c@(npy_sin@c@(x) * npy_cosh@c@(y), npy_cos@c@(x) * npy_sinh@c@(y));
+    /* csin(z) = -I * csinh(I * z) */
+    z = npy_csinh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z)));
+    return npy_cpack@c@(npy_cimag@c@(z), -npy_creal@c@(z));
 }
 #endif
 
@@ -627,13 +626,132 @@ static @ctype@ _npy_scaled_cexp@c@(@type@ x, @type@ y, npy_int expt)
 #endif
 
 #ifndef HAVE_CSINH@C@
+/*
+ * Taken from the msun library in FreeBSD, rev 226599.
+ *
+ * Hyperbolic sine of a complex argument z = x + i y.
+ *
+ * sinh(z) = sinh(x+iy)
+ *         = sinh(x) cos(y) + i cosh(x) sin(y).
+ *
+ * Exceptional values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#if @precision@ == 1
+#define CSINH_BIG 9.0f
+#define CSINH_HUGE 1.70141183e+38f
+#endif
+#if @precision@ == 2
+#define CSINH_BIG 22.0
+#define CSINH_HUGE 8.9884656743115795e+307
+#endif
+#if @precision@ >= 3
+#define CSINH_BIG 24.0L
+#define CSINH_HUGE 5.94865747678615882543e+4931L
+#endif
+
 @ctype@ npy_csinh@c@(@ctype@ z)
 {
-    @type@ x, y;
+    @type@ x, y, h, absx;
+    npy_int xfinite, yfinite;
+
     x = npy_creal@c@(z);
     y = npy_cimag@c@(z);
-    return npy_cpack@c@(npy_cos@c@(y)*npy_sinh@c@(x), npy_sin@c@(y)*npy_cosh@c@(x));
+
+    xfinite = npy_isfinite(x);
+    yfinite = npy_isfinite(y);
+
+    /* Handle the nearly-non-exceptional cases where x and y are finite. */
+    if (xfinite && yfinite) {
+        if (y == 0)
+            return npy_cpack@c@(npy_sinh@c@(x), y);
+        absx = npy_fabs@c@(x);
+        if (absx < CSINH_BIG) /* small x: normal case */
+            return npy_cpack@c@(npy_sinh@c@(x) * npy_cos@c@(y), npy_cosh@c@(x) * npy_sin@c@(y));
+
+        /* |x| >= 22, so cosh(x) ~= exp(|x|) */
+        if (absx < SCALED_CEXP_LOWER@C@) {
+            /* x < 710: exp(|x|) won't overflow */
+            h = npy_exp@c@(npy_fabs@c@(x)) * 0.5;
+            return npy_cpack@c@(npy_copysign@c@(h, x) * npy_cos@c@(y), h * npy_sin@c@(y));
+        } else if (x < SCALED_CEXP_UPPER@C@) {
+            /* x < 1455: scale to avoid overflow */
+            z = _npy_scaled_cexp@c@(absx, y, -1);
+            return npy_cpack@c@(npy_creal@c@(z) * npy_copysign@c@(1, x), npy_cimag@c@(z));
+        } else {
+            /* x >= 1455: the result always overflows */
+            h = CSINH_HUGE * x;
+            return npy_cpack@c@(h * npy_cos@c@(y), h * h * npy_sin@c@(y));
+        }
+    }
+
+    /*
+     * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
+     * The sign of 0 in the result is unspecified.  Choice = normally
+     * the same as dNaN.  Raise the invalid floating-point exception.
+     *
+     * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
+     * The sign of 0 in the result is unspecified.  Choice = normally
+     * the same as d(NaN).
+     */
+    if (x == 0 && !yfinite)
+        return npy_cpack@c@(npy_copysign@c@(0, x * (y - y)), y - y);
+
+    /*
+     * sinh(+-Inf +- I 0) = +-Inf + I +-0.
+     *
+     * sinh(NaN +- I 0)   = d(NaN) + I +-0.
+     */
+    if (y == 0 && !xfinite) {
+        if (npy_isnan(x))
+            return z;
+        return npy_cpack@c@(x, npy_copysign@c@(0, y));
+    }
+
+    /*
+     * sinh(x +- I Inf) = dNaN + I dNaN.
+     * Raise the invalid floating-point exception for finite nonzero x.
+     *
+     * sinh(x + I NaN) = d(NaN) + I d(NaN).
+     * Optionally raises the invalid floating-point exception for finite
+     * nonzero x.  Choice = don't raise (except for signaling NaNs).
+     */
+    if (xfinite && !yfinite)
+        return npy_cpack@c@(y - y, x * (y - y));
+
+    /*
+     * sinh(+-Inf + I NaN)  = +-Inf + I d(NaN).
+     * The sign of Inf in the result is unspecified.  Choice = normally
+     * the same as d(NaN).
+     *
+     * sinh(+-Inf +- I Inf) = +Inf + I dNaN.
+     * The sign of Inf in the result is unspecified.  Choice = always +.
+     * Raise the invalid floating-point exception.
+     *
+     * sinh(+-Inf + I y)   = +-Inf cos(y) + I Inf sin(y)
+     */
+    if (!xfinite && !npy_isnan(x)) {
+        if (!yfinite)
+            return npy_cpack@c@(x * x, x * (y - y));
+        return npy_cpack@c@(x * npy_cos@c@(y), NPY_INFINITY@C@ * npy_sin@c@(y));
+    }
+
+    /*
+     * sinh(NaN + I NaN)  = d(NaN) + I d(NaN).
+     *
+     * sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
+     * Optionally raises the invalid floating-point exception.
+     * Choice = raise.
+     *
+     * sinh(NaN + I y)    = d(NaN) + I d(NaN).
+     * Optionally raises the invalid floating-point exception for finite
+     * nonzero y.  Choice = don't raise (except for signaling NaNs).
+     */
+    return npy_cpack@c@((x * x) * (y - y), (x + x) * (y - y));
 }
+#undef CSINH_BIG
+#undef CSINH_HUGE
 #endif
 
 #ifndef HAVE_CTANH@C@