8000 BUG: Can't rely on C99 complex arithmetic. · numpy/numpy@050dca3 · GitHub
[go: up one dir, main page]
Skip to content

Commit 050dca3

Browse files
WarrenWeckesserWarrenWeckesser
committed
BUG: Can't rely on C99 complex arithmetic.
* Add utility functions npy_cmul@c@(z, w) and npy_cdiv@c@(z, w). * Use them in log1p for complex arithmetic.
1 parent 2511f87 commit 050dca3

File tree

3 files changed

+29
-2
lines changed

3 files changed

+29
-2
lines changed

‎numpy/core/include/numpy/npy_math.h

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -457,6 +457,9 @@ static inline npy_clongdouble npy_cpackl(npy_longdouble x, npy_longdouble y)
457457
double npy_cabs(npy_cdouble z);
458458
double npy_carg(npy_cdouble z);
459459

460+
npy_cdouble npy_cmul(npy_cdouble z, npy_cdouble w);
461+
npy_cdouble npy_cdiv(npy_cdouble z, npy_cdouble w);
462+
460463
npy_cdouble npy_cexp(npy_cdouble z);
461464
npy_cdouble npy_clog(npy_cdouble z);
462465
npy_cdouble npy_cpow(npy_cdouble x, npy_cdouble y);
@@ -485,6 +488,9 @@ npy_cdouble npy_catanh(npy_cdouble z);
485488
float npy_cabsf(npy_cfloat z);
486489
float npy_cargf(npy_cfloat z);
487490

491+
npy_cfloat npy_cmulf(npy_cfloat z, npy_cfloat w);
492+
npy_cfloat npy_cdivf(npy_cfloat z, npy_cfloat w);
493+
488494
npy_cfloat npy_cexpf(npy_cfloat z);
489495
npy_cfloat npy_clogf(npy_cfloat z);
490496
npy_cfloat npy_cpowf(npy_cfloat x, npy_cfloat y);
@@ -514,6 +520,9 @@ npy_cfloat npy_catanhf(npy_cfloat z);
514520
npy_longdouble npy_cabsl(npy_clongdouble z);
515521
npy_longdouble npy_cargl(npy_clongdouble z);
516522

523+
npy_clongdouble npy_cmull(npy_clongdouble z, npy_clongdouble w);
524+
npy_clongdouble npy_cdivl(npy_clongdouble z, npy_clongdouble w);
525+
517526
npy_clongdouble npy_cexpl(npy_clongdouble z);
518527
npy_clongdouble npy_clogl(npy_clongdouble z);
519528
npy_clongdouble npy_cpowl(npy_clongdouble x, npy_clongdouble y);

‎numpy/core/src/npymath/npy_math_complex.c.src

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -123,6 +123,18 @@ cdiv@c@(@ctype@ a, @ctype@ b)
123123
* Custom implementation of missing complex C99 functions
124124
*=========================================================*/
125125

126+
@ctype@
127+
npy_cmul@c@(@ctype@ z, @ctype@ w)
128+
{
129+
return cmul@c@(z, w);
130+
}
131+
132+
@ctype@
133+
npy_cdiv@c@(@ctype@ z, @ctype@ w)
134+
{
135+
return cdiv@c@(z, w);
136+
}
137+
126138
#ifndef HAVE_CABS@C@
127139
@type@
128140
npy_cabs@c@(@ctype@ z)

‎numpy/core/src/umath/funcs.inc.src

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -366,7 +366,10 @@ nc_log1p@c@(@ctype@ *x, @ctype@ *r)
366366
* every computer scientist should know about floating-point
367367
* arithmetic".
368368
*/
369-
@ctype@ u = *x + 1.0;
369+
370+
/* u = 1 + *x */
371+
@ctype@ u = npy_cpack@c@(1 + x_re, x_im);
372+
370373
if (npy_creal@c@(u) - 1.0 == x_re) {
371374
/*
372375
* Don't bother multiplying by *x / (u - 1), because that quotient
@@ -375,7 +378,10 @@ nc_log1p@c@(@ctype@ *x, @ctype@ *r)
375378
*r = npy_clog@c@(u);
376379
}
377380
else {
378-
*r = npy_clog@c@(u) * (*x / (u - 1.0));
381+
/* um1 = u - 1 = (1 + *x) - 1 */
382+
@ctype@ um1 = npy_cpack@c@(npy_creal@c@(u) - 1.0, x_im);
383+
/* *r = log(u) * (*x / (u - 1)) */
384+
*r = npy_cmul@c@(npy_clog@c@(u), npy_cdiv@c@(*x, um1));
379385
}
380386
}
381387
return;

0 commit comments

Comments
 (0)
0