8000 BUG: umath: Fix log1p for complex inputs. by WarrenWeckesser · Pull Request #24416 · numpy/numpy · GitHub
[go: up one dir, main page]
Skip to content

BUG: umath: Fix log1p for complex inputs. #24416

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 4 commits into from
Closed

BUG: umath: Fix log1p for complex inputs. #24416

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2b1ca5e
BUG: umath: Fix log1p for complex inputs.
WarrenWeckesser Aug 14, 2023
91568d6
BUG: Can't rely on C99 complex arithmetic.
WarrenWeckesser Aug 14, 2023
5c010e6
MAINT: Use NPY_NAN* instead of NAN.
WarrenWeckesser Jun 6, 2024
ef231fc
Update tests.
WarrenWeckesser Jun 22, 2024
File filter

Filter by extension

Filter by extension
8000
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions numpy/_core/include/numpy/npy_math.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -457,6 +457,9 @@ static inline npy_clongdouble npy_cpackl(npy_longdouble x, npy_longdouble y)
double npy_cabs(npy_cdouble z);
double npy_carg(npy_cdouble z);

npy_cdouble npy_cmul(npy_cdouble z, npy_cdouble w);
npy_cdouble npy_cdiv(npy_cdouble z, npy_cdouble w);

npy_cdouble npy_cexp(npy_cdouble z);
npy_cdouble npy_clog(npy_cdouble z);
npy_cdouble npy_cpow(npy_cdouble x, npy_cdouble y);
Expand Down Expand Up @@ -485,6 +488,9 @@ npy_cdouble npy_catanh(npy_cdouble z);
float npy_cabsf(npy_cfloat z);
float npy_cargf(npy_cfloat z);

npy_cfloat npy_cmulf(npy_cfloat z, npy_cfloat w);
npy_cfloat npy_cdivf(npy_cfloat z, npy_cfloat w);

npy_cfloat npy_cexpf(npy_cfloat z);
npy_cfloat npy_clogf(npy_cfloat z);
npy_cfloat npy_cpowf(npy_cfloat x, npy_cfloat y);
Expand Down Expand Up @@ -514,6 +520,9 @@ npy_cfloat npy_catanhf(npy_cfloat z);
npy_longdouble npy_cabsl(npy_clongdouble z);
npy_longdouble npy_cargl(npy_clongdouble z);

npy_clongdouble npy_cmull(npy_clongdouble z, npy_clongdouble w);
npy_clongdouble npy_cdivl(npy_clongdouble z, npy_clongdouble w);

npy_clongdouble npy_cexpl(npy_clongdouble z);
npy_clongdouble npy_clogl(npy_clongdouble z);
npy_clongdouble npy_cpowl(npy_clongdouble x, npy_clongdouble y);
Expand Down
12 changes: 12 additions & 0 deletions numpy/_core/src/npymath/npy_math_complex.c.src
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,18 @@ cdiv@c@(@ctype@ a, @ctype@ b)
* Custom implementation of missing complex C99 functions
*=========================================================*/

@ctype@
npy_cmul@c@(@ctype@ z, @ctype@ w)
{
return cmul@c@(z, w);
}

@ctype@
npy_cdiv@c@(@ctype@ z, @ctype@ w)
{
return cdiv@c@(z, w);
}

#ifndef HAVE_CABS@C@
@type@
npy_cabs@c@(@ctype@ z)
Expand Down
59 changes: 56 additions & 3 deletions numpy/_core/src/umath/funcs.inc.src
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,6 +268,7 @@ npy_ObjectClip(PyObject *arr, PyObject *min, PyObject *max) {
* #ctype = npy_cfloat, npy_cdouble, npy_clongdouble#
* #ftype = npy_float, npy_double, npy_longdouble#
* #c = f, ,l#
* #C = F, ,L#
*/

static void
Expand Down Expand Up @@ -310,9 +311,61 @@ nc_log@c@(@ctype@ *x, @ctype@ *r)
static void
nc_log1p@c@(@ctype@ *x, @ctype@ *r)
{
@ftype@ l = npy_hypot@c@(npy_creal@c@(*x) + 1,npy_cimag@c@(*x));
npy_csetimag@c@(r, npy_atan2@c@(npy_cimag@c@(*x), npy_creal@c@(*x) + 1));
npy_csetreal@c@(r, npy_log@c@(l));
@ftype@ x_re = npy_creal@c@(*x);
@ftype@ x_im = npy_cimag@c@(*x);

if (x_im == 0.0) {
/*
* Imaginary part of *x is +/- 0.0. Use the real-valued function
* log1p or log to compute the real part of the result. If the
* input is on the branch cut, the imaginary part of the result is
* +/- pi, otherwise it is +/- 0.0 (i.e. same as x_im).
*/
if (npy_isnan(x_re)) {
npy_csetreal@c@(r, NPY_NAN@C@);
npy_csetimag@c@(r, NPY_NAN@C@);
}
else if (x_re >= -1) {
/*
* Note that if x_re == -1, this will generate a "divide by zero"
* RuntimeWarning and set the real part of the result to -inf
* (consistent with log(0.0 +/- 0.0j)).
*/
npy_csetreal@c@(r, npy_log1p@c@(x_re));
npy_csetimag@c@(r, x_im);
}
else {
/*
* On the branch cut.
*/
npy_csetreal@c@(r, npy_log@c@(-(1 + x_re)));
npy_csetimag@c@(r, npy_copysign@c@(NPY_PI@c@, x_im));
}
}
else {
/*
* Use the log1p trick given as Theorem 4 of Goldberg's paper "What
* every computer scientist should know about floating-point
* arithmetic".
*/

/* u = 1 + *x */
@ctype@ u = npy_cpack@c@(1 + x_re, x_im);

if (npy_creal@c@(u) - 1.0 == x_re) {
/*
* Don't bother multiplying by *x / (u - 1), because that quotient
* is 1, and the complex division might introduce a (small) error.
*/
*r = npy_clog@c@(u);
}
else {
/* um1 = u - 1 = (1 + *x) - 1 */
@ctype@ um1 = npy_cpack@c@(npy_creal@c@(u) - 1.0, x_im);
/* *r = log(u) * (*x / (u - 1)) */
*r = npy_cmul@c@(npy_clog@c@(u), npy_cdiv@c@(*x, um1));
}
}
return;
}

Expand Down
184 changes: 178 additions & 6 deletions numpy/_core/tests/test_umath.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2102,6 +2102,7 @@ def test_strided_float32(self):
assert_array_almost_equal_nulp(np.sin(x_f32_large[::jj]), sin_true[::jj], nulp=2)
assert_array_almost_equal_nulp(np.cos(x_f32_large[::jj]), cos_true[::jj], nulp=2)


class TestLogAddExp(_FilterInvalids):
def test_logaddexp_values(self):
x = [1, 2, 3, 4, 5]
Expand Down Expand Up @@ -2152,13 +2153,184 @@ def test_log1p(self):
assert_almost_equal(ncu.log1p(0.2), ncu.log(1.2))
assert_almost_equal(ncu.log1p(1e-6), ncu.log(1+1e-6))

def test_special(self):
# Special cases that we test for equality. The floating point warnings
# triggered by some of these values are ignored.
@pytest.mark.parametrize(
'z, wref',
[(np.nan, np.nan),
(np.inf, np.inf),
(-1.0, -np.inf),
(-2.0, np.nan),
(-np.inf, np.nan),
(complex(np.nan, 0.0), complex(np.nan, np.nan)),
(complex(np.nan, 1), complex(np.nan, np.nan)),
(complex(1, np.nan), complex(np.nan, np.nan)),
(complex(np.inf, 1), complex(np.inf, 0.0)),
(complex(np.inf, -1), complex(np.inf, -0.0)),
(complex(-np.inf, 1), complex(np.inf, np.pi)),
(complex(-np.inf, -1), complex(np.inf, -np.pi)),
(complex(-np.inf, 0.0), complex(np.inf, np.pi)),
(complex(-np.inf, -0.0), complex(np.inf, -np.pi)),
(complex(0, np.inf), complex(np.inf, np.pi/2)),
(complex(0, -np.inf), complex(np.inf, -np.pi/2)),
(complex(-1, 0), complex(-np.inf, 0))]
)
def test_special(self, z, wref):
with np.errstate(invalid="ignore", divide="ignore"):
assert_equal(ncu.log1p(np.nan), np.nan)
assert_equal(ncu.log1p(np.inf), np.inf)
assert_equal(ncu.log1p(-1.), -np.inf)
assert_equal(ncu.log1p(-2.), np.nan)
assert_equal(ncu.log1p(-np.inf), np.nan)
w = np.log1p(z)
assert_equal(w, wref)

# Test w = log1p(z) for complex double z (np.complex128).
# Reference values were computed with mpmath, e.g.
# from mpmath import mp
# mp.dps = 200
# wref = complex(mp.log1p(z))
@pytest.mark.parametrize(
'z, wref',
[
# Near 0+0j:
(0.0, 0.0 + 0.0j),
(3.4259e-13 + 6.71894e-08j,
3.448472077361198e-13 + 6.718939999997688e-08j),
(1e-280 + 0j, 1e-280 + 0j),
(1e-18 + 0j, 1e-18 + 0j),
(1e-18 + 1e-12j, 1.0000005e-18 + 1e-12j),
(-1e-15 + 3e-8j, -5.5e-16 + 3.000000000000002e-08j),
(1e-15 + 3e-8j, 1.4499999999999983e-15 + 2.999999999999996e-08j),
(1e-50 + 1e-200j, 1e-50 + 1e-200j),
(1e-200 - 1e-200j, 1e-200 - 1e-200j),
(1e-18j, 5.0000000000000005e-37 + 1e-18j),
(-1.25e-8 + 5e-12j, -1.2500000078124988e-08 + 5.0000000625e-12j),
(-8e-8 - 0.0004j, 3.200000005993663e-15-0.0004000000106666662j),
# Close to or on the unit circle |z+1| = 1:
(-4.999958e-05 - 0.009999833j,
-7.0554155328678184e-15 - 0.009999999665816696j),
(-0.57113 - 0.90337j, 3.4168883248419116e-06 - 1.1275564209486122j),
(-0.011228922063957758 + 0.14943813247359922j,
-4.499779370610757e-17 + 0.15j),
(-1.9999999995065196 - 3.141592653092555e-05j,
-1.4766035946815242e-16 - 3.1415612376632573j),
(1e-18 + 0.1j, 0.0049751654265840425 + 0.09966865249116204j),
(0.1 + 1e-18j, 0.09531017980432487+9.090909090909091e-19j),
(-1.25e-05 + 0.005j, 7.812499999381789e-11 + 0.005000020833177082j),
(-0.01524813-0.173952j, -2.228118716056777e-06-0.1748418364650139j),
# Very large real or imaginary part:
(1e200 + 1e200j, 460.8635921890891 + 0.7853981633974483j),
(-1 + 1e250j, 575.6462732485114 + 1.5707963267948966j),
(1e250 + 1j, 575.6462732485114 + 1e-250j),
(1e275 + 1e-225j, 633.2109005733626 + 0j),
# Other:
(1.0, 0.6931471805599453 + 0j),
(-0.75 + 0j, -1.3862943611198906 + 0j),
(0.2 + 0.3j, 0.21263386770217205 + 0.24497866312686414j),
(-0.978249978133276 - 0.015379941691430407j,
-3.6254002951579887 - 0.6154905511361795j),
]
)
def test_complex_double(self, z, wref):
w = np.log1p(z)
assert_allclose(w.real, wref.real, rtol=1e-15)
assert_allclose(w.imag, wref.imag, rtol=1e-15)

# Test w = log1p(z) for complex float z (np.complex64).
# Reference values were computed with mpmath, e.g.
# from mpmath import mp
# mp.dps = 200
# wref = np.complex64(mp.log1p(z))
@pytest.mark.parametrize(
'z, wref',
[(np.complex64(0.0), np.complex64(0.0)),
(np.complex64(1e-10 + 3e-6j), np.complex64(1.045e-10 + 3e-06j)),
(np.complex64(-1e-8 - 2e-5j), np.complex64(-9.8e-09 - 2e-05j)),
(np.complex64(-2e-32 + 3e-32j), np.complex64(-2e-32 + 3e-32j)),
(np.complex64(-0.57113 - 0.90337j),
np.complex64(3.4470238e-06 - 1.1275564j)),
(np.complex64(-0.5 + 0.866025j),
np.complex64(-3.7479518e-07+1.0471973j)),
(np.complex64(-3 - 0.5j), np.complex64(0.7234595-2.896614j)),
(np.complex64(2.0 - 0.5j), np.complex64(1.1123117-0.16514868j)),
(np.complex64(3e31 - 4e31j), np.complex64(72.98958 - 0.9272952j))]
)
def test_complex_single(self, z, wref):
w = np.log1p(z)
rtol = 5*np.finfo(np.float32).eps
assert_allclose(w.real, wref.real, rtol=rtol)
assert_allclose(w.imag, wref.imag, rtol=rtol)

@pytest.mark.skipif(np.finfo(np.longdouble).machep != -63,
reason='reference values are for 80 bit '
'extended precision')
def test_longdouble_is_80_bit_extended(self):
# Passing strings to np.longdouble() here is essential for
# this test.
z = np.longdouble('-0.2') + np.longdouble('0.6')*1j
w = np.log1p(z)
wref = (np.longdouble('1.084202172485504434e-20') +
np.longdouble('0.6435011087932843868')*1j)
rtol = 5*np.finfo(np.longdouble).eps
assert_allclose(w.real, wref.real, rtol=rtol)
assert_allclose(w.imag, wref.imag, rtol=rtol)

# The reference values for the IEEE float128 test were computed
# with mpmath, e.g.:
#
# from mpmath import mp
# mp.prec = 1024 # bits of precision
# with mp.workprec(113): # Emulate IEEE float128 input
# z = mp.mpc('-0.2', '0.6')
# w = mp.log1p(z) # Compute with high precision.
# with mp.workprec(113):
# print(w.real)
# print(w.imag)
#
@pytest.mark.skipif(np.finfo(np.longdouble).machep != -112,
reason='reference values are for IEEE float128')
@pytest.mark.skipif(sys.platform == 'emscripten',
reason="emscripten's math library does not support "
"IEEE float128")
@pytest.mark.parametrize(
'z, wref',
[(np.longdouble('-0.2') + np.longdouble('0.6')*1j,
np.longdouble('-1.92592994438723585305597794258493e-35') +
np.longdouble('0.643501108793284386802809228717323')*1j),
(np.longdouble('3e-23') + np.longdouble('-5e-22')*1j,
np.longdouble('3.000000000000000000012455e-23') +
np.longdouble('-4.99999999999999999999985e-22')*1j),
(np.longdouble('-1.3') + np.longdouble('1.9')*1j,
np.longdouble('0.654166409825089380175052108173541') +
np.longdouble('1.72739820377691197435359801635447')*1j)]
)
def test_longdouble_is_ieee_float128(self, z, wref):
w = np.log1p(z)
rtol = 5*np.finfo(np.longdouble).eps
assert_allclose(w.real, wref.real, rtol=rtol)
assert_allclose(w.imag, wref.imag, rtol=rtol)

@pytest.mark.parametrize('typ', [np.complex64, np.complex128])
def test_branch_cut(self, typ):
x = -1.5
zpos = typ(x, 0.0)
wpos = np.log1p(zpos)
realtype = zpos.real.dtype.type
eps = np.finfo(realtype).eps
# On some platforms, atan2(0.0, -0.5) is not exactly equal to np.pi.
assert_allclose(wpos.imag, realtype(np.pi), rtol=2*eps)
zneg = typ(x, -0.0)
wneg = np.log1p(zneg)
assert_allclose(wneg.imag, realtype(-np.pi), rtol=2*eps)
assert wpos.real == wneg.real

@pytest.mark.parametrize('x', [-0.5, -1e-12, -1e-18, 0.0, 1e-18, 1e-12])
@pytest.mark.parametrize('typ', [np.complex64, np.complex128])
def test_imag_zero(self, x, typ):
# Test real inputs with x > -1 and the imaginary part +/- 0.0.
zpos = typ(x, 0.0)
wpos = np.log1p(zpos)
assert_equal(wpos.imag, 0.0)
zneg = typ(x, -0.0)
wneg = np.log1p(zneg)
assert_equal(wneg.imag, -0.0)
assert wpos.real == wneg.real


class TestExpm1:
Expand Down
Loading
0