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

Jump to bottom
Closed
wants to merge 4 commits into from
Closed

BUG: umath: Fix log1p for complex inputs. #24416

wants to merge 4 commits into from

Conversation

WarrenWeckesser
Copy link
Member
@WarrenWeckesser WarrenWeckesser commented Aug 14, 2023 •
edited
Loading

Reimplement the complex log1p function. Use the log1p trick from Theorem 4 of Goldberg's paper "What every computer scientist should know about floating-point arithmetic". Include special handling of an input with imaginary part ±0.0 to ensure the sign of the imaginary part of the result is correct and consistent with the complex log function.

Closes gh-22609.

< 8000 div class="TimelineItem js-targetable-element" data-team-hovercards-enabled id="event-10089661030">
@github-actions github-actions bot added the 00 - Bug label Aug 14, 2023
@WarrenWeckesser WarrenWeckesser mentioned this pull request Aug 14, 2023
Merged
@WarrenWeckesser
Copy link
Member Author

It looks like the failures on Emscripten/Pyodide, musllinux_x86_64 and Cygwin are because the underlying complex log function clog provided by the standard math library is not very accurate.

If I run a C program that computes w = clog(CMPLX(1.0, 1e-12)) and compile it with musl-gcc on my Linux machine, clog returns w = 1e-12j. The correct value is wtrue = 5e-25 + 1e-12j, which gives a relative error of abs(w - wtrue)/abs(wtrue) = 5e-13. (musl-gcc --version shows x86_64-linux-gnu-gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0.) With clang on a Mac and with gcc on Linux, the value is computed correctly. For some other inputs, the relative error of the musl-based clog can be as high as 5e-10.

The unit test tolerances can be loosened a bit, but I'm tempted to make the tolerances of the tests platform-dependent.

@mattip
Copy link
Member
mattip commented Aug 15, 2023

the underlying complex log function clog provided by the standard math library is not very accurate.

We can blocklist clog on more platforms (currently only MSVC), and use our implementation. I guess we should be reproducing the cpython cmath tests to catch cases like this.

@WarrenWeckesser WarrenWeckesser marked this pull request as draft June 5, 2024 00:56
@WarrenWeckesser
Copy link
Member Author
WarrenWeckesser commented Jun 5, 2024 •
edited
Loading

After gh-24448 was merged, I expected all the tests to pass after I rebased, and in fact, the old failures no longer occur.

The failure in the Windows tests / x86-64, LP64 OpenBLAS (Clang-cl) job is weird. It shows an overflow warning being generated when log1p(nan + 0j) is computed. But for that input, the result is simply set to nan + nanj and returned, with no floating point arithmetic taking place, so there should be no overflow.

@WarrenWeckesser
Copy link
Member Author

Somehow these two lines

            npy_csetreal@c@(r, NAN);
            npy_csetimag@c@(r, NAN);

are causing FE_OVERFLOW to be set in the job Windows tests / x86-64, LP64 OpenBLAS (Clang-cl).

The block of code where this happens, with my print-debugging included, is

        if (npy_isnan(x_re)) {
            if (fetestexcept(FE_OVERFLOW)) {
                fprintf(stderr, "***** FE_OVERFLOW is now set, after calling npy_isnan(x_re)!\n");
                fflush(stderr);
            }
            fprintf(stderr, "nc_log1p@c@: setting re and im parts of r to NAN.\n");
            fflush(stderr);

            npy_csetreal@c@(r, NAN);
            npy_csetimag@c@(r, NAN);

            if (fetestexcept(FE_OVERFLOW)) {
                fprintf(stderr, "***** FE_OVERFLOW is now set, after setting r to (NAN, NAN)!\n");
                fflush(stderr);
            }

        }

The output captured in the test log is

nc_log1p: x_re = nan  x_im = 0.000000
nc_log1p: setting re and im parts of r to NAN.
***** FE_OVERFLOW is now set, after setting r to (NAN, NAN)!
nc_log1p: returning, r = (nan, nan)
***** FE_OVERFLOW is now set, immediately before returning!

The first output line is from a print earlier in the code, and the last two output lines are from prints that occur later.

In the code shown, the first call of fetestexcept(FE_OVERFLOW) does not succeed. The message about setting r is printed, then the functions npy_csetreal and npy_csetimag are called. After those calls, the call of fetestexcept(FE_OVERFLOW) succeeds.

Those functions just assign NAN to memory locations. How could that result in FE_OVERFLOW being set?

@WarrenWeckesser
Copy link
Member Author
WarrenWeckesser commented Jun 6, 2024 •
edited
Loading

Of course, I search the web after my making my previous comment: https://stackoverflow.com/questions/69929589/is-it-considered-normal-that-f-nan-may-cause-raising-floating-point-exceptions

TL;DR: On this particular test platform (Windows+clang), the NAN macro is implemented as an expression that evaluates to a nan, e.g. something like INFINITY * 0.0. That expression is evaluated at runtime, resulting in various floating point exception flags being set, including FE_OVERFLOW.

@WarrenWeckesser
BUG: umath: Fix log1p for complex inputs.
2b1ca5e 
Reimplement the complex log1p function.  Use the log1p trick from Theorem 4
of Goldberg's paper "What every computer scientist should know about
floating-point arithmetic".  Include special handling of an input with
imaginary part 0.0 to ensure the sign of the imaginary part of the result
is correct and consistent with the complex log function.

Closes numpygh-22609.
@WarrenWeckesser
BUG: Can't rely on C99 complex arithmetic.
91568d6 
* Add utility functions npy_cmul@c@(z, w) and npy_cdiv@c@(z, w).
* Use them in log1p for complex arithmetic.
@WarrenWeckesser
MAINT: Use NPY_NAN* instead of NAN.
5c010e6 
With some compilers, NAN is expanded to an expression that
is evaluated at runtime. That can cause various floating
point exception flags to be set, including FE_OVERFLOW.
@WarrenWeckesser
Update tests.
ef231fc
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
00 - Bug
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

BUG: Inaccurate log1p for small complex input
2 participants
@WarrenWeckesser @mattip
0