BUG: umath: Fix log1p for complex inputs.

[WarrenWeckesser]

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.

[WarrenWeckesser]

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]

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]

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]

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]

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.