diff options
Diffstat (limited to 'numpy/core/src')
| -rw-r--r-- | numpy/core/src/scalarmathmodule.c.src | 15 | ||||
| -rw-r--r-- | numpy/core/src/umath/loops.c.src | 19 |
2 files changed, 25 insertions, 9 deletions
diff --git a/numpy/core/src/scalarmathmodule.c.src b/numpy/core/src/scalarmathmodule.c.src index 2bcc516b1..56f1bc238 100644 --- a/numpy/core/src/scalarmathmodule.c.src +++ b/numpy/core/src/scalarmathmodule.c.src @@ -382,10 +382,17 @@ static npy_half (*_basic_half_fmod)(npy_half, npy_half); (outp)->real = (a).real * (b).real - (a).imag * (b).imag; \ (outp)->imag = (a).real * (b).imag + (a).imag * (b).real; \ } while(0) -#define @name@_ctype_divide(a, b, outp) do{ \ - @rtype@ d = (b).real*(b).real + (b).imag*(b).imag; \ - (outp)->real = ((a).real*(b).real + (a).imag*(b).imag)/d; \ - (outp)->imag = ((a).imag*(b).real - (a).real*(b).imag)/d; \ +/* Note: complex division by zero must yield some complex inf */ +#define @name@_ctype_divide(a, b, outp) do{ \ + @rtype@ d = (b).real*(b).real + (b).imag*(b).imag; \ + if (d != 0) { \ + (outp)->real = ((a).real*(b).real + (a).imag*(b).imag)/d; \ + (outp)->imag = ((a).imag*(b).real - (a).real*(b).imag)/d; \ + } \ + else { \ + (outp)->real = (a).real/d; \ + (outp)->imag = (a).imag/d; \ + } \ } while(0) #define @name@_ctype_true_divide @name@_ctype_divide #define @name@_ctype_floor_divide(a, b, outp) do { \ diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src index 5212207da..54e5ac984 100644 --- a/numpy/core/src/umath/loops.c.src +++ b/numpy/core/src/umath/loops.c.src @@ -1804,11 +1804,20 @@ C@TYPE@_divide(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func const @type@ in1i = ((@type@ *)ip1)[1]; const @type@ in2r = ((@type@ *)ip2)[0]; const @type@ in2i = ((@type@ *)ip2)[1]; - if (npy_fabs@c@(in2r) >= npy_fabs@c@(in2i)) { - const @type@ rat = in2i/in2r; - const @type@ scl = 1.0@c@/(in2r + in2i*rat); - ((@type@ *)op1)[0] = (in1r + in1i*rat)*scl; - ((@type@ *)op1)[1] = (in1i - in1r*rat)*scl; + const @type@ in2r_abs = npy_fabs@c@(in2r); + const @type@ in2i_abs = npy_fabs@c@(in2i); + if (in2r_abs >= in2i_abs) { + if (in2r_abs == 0 && in2i_abs == 0) { + /* divide by zero should yield a complex inf or nan */ + ((@type@ *)op1)[0] = in1r/in2r_abs; + ((@type@ *)op1)[1] = in1i/in2i_abs; + } + else { + const @type@ rat = in2i/in2r; + const @type@ scl = 1.0@c@/(in2r + in2i*rat); + ((@type@ *)op1)[0] = (in1r + in1i*rat)*scl; + ((@type@ *)op1)[1] = (in1i - in1r*rat)*scl; + } } else { const @type@ rat = in2r/in2i; |