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().