summaryrefslogtreecommitdiff
path: root/numpy/core/src/npymath
diff options
context:
space:
mode:
authormattip <matti.picus@gmail.com>2022-08-22 10:25:28 +0300
committerSebastian Berg <sebastian@sipsolutions.net>2022-09-27 18:11:07 +0200
commitc6366c16844b162749dca5aab5f8541b319d99c9 (patch)
tree0203bdf67d4696dd23f005fcffa20742c06d04f4 /numpy/core/src/npymath
parent9c5a2b41cc0ade160a01b7423630f7ee1b8e63dc (diff)
downloadnumpy-c6366c16844b162749dca5aab5f8541b319d99c9.tar.gz
BUILD: move some never-blocked complex math functions to 'mandatory'
Diffstat (limited to 'numpy/core/src/npymath')
-rw-r--r--numpy/core/src/npymath/npy_math_complex.c.src462
1 files changed, 18 insertions, 444 deletions
diff --git a/numpy/core/src/npymath/npy_math_complex.c.src b/numpy/core/src/npymath/npy_math_complex.c.src
index ce2772273..696495dab 100644
--- a/numpy/core/src/npymath/npy_math_complex.c.src
+++ b/numpy/core/src/npymath/npy_math_complex.c.src
@@ -158,9 +158,7 @@ npy_carg@c@(@ctype@ z)
#define SCALED_CEXP_LOWERL 11357.216553474703895L
#define SCALED_CEXP_UPPERL 22756.021937783004509L
-#if !defined(HAVE_CSINH@C@) || \
- !defined(HAVE_CCOSH@C@) || \
- !defined(HAVE_CEXP@C@)
+#if !defined(HAVE_CEXP@C@)
static
@ctype@
@@ -523,443 +521,6 @@ npy_cpow@c@ (@ctype@ a, @ctype@ b)
#endif
}
-
-#ifndef HAVE_CCOS@C@
-@ctype@
-npy_ccos@c@(@ctype@ z)
-{
- /* ccos(z) = ccosh(I * z) */
- return npy_ccosh@c@(npy_cpack@c@(-npy_cimag@c@(z), npy_creal@c@(z)));
-}
-#endif
-
-#ifndef HAVE_CSIN@C@
-@ctype@
-npy_csin@c@(@ctype@ z)
-{
- /* 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
-
-#ifndef HAVE_CTAN@C@
-@ctype@
-npy_ctan@c@(@ctype@ z)
-{
- /* ctan(z) = -I * ctanh(I * z) */
- z = npy_ctanh@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
-
-#ifndef HAVE_CCOSH@C@
-/*
- * Taken from the msun library in FreeBSD, rev 226599.
- *
- * Hyperbolic cosine of a complex argument z = x + i y.
- *
- * cosh(z) = cosh(x+iy)
- * = cosh(x) cos(y) + i sinh(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.
- *
- * CCOSH_BIG is chosen such that
- * spacing(0.5 * exp(CCOSH_BIG)) > 0.5*exp(-CCOSH_BIG)
- * although the exact value assigned to CCOSH_BIG is not so important
- */
-@ctype@
-npy_ccosh@c@(@ctype@ z)
-{
-#if @precision@ == 1
- const npy_float CCOSH_BIG = 9.0f;
- const npy_float CCOSH_HUGE = 1.70141183e+38f;
-#endif
-#if @precision@ == 2
- const npy_double CCOSH_BIG = 22.0;
- const npy_double CCOSH_HUGE = 8.9884656743115795e+307;
-#endif
-#if @precision@ >= 3
-#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
- const npy_longdouble CCOSH_BIG = 22.0L;
- const npy_longdouble CCOSH_HUGE = 8.9884656743115795e+307L;
-#else
- const npy_longdouble CCOSH_BIG = 24.0L;
- const npy_longdouble CCOSH_HUGE = 5.94865747678615882543e+4931L;
-#endif
-#endif
-
- @type@ x, y, h, absx;
- npy_int xfinite, yfinite;
-
- x = npy_creal@c@(z);
- y = npy_cimag@c@(z);
-
- 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_cosh@c@(x), x * y);
- }
- absx = npy_fabs@c@(x);
- if (absx < CCOSH_BIG) {
- /* small x: normal case */
- return npy_cpack@c@(npy_cosh@c@(x) * npy_cos@c@(y),
- npy_sinh@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@(absx) * 0.5@c@;
- return npy_cpack@c@(h * npy_cos@c@(y),
- npy_copysign@c@(h, x) * npy_sin@c@(y));
- }
- else if (absx < 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_cimag@c@(z) * npy_copysign@c@(1, x));
- }
- else {
- /* x >= 1455: the result always overflows */
- h = CCOSH_HUGE * x;
- return npy_cpack@c@(h * h * npy_cos@c@(y), h * npy_sin@c@(y));
- }
- }
-
- /*
- * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
- * The sign of 0 in the result is unspecified. Choice = normally
- * the same as dNaN. Raise the invalid floating-point exception.
- *
- * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
- * The sign of 0 in the result is unspecified. Choice = normally
- * the same as d(NaN).
- */
- if (x == 0 && !yfinite) {
- return npy_cpack@c@(y - y, npy_copysign@c@(0, x * (y - y)));
- }
-
- /*
- * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
- *
- * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
- * The sign of 0 in the result is unspecified.
- */
- if (y == 0 && !xfinite) {
- return npy_cpack@c@(x * x, npy_copysign@c@(0, x) * y);
- }
-
- /*
- * cosh(x +- I Inf) = dNaN + I dNaN.
- * Raise the invalid floating-point exception for finite nonzero x.
- *
- * cosh(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));
- }
-
- /*
- * cosh(+-Inf + I NaN) = +Inf + I d(NaN).
- *
- * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
- * The sign of Inf in the result is unspecified. Choice = always +.
- * Raise the invalid floating-point exception.
- *
- * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
- */
- if (npy_isinf(x)) {
- if (!yfinite) {
- return npy_cpack@c@(x * x, x * (y - y));
- }
- return npy_cpack@c@((x * x) * npy_cos@c@(y), x * npy_sin@c@(y));
- }
-
- /*
- * cosh(NaN + I NaN) = d(NaN) + I d(NaN).
- *
- * cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
- * Optionally raises the invalid floating-point exception.
- * Choice = raise.
- *
- * cosh(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 CCOSH_BIG
-#undef CCOSH_HUGE
-#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.
- */
-@ctype@
-npy_csinh@c@(@ctype@ z)
-{
-#if @precision@ == 1
- const npy_float CSINH_BIG = 9.0f;
- const npy_float CSINH_HUGE = 1.70141183e+38f;
-#endif
-#if @precision@ == 2
- const npy_double CSINH_BIG = 22.0;
- const npy_double CSINH_HUGE = 8.9884656743115795e+307;
-#endif
-#if @precision@ >= 3
-#if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
- const npy_longdouble CSINH_BIG = 22.0L;
- const npy_longdouble CSINH_HUGE = 8.9884656743115795e+307L;
-#else
- const npy_longdouble CSINH_BIG = 24.0L;
- const npy_longdouble CSINH_HUGE = 5.94865747678615882543e+4931L;
-#endif
-#endif
-
- @type@ x, y, h, absx;
- npy_int xfinite, yfinite;
-
- x = npy_creal@c@(z);
- y = npy_cimag@c@(z);
-
- 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@c@;
- 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@
-/*
- * Taken from the msun library in FreeBSD, rev 226600.
- *
- * Hyperbolic tangent of a complex argument z = x + i y.
- *
- * The algorithm is from:
- *
- * W. Kahan. Branch Cuts for Complex Elementary Functions or Much
- * Ado About Nothing's Sign Bit. In The State of the Art in
- * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987.
- *
- * Method:
- *
- * Let t = tan(x)
- * beta = 1/cos^2(y)
- * s = sinh(x)
- * rho = cosh(x)
- *
- * We have:
- *
- * tanh(z) = sinh(z) / cosh(z)
- *
- * sinh(x) cos(y) + i cosh(x) sin(y)
- * = ---------------------------------
- * cosh(x) cos(y) + i sinh(x) sin(y)
- *
- * cosh(x) sinh(x) / cos^2(y) + i tan(y)
- * = -------------------------------------
- * 1 + sinh^2(x) / cos^2(y)
- *
- * beta rho s + i t
- * = ----------------
- * 1 + beta s^2
- *
- * Modifications:
- *
- * I omitted the original algorithm's handling of overflow in tan(x) after
- * verifying with nearpi.c that this can't happen in IEEE single or double
- * precision. I also handle large x differently.
- */
-
-#define TANH_HUGE 22.0
-#define TANHF_HUGE 11.0F
-#define TANHL_HUGE 42.0L
-
-@ctype@
-npy_ctanh@c@(@ctype@ z)
-{
- @type@ x, y;
- @type@ t, beta, s, rho, denom;
-
- x = npy_creal@c@(z);
- y = npy_cimag@c@(z);
-
- /*
- * ctanh(NaN + i 0) = NaN + i 0
- *
- * ctanh(NaN + i y) = NaN + i NaN for y != 0
- *
- * The imaginary part has the sign of x*sin(2*y), but there's no
- * special effort to get this right.
- *
- * ctanh(+-Inf +- i Inf) = +-1 +- 0
- *
- * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite
- *
- * The imaginary part of the sign is unspecified. This special
- * case is only needed to avoid a spurious invalid exception when
- * y is infinite.
- */
- if (!npy_isfinite(x)) {
- if (npy_isnan(x)) {
- return npy_cpack@c@(x, (y == 0 ? y : x * y));
- }
- return npy_cpack@c@(npy_copysign@c@(1,x),
- npy_copysign@c@(0,
- npy_isinf(y) ?
- y : npy_sin@c@(y) * npy_cos@c@(y)));
- }
-
- /*
- * ctanh(x + i NAN) = NaN + i NaN
- * ctanh(x +- i Inf) = NaN + i NaN
- */
- if (!npy_isfinite(y)) {
- return (npy_cpack@c@(y - y, y - y));
- }
-
- /*
- * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
- * approximation sinh^2(huge) ~= exp(2*huge) / 4.
- * We use a modified formula to avoid spurious overflow.
- */
- if (npy_fabs@c@(x) >= TANH@C@_HUGE) {
- @type@ exp_mx = npy_exp@c@(-npy_fabs@c@(x));
- return npy_cpack@c@(npy_copysign@c@(1, x),
- 4 * npy_sin@c@(y) * npy_cos@c@(y) *
- exp_mx * exp_mx);
- }
-
- /* Kahan's algorithm */
- t = npy_tan@c@(y);
- beta = 1 + t * t; /* = 1 / cos^2(y) */
- s = npy_sinh@c@(x);
- rho = npy_sqrt@c@(1 + s * s); /* = cosh(x) */
- denom = 1 + beta * s * s;
- return (npy_cpack@c@((beta * rho * s) / denom, t / denom));
-}
-#undef TANH_HUGE
-#undef TANHF_HUGE
-#undef TANHL_HUGE
-#endif
-
#if !defined (HAVE_CACOS@C@) || !defined (HAVE_CASINH@C@)
/*
* Complex inverse trig functions taken from the msum library in FreeBSD
@@ -1752,10 +1313,8 @@ npy_@kind@@c@(@ctype@ z)
/**end repeat1**/
/**begin repeat1
- * #kind = cexp,clog,csqrt,ccos,csin,ctan,ccosh,csinh,ctanh,
- * cacos,casin,catan,cacosh,casinh,catanh#
- * #KIND = CEXP,CLOG,CSQRT,CCOS,CSIN,CTAN,CCOSH,CSINH,CTANH,
- * CACOS,CASIN,CATAN,CACOSH,CASINH,CATANH#
+ * #kind = cexp,clog,csqrt,cacos,casin,catan,cacosh,casinh,catanh#
+ * #KIND = CEXP,CLOG,CSQRT,CACOS,CASIN,CATAN,CACOSH,CASINH,CATANH#
*/
#ifdef HAVE_@KIND@@C@
@ctype@
@@ -1770,5 +1329,20 @@ npy_@kind@@c@(@ctype@ z)
#endif
/**end repeat1**/
+/* Mandatory complex functions */
+
+/**begin repeat1
+ * #kind = csin,csinh,ccos,ccosh,ctan,ctanh#
+ */
+@ctype@
+npy_@kind@@c@(@ctype@ z)
+{
+ __@ctype@_to_c99_cast z1;
+ __@ctype@_to_c99_cast ret;
+ z1.npy_z = z;
+ ret.c99_z = @kind@@c@(z1.c99_z);
+ return ret.npy_z;
+}
+/**end repeat1**/
/**end repeat**/