Add complex number support to expm1

[kgryte]

This PR

• adds complex number support to expm1 by documenting special cases. The exponential function is an entire function in the complex plane. Thus, the function does not have branch cuts.
• updates the input and output array data types to be any floating-point data type, not just real-valued floating-point data types.
• derives special cases from C99 exp and tested against NumPy (script found below).

Details

importnumpyasnpimportmathdefis_equal_float(x, y):
"""Test whether two floating-point numbers are equal with special consideration for zeros and NaNs. Parameters ---------- x : float First input number. y : float Second input number. Returns ------- bool Boolean indicating whether two floating-point numbers are equal. Examples -------- >>> is_equal_float(0.0, -0.0) False >>> is_equal_float(-0.0, -0.0) True """# Handle +-0:ifx==0.0andy==0.0:
returnmath.copysign(1.0, x) ==math.copysign(1.0, y)
# Handle NaNs:ifx!=x:
returny!=y# Everything else, including infinities:returnx==ydefis_equal(x, y):
"""Test whether two complex numbers are equal with special consideration for zeros and NaNs. Parameters ---------- x : complex First input number. y : complex Second input number. Returns ------- bool Boolean indicating whether two complex numbers are equal. Examples -------- >>> import numpy as np >>> is_equal(complex(np.nan, np.nan), complex(np.nan, np.nan)) True """returnis_equal_float(x.real, y.real) andis_equal_float(x.imag, y.imag)
# Strided array consisting of input values and expected values:values= [
complex(0.0, 0.0), # 0complex(0.0, 0.0), # 0complex(-0.0, 0.0), # 1complex(-0.0, 0.0), # 1complex(1.0, np.inf), # 2complex(np.nan, np.nan), # 2complex(1.0, np.nan), # 3complex(np.nan, np.nan), # 3complex(np.inf, 0.0), # 4complex(np.inf, 0.0), # 4, seems to be a bug in NumPy, as it returns (inf+nanj), vs np.exp(complex(np.inf, 0.0))-1.0 == (inf+0j)complex(-np.inf, 1.0), # 5complex(-1.0, 0.0), # 5complex(np.inf, 1.0), # 6complex(np.inf, np.inf), # 6complex(-np.inf, np.inf), # 7complex(-1.0, 0.0), # 7, seems to be a bug in NumPy, as it returns (nan+nanj), vs np.exp(complex(-np.inf, np.inf))-1.0 == (-1+0j)complex(np.inf, np.inf), # 8complex(np.inf, np.nan), # 8, seems to be a bug in NumPy, as it returns (nan+nanj), vs np.exp(complex(np.inf, np.inf))-1.0 == (inf+nanj)complex(-np.inf, np.nan), # 9complex(-1.0, 0.0), # 9, seems to be a bug in NumPy, as it returns (nan+nanj), vs np.exp(complex(-np.inf, np.nan))-1.0 == (-1+0j)complex(np.inf, np.nan), # 10complex(np.inf, np.nan), # 10, seems to be a bug in NumPy, as it returns (nan+nanj), vs np.exp(complex(np.inf, np.nan))-1.0 == (inf+nanj)complex(np.nan, 0.0), # 11complex(np.nan, 0.0), # 11, seems to be a bug in NumPy, as it returns (nan+nanj), vs np.exp(complex(np.nan, 0.0))-1.0 == (nan+0j)complex(np.nan, 1.0), # 12complex(np.nan, np.nan), # 12complex(np.nan, np.nan), # 13complex(np.nan, np.nan) # 13
]
foriinrange(len(values)//2):
j=i*2v=values[j]
e=values[j+1]
actual=np.expm1(v)
print('Value: {value}'.format(value=str(v)))
print('Actual: {actual}'.format(actual=str(actual)))
print('Naive: {naive}'.format(naive=str(np.exp(v)-1.0)))
print('Expected: {expected}'.format(expected=str(e)))
print('Equal: {is_equal}'.format(is_equal=str(is_equal(actual, e))))
print('\n')

Value: 0j
Actual: 0j
Naive: 0j
Expected: 0j
Equal: True
Value: (-0+0j)
Actual: (-0+0j)
Naive: 0j
Expected: (-0+0j)
Equal: True
/path/to/cexpm1.py:113: RuntimeWarning: invalid value encountered in expm1
actual = np.expm1(v)
Value: (1+infj)
Actual: (nan+nanj)
/path/to/cexpm1.py:116: RuntimeWarning: invalid value encountered in exp
print('Naive: {naive}'.format(naive=str(np.exp(v)-1.0)))
Naive: (nan+nanj)
Expected: (nan+nanj)
Equal: True
Value: (1+nanj)
Actual: (nan+nanj)
Naive: (nan+nanj)
Expected: (nan+nanj)
Equal: True
Value: (inf+0j)
Actual: (inf+nanj)
Naive: (inf+0j)
Expected: (inf+0j)
Equal: False
Value: (-inf+1j)
Actual: (-1+0j)
Naive: (-1+0j)
Expected: (-1+0j)
Equal: True
Value: (inf+1j)
Actual: (inf+infj)
Naive: (inf+infj)
Expected: (inf+infj)
Equal: True
Value: (-inf+infj)
Actual: (nan+nanj)
Naive: (-1+0j)
Expected: (-1+0j)
Equal: False
Value: (inf+infj)
Actual: (nan+nanj)
Naive: (inf+nanj)
Expected: (inf+nanj)
Equal: False
Value: (-inf+nanj)
Actual: (nan+nanj)
Naive: (-1+0j)
Expected: (-1+0j)
Equal: False
Value: (inf+nanj)
Actual: (nan+nanj)
Naive: (inf+nanj)
Expected: (inf+nanj)
Equal: False
Value: (nan+0j)
Actual: (nan+nanj)
Naive: (nan+0j)
Expected: (nan+0j)
Equal: False
Value: (nan+1j)
Actual: (nan+nanj)
Naive: (nan+nanj)
Expected: (nan+nanj)
Equal: True
Value: (nan+nanj)
Actual: (nan+nanj)
Naive: (nan+nanj)
Expected: (nan+nanj)
Equal: True

Notes

• NumPy currently fails for 46 complex number special cases. It's behavior is inconsistent with np.exp(z)-1, as documented in the script above.

[kgryte]

As no objections have been raised to the changes introduced in this PR and the changes follow established conventions (C99 and equivalent exp(x)-1 behavior), will merge. Revisions to special cases can be made in follow-up PRs.

[honno]

Seems torch doesn't support complex inputs for expm1

>>>torch.expm1(torch.as_tensor(1.+1.j))
RuntimeError: "expm1_vml_cpu"notimplementedfor'ComplexFloat'

exp however does take complex inputs. Couldn't see any relevant issues to this on PyTorch's tracker.