numpy  2.0.0
src/npymath/npy_math_complex.c.src File Reference
#include "npy_math_common.h"
#include "npy_math_private.h"
#include <numpy/utils.h>

Defines

#define raise_inexact()   do { volatile npy_float junk = 1 + tiny; } while(0)
#define SCALED_CEXP_LOWERF   88.722839f
#define SCALED_CEXP_UPPERF   192.69492f
#define SCALED_CEXP_LOWER   710.47586007394386
#define SCALED_CEXP_UPPER   1454.9159319953251
#define SCALED_CEXP_LOWERL   11357.216553474703895L
#define SCALED_CEXP_UPPERL   22756.021937783004509L
#define THRESH   (@TMAX@ / (1 + NPY_SQRT2@c@))
#define TANH_HUGE   22.0
#define TANHF_HUGE   11.0F
#define TANHL_HUGE   42.0L

Functions

static NPY_INLINE ctype cadd c (@ctype @a,@ctype @b)
static NPY_INLINE ctype cneg c (@ctype @a)
static ctype _npy_scaled_cexp c (@type @x,@type @y, npy_int expt)
static NPY_INLINE type _f c (@type @a,@type @b,@type @hypot_a_b)
static NPY_INLINE void
_do_hard_work 
c (@type @x,@type @y,@type @*rx, npy_int *B_is_usable,@type @*B,@type @*sqrt_A2my2,@type @*new_y)
static NPY_INLINE void
_clog_for_large_values 
c (@type @x,@type @y,@type @*rr,@type @*ri)
static NPY_INLINE type _sum_squares c (@type @x,@type @y)

Variables

static __COMP_NPY_UNUSED npy_float tiny = 3.9443045e-31f
static const ctype c_1 c = {1.0@C@, 0.0}

Define Documentation

#define raise_inexact ( )    do { volatile npy_float junk = 1 + tiny; } while(0)
#define SCALED_CEXP_LOWER   710.47586007394386
#define SCALED_CEXP_LOWERF   88.722839f
cexp and (ccos, csin)h functions need to calculate exp scaled by another number. This can be difficult if exp(x) overflows. By doing this way, we don't risk overflowing exp. This likely raises floating-point exceptions, if we decide that we care.
This is only useful over a limited range, (see below) an expects that the input values are in this range.
This is based on the technique used in FreeBSD's __frexp_exp and __ldexp_(c)exp functions by David Schultz.
SCALED_CEXP_LOWER = log(FLT_MAX) SCALED_CEXP_UPPER = log(2) + log(FLT_MAX) - log(FLT_TRUE_MIN), where FLT_TRUE_MIN is the smallest possible subnormal number.
#define SCALED_CEXP_LOWERL   11357.216553474703895L
#define SCALED_CEXP_UPPER   1454.9159319953251
#define SCALED_CEXP_UPPERF   192.69492f
#define SCALED_CEXP_UPPERL   22756.021937783004509L
#define TANH_HUGE   22.0
Taken from the msun library in FreeBSD, rev 226600.
Hyperbolic tangent of a complex argument z = x + i y.
The algorithm is from: <blockquote> 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.</blockquote>
Method: <blockquote>

Let t = tan(x)
beta = 1/cos^2(y) s = sinh(x) rho = cosh(x)
We have:
tanh(z) = sinh(z) / cosh(z) <blockquote> <blockquote> sinh(x) cos(y) + i cosh(x) sin(y)</blockquote>

System Message: WARNING/2 (<string>, line 23) Block quote ends without a blank line; unexpected unindent.
= ---------------------------------

cosh(x) cos(y) + i sinh(x) sin(y)

cosh(x) sinh(x) / cos^2(y) + i tan(y)

= -------------------------------------

<blockquote class="first"> 1 + sinh^2(x) / cos^2(y)</blockquote>

beta rho s + i t

= ----------------
1 + beta s^2

</blockquote> </blockquote>

Modifications: <blockquote> 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.</blockquote>
#define TANHF_HUGE   11.0F
#define TANHL_HUGE   42.0L
#define THRESH   (@TMAX@ / (1 + NPY_SQRT2@c@))
We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)).

Function Documentation

static NPY_INLINE ctype cadd c ( @ctype @  a,
@ctype @  b 
) [static]
System Message: SEVERE/4 (<string>, line 1)
Missing matching underline for section title overline.

==========================================================
 Helper functions
 

<blockquote> These are necessary because we do not count on using a C99 compiler.</blockquote>

System Message: WARNING/2 (<string>, line 6) Block quote ends without a blank line; unexpected unindent.

System Message: ERROR/3 (<string>, line 6) Document may not end with a transition.
Always use this function because of the multiplication for small integer powers, but in the body use cpow if it is available.
private function for use in npy_pow{f, ,l}
divide by zeros should yield a complex inf or nan
NB: there are four complex zeros; c0 = (+-0, +-0), so that unlike for reals, c0**p, with <cite>p</cite> negative is in general ill-defined. <blockquote> c0**z with z complex is also ill-defined.</blockquote>
Raise invalid
unroll: handle inf better
unroll: handle inf better
unroll: handle inf better

References c, C, npy_cimag(), npy_cpack(), npy_creal(), and npy_fabs().

static NPY_INLINE ctype cneg c ( @ctype @  z) [static]
System Message: SEVERE/4 (<string>, line 1)
Title overline & underline mismatch.

==========================================================
 Custom implementation of missing complex C99 functions
=========================================================
algorithm from cpython, rev. d86f5686cef9
The usual formula for the real part is log(hypot(z.real, z.imag)). There are four situations where this formula is potentially problematic:
(1) the absolute value of z is subnormal. Then hypot is subnormal, so has fewer than the usual number of bits of accuracy, hence may have large relative error. This then gives a large absolute error in the log. This can be solved by rescaling z by a suitable power of 2.
(2) the absolute value of z is greater than DBL_MAX (e.g. when both z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX) Again, rescaling solves this.
(3) the absolute value of z is close to 1. In this case it's difficult to achieve good accuracy, at least in part because a change of 1ulp in the real or imaginary part of z can result in a change of billions of ulps in the correctly rounded answer.
(4) z = 0. The simplest thing to do here is to call the floating-point log with an argument of 0, and let its behaviour (returning -infinity, signaling a floating-point exception, setting errno, or whatever) determine that of c_log. So the usual formula is fine here.
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
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.
casinh(z) = z + O(z^3) as z -> 0
casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity The above formula works for the imaginary part as well, because Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)

System Message: ERROR/3 (<string>, line 6) Unexpected indentation.

<blockquote> as z -> infinity, uniformly in y</blockquote>

real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y). Assumes x and y are not NaN, and one of x and y is larger than RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use the code creal(1/z), because the imaginary part may produce an unwanted underflow. This is only called in a context where inexact is always raised before the call, so no effort is made to avoid or force inexact.
r is nan
r is +- inf
x = +inf, y = +-inf | nan
x = -inf, y = nan | +i inf
catch cases where hypot(ax, ay) is subnormal
log(+/-0 +/- 0i)
raise divide-by-zero floating point exception

<

max(ax, ay)

<

min(ax, ay)
Handle special cases.

<

raise invalid if b is not a NaN

<

return NaN + NaN i
csqrt(inf + NaN i) = inf + NaN i csqrt(inf + y i) = inf + 0 i csqrt(-inf + NaN i) = NaN +- inf i csqrt(-inf + y i) = 0 + inf i
The remaining special case (b is NaN) is handled just fine by the normal code path below.
Scale to avoid overflow.
Algorithm 312, CACM vol 10, Oct 1967.
Rescale.
ccos(z) = ccosh(I * z)
csin(z) = -I * csinh(I * z)
ctan(z) = -I * ctanh(I * z)
Handle the nearly-non-exceptional cases where x and y are finite.
small x: normal case
|x| >= 22, so cosh(x) ~= exp(|x|)

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "x".
System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "x".
x < 710: exp(|x|) won't overflow

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "x".
x < 1455: scale to avoid overflow
x >= 1455: the result always overflows
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).
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.
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).
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)
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).
Handle the nearly-non-exceptional cases where x and y are finite.
small x: normal case
|x| >= 22, so cosh(x) ~= exp(|x|)

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "x".
System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "x".
x < 710: exp(|x|) won't overflow

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "x".
x < 1455: scale to avoid overflow
x >= 1455: the result always overflows
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).
sinh(+-Inf +- I 0) = +-Inf + I +-0.
sinh(NaN +- I 0) = d(NaN) + I +-0.
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).
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)
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).
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.
ctanh(x + i NAN) = NaN + i NaN ctanh(x +- i Inf) = NaN + i NaN
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.
Kahan's algorithm

<

= 1 / cos^2(y)

<

= cosh(x)
cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf
cacos(NaN + I*+-Inf) = NaN + I*-+Inf
cacos(0 + I*NaN) = PI/2 + I*NaN with inexact
All other cases involving NaN return NaN + I*NaN. C99 leaves it optional whether to raise invalid if one of the arguments is not NaN, so we opt not to raise it.
clog...() will raise inexact unless x or y is infinite.
Avoid spuriously raising inexact for z = 1.
All remaining cases are inexact.
casin(z) = I * conj( casinh(I * conj(z)) )
catan(z) = I * conj( catanh(I * conj(z)) )
cacosh(z) = I*cacos(z) or -I*cacos(z) where the sign is chosen so Re(cacosh(z)) >= 0.
cacosh(NaN + I*NaN) = NaN + I*NaN
cacosh(NaN + I*+-Inf) = +Inf + I*NaN
cacosh(+-Inf + I*NaN) = +Inf + I*NaN
cacosh(0 + I*NaN) = NaN + I*NaN
casinh(+-Inf + I*NaN) = +-Inf + I*NaN
casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN
casinh(NaN + I*0) = NaN + I*0
All other cases involving NaN return NaN + I*NaN. C99 leaves it optional whether to raise invalid if one of the arguments is not NaN, so we opt not to raise it.
clog...() will raise inexact unless x or y is infinite.
Avoid spuriously raising inexact for z = 0.
All remaining cases are inexact.
This helps handle many cases.
To ensure the same accuracy as atan(), and to filter out z = 0.
catanh(+-Inf + I*NaN) = +-0 + I*NaN
catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2
All other cases involving NaN return NaN + I*NaN. C99 leaves it optional whether to raise invalid if one of the arguments is not NaN, so we opt not to raise it.
z = 0 was filtered out above. All other cases must raise inexact, but this is the only only that needs to do it explicitly.

References c, npy_cos(), npy_cpack(), npy_exp(), npy_frexp(), npy_ldexp(), NPY_LOGE2, and npy_sin().

static ctype _npy_scaled_cexp c ( @type @  x,
@type @  y,
npy_int  expt 
) [static]

References c, and npy_cpack().

static NPY_INLINE type _f c ( @type @  a,
@type @  b,
@type @  hypot_a_b 
) [static]
Complex inverse trig functions taken from the msum library in FreeBSD revision 251404
The algorithm is very close to that in "Implementing the complex arcsine and arccosine functions using exception handling" by T. E. Hull, Thomas F. Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, http://dl.acm.org/citation.cfm?id=275324.
Throughout we use the convention z = x + I*y.
casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B) where A = (|z+I| + |z-I|) / 2 B = (|z+I| - |z-I|) / 2 = y/A

These formulas become numerically unstable:
  1. for Re(casinh(z)) when z is close to the line segment [-I, I] (that is, Re(casinh(z)) is close to 0);
  2. for Im(casinh(z)) when z is close to either of the intervals [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is close to PI/2).
These numerical problems are overcome by defining f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2 Then if A < A_crossover, we use

System Message: ERROR/3 (<string>, line 27) Unexpected indentation.

<blockquote> log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1))) A-1 = f(x, 1+y) + f(x, 1-y)</blockquote>

System Message: WARNING/2 (<string>, line 29) Block quote ends without a blank line; unexpected unindent.
and if B > B_crossover, we use
asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y))) A-y = f(x, y+1) + f(x, y-1)
System Message: WARNING/2 (<string>, line 32) Definition list ends without a blank line; unexpected unindent.
where without loss of generality we have assumed that x and y are non-negative.
Much of the difficulty comes because the intermediate computations may produce overflows or underflows. This is dealt with in the paper by Hull et al by using exception handling. We do this by detecting when computations risk underflow or overflow. The hardest part is handling the underflows when computing f(a, b).
Note that the function f(a, b) does not appear explicitly in the paper by Hull et al, but the idea may be found on pages 308 and 309. Introducing the function f(a, b) allows us to concentrate many of the clever tricks in this paper into one function.

Docutils System Messages

System Message: ERROR/3 (<string>, line 12); backlink Undefined substitution referenced: "z+I".
System Message: ERROR/3 (<string>, line 12); backlink Undefined substitution referenced: "z-I".
System Message: ERROR/3 (<string>, line 12); backlink Undefined substitution referenced: "z+I".
System Message: ERROR/3 (<string>, line 12); backlink Undefined substitution referenced: "z-I".
System Message: ERROR/3 (<string>, line 20); backlink Undefined substitution referenced: "Im(casinh(z))".
Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2. Pass hypot(a, b) as the third argument.
static NPY_INLINE void _do_hard_work c ( @type @  x,
@type @  y,
@type @*  rx,
npy_int B_is_usable,
@type @*  B,
@type @*  sqrt_A2my2,
@type @*  new_y 
) [static]
All the hard work is contained in this function. x and y are assumed positive or zero, and less than RECIP_EPSILON. Upon return: rx = Re(casinh(z)) = -Im(cacos(y + I*x)). B_is_usable is set to 1 if the value of B is usable. If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y. If returning sqrt_A2my2 has potential to result in an underflow, it is rescaled, and new_y is similarly rescaled.

<

A, B, R, and S are as in Hull et al.

<

A-1, A-y.

<

|z+I|

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z+I".

<

|z-I|

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z-I".
A = (|z+I| + |z-I|) / 2

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z+I".
System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z-I".
Mathematically A >= 1. There is a small chance that this will not be so because of rounding errors. So we will make certain it is so.
Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y). rx = log1p(Am1 + sqrt(Am1*(A+1)))
fp is of order x^2, and fm = x/2. A = 1 (inexactly).
Underflow will not occur because x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and A = 1 (inexactly).

<

if (y > 1)
A-1 = y-1 (inexactly).
Avoid a possible underflow caused by y/A. For casinh this would be legitimate, but will be picked up by invoking atan2 later on. For cacos this would not be legitimate.
B = (|z+I| - |z-I|) / 2 = y/A

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z+I".
System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z-I".
Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1). sqrt_A2my2 = sqrt(Amy*(A+y))
fp is of order x^2, and fm = x/2. A = 1 (inexactly).
Underflow will not occur because x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN and x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and A = y (inexactly).
y < RECIP_EPSILON. So the following scaling should avoid any underflow problems.

<

if (y < 1)
fm = 1-y >= DBL_EPSILON, fp is of order x^2, and A = 1 (inexactly).
static NPY_INLINE void _clog_for_large_values c ( @type @  x,
@type @  y,
@type @*  rr,
@type @*  ri 
) [static]
Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.

Docutils System Messages

System Message: ERROR/3 (<string>, line 1); backlink Undefined substitution referenced: "z".
Avoid overflow in hypot() when x and y are both very large. Divide x and y by E, and then add 1 to the logarithm. This depends on E being larger than sqrt(2). Dividing by E causes an insignificant loss of accuracy; however this method is still poor since it is uneccessarily slow.
Avoid overflow when x or y is large. Avoid underflow when x or y is small.
static NPY_INLINE type _sum_squares c ( @type @  x,
@type @  y 
) [static]
sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow). Assumes x*x and y*y will not overflow. Assumes x and y are finite. Assumes y is non-negative. Assumes fabs(x) >= DBL_EPSILON.
Avoid underflow when y is small.

Variable Documentation

ctype npy_catanh c = {1.0@C@, 0.0} [static]
begin repeat
#type = npy_float, npy_double, npy_longdouble# #ctype = npy_cfloat,npy_cdouble,npy_clongdouble# c = f, , l# C = F, , L# #TMAX = FLT_MAX, DBL_MAX, LDBL_MAX# #TMIN = FLT_MIN, DBL_MIN, LDBL_MIN# #TMANT_DIG = FLT_MANT_DIG, DBL_MANT_DIG, LDBL_MANT_DIG# #TEPS = FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON# #precision = 1, 2, 3#
System Message: SEVERE/4 (<string>, line 1)
Title overline & underline mismatch.

==========================================================
 Constants
=========================================================
System Message: SEVERE/4 (<string>, line 1)
Title overline & underline mismatch.

==========================================================
 Custom implementation of missing complex C99 functions
=========================================================
algorithm from cpython, rev. d86f5686cef9
The usual formula for the real part is log(hypot(z.real, z.imag)). There are four situations where this formula is potentially problematic:
(1) the absolute value of z is subnormal. Then hypot is subnormal, so has fewer than the usual number of bits of accuracy, hence may have large relative error. This then gives a large absolute error in the log. This can be solved by rescaling z by a suitable power of 2.
(2) the absolute value of z is greater than DBL_MAX (e.g. when both z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX) Again, rescaling solves this.
(3) the absolute value of z is close to 1. In this case it's difficult to achieve good accuracy, at least in part because a change of 1ulp in the real or imaginary part of z can result in a change of billions of ulps in the correctly rounded answer.
(4) z = 0. The simplest thing to do here is to call the floating-point log with an argument of 0, and let its behaviour (returning -infinity, signaling a floating-point exception, setting errno, or whatever) determine that of c_log. So the usual formula is fine here.
Always use this function because of the multiplication for small integer powers, but in the body use cpow if it is available.
private function for use in npy_pow{f, ,l}
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
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.
casinh(z) = z + O(z^3) as z -> 0
casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity The above formula works for the imaginary part as well, because Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)

System Message: ERROR/3 (<string>, line 6) Unexpected indentation.

<blockquote> as z -> infinity, uniformly in y</blockquote>

real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y). Assumes x and y are not NaN, and one of x and y is larger than RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use the code creal(1/z), because the imaginary part may produce an unwanted underflow. This is only called in a context where inexact is always raised before the call, so no effort is made to avoid or force inexact.
divide by zeros should yield a complex inf or nan

Referenced by _BOOL(), _ctype_power(), array_diagonal(), c(), extint_safe_binop(), extint_shl_128(), HALF_copysign(), HALF_sign(), NumPyOS_ascii_isspace(), PyArray_Diagonal(), and TIMEDELTA_mq_m_divide().

__COMP_NPY_UNUSED npy_float tiny = 3.9443045e-31f [static]