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

Changes from 1 commit
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
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Update tests.
  • Loading branch information
@WarrenWeckesser
WarrenWeckesser committed Jun 22, 2024
commit ef231fc186e1f49035d76c516d4dc6826e6e444c
120 changes: 100 additions & 20 deletions numpy/_core/tests/test_umath.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2180,74 +2180,154 @@ def test_special(self, z, wref):
w = np.log1p(z)
assert_equal(w, wref)

# Test w = log1p(z) for complex z (np.complex128).
# 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',
[(1e-280 + 0j, 1e-280 + 0j),
[
# 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-18 + 0.1j, 0.0049751654265840425 + 0.09966865249116204j),
(-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),
(3.4259e-13 + 6.71894e-08j,
3.448472077361198e-13 + 6.718939999997688e-08j),
(0.1 + 1e-18j, 0.09531017980432487+9.090909090909091e-19j),
(-0.57113 - 0.90337j, 3.4168883248419116e-06 - 1.1275564209486122j),
(0.2 + 0.3j, 0.21263386770217205 + 0.24497866312686414j),
(-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),
(-0.75 + 0j, -1.3862943611198906 + 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, wref, rtol=1e-15)
assert_allclose(w.real, wref.real, rtol=1e-15)
assert_allclose(w.imag, wref.imag, rtol=1e-15)

# Test w = log1p(z) for complex z (np.complex64).
# 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(1e-10 + 3e-6j), np.complex64(1.045e-10 + 3e-06j)),
[(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)
assert_allclose(w, wref, rtol=1e-6)
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)

def test_branch_cut(self):
@pytest.mark.parametrize('typ', [np.complex64, np.complex128])
def test_branch_cut(self, typ):
x = -1.5
zpos = complex(x, 0.0)
zpos = typ(x, 0.0)
wpos = np.log1p(zpos)
assert wpos.imag == np.pi
zneg = complex(x, -0.0)
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 wneg.imag == -np.pi
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])
def test_imag_zero(self, x):
@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 = complex(x, 0.0)
zpos = typ(x, 0.0)
wpos = np.log1p(zpos)
assert_equal(wpos.imag, 0.0)
zneg = complex(x, -0.0)
zneg = typ(x, -0.0)
wneg = np.log1p(zneg)
assert_equal(wneg.imag, -0.0)
assert wpos.real == wneg.real
Expand Down
Loading
0