diff options
| author | Charles Harris <charlesr.harris@gmail.com> | 2009-07-07 01:44:43 +0000 |
|---|---|---|
| committer | Charles Harris <charlesr.harris@gmail.com> | 2009-07-07 01:44:43 +0000 |
| commit | b2c38b0c7d4f4740e5738233ed55cdfe5c38b512 (patch) | |
| tree | f22a7e2d81e7c063475d8b86ccd0321336759690 /numpy | |
| parent | fc21cc0f5e6f392a64732e925ba7194a926b3e74 (diff) | |
| download | numpy-b2c38b0c7d4f4740e5738233ed55cdfe5c38b512.tar.gz | |
Make complex division more robust against over/under flow. Fixes ticket #1159.
Diffstat (limited to 'numpy')
| -rw-r--r-- | numpy/core/src/umath/loops.c.src | 15 |
1 files changed, 12 insertions, 3 deletions
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src index 929b0bb49..8c21fefc0 100644 --- a/numpy/core/src/umath/loops.c.src +++ b/numpy/core/src/umath/loops.c.src @@ -1108,9 +1108,18 @@ 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]; - @type@ d = in2r*in2r + in2i*in2i; - ((@type@ *)op1)[0] = (in1r*in2r + in1i*in2i)/d; - ((@type@ *)op1)[1] = (in1i*in2r - in1r*in2i)/d; + if (abs(in2r) >= abs(in2i)) { + const @type@ rat = in2i/in2r; + const @type@ scl = 1/(in2r + in2i*rat); + ((@type@ *)op1)[0] = (in1r + in1i*rat)*scl; + ((@type@ *)op1)[1] = (in1i - in1r*rat)*scl; + } + else { + const @type@ rat = in2r/in2i; + const @type@ scl = 1/(in2i + in2r*rat); + ((@type@ *)op1)[0] = (in1r*rat + in1i)*scl; + ((@type@ *)op1)[1] = (in1i*rat - in1r)*scl; + } } } |