`std::log1p(x)` with long double appears to use double precision

[WarrenWeckesser]

I think I can guess the source of the problem, but I want to be sure.

check_log1p.cpp

#include <cstdio>
#include <cmath>
#include <limits>

int main()
{
    printf("numeric_limits<long double>::digits = %d\n",
           std::numeric_limits<long double>::digits);
    long double x = -5.0000000000000000000000005e-35L;
    printf("%40.33Le\n", x);
    long double yld = std::log1p(x);
    printf("%40.33Le\n", yld);
    long double yd = static_cast<long double>(std::log1p(static_cast<double>(x)));
    printf("%40.33Le\n", yd);
}

Compile and run:

% em++ check_log1p.cpp -sPRINTF_LONG_DOUBLE
% node a.out.js                            
numeric_limits<long double>::digits = 113
-5.000000000000000000000000500000000e-35
-4.999999999999999638373019459082555e-35
-4.999999999999999638373019459082555e-35

The value of x is small enough that log1p(x) should equal x. The last two lines being the same shows that the calculation of std::log1p(x) is being done with double, not long double.

My guess is that, even though you provide 113 bits precision in long double, your std math library does not support that precision. Is that correct?

The silent downcast to double is a nasty wart in the API. Is there any way to enable an error or warning when that happens?

Version of emscripten/emsdk:

emcc (Emscripten gcc/clang-like replacement + linker emulating GNU ld) 3.1.61 (67fa4c16496b157a7fc3377afd69ee0445e8a6e3)
clang version 19.0.0git (https:/github.com/llvm/llvm-project 7cfffe74eeb68fbb3fb9706ac7071f8caeeb6520)
Target: wasm32-unknown-emscripten
Thread model: posix
InstalledDir: /Users/warren/repos/git/forks/emsdk/upstream/bin

Full link command and output with -v appended:

Here's the compile command with -v added:

% em++ -v check_log1p.cpp -sPRINTF_LONG_DOUBLE
 "/Users/warren/repos/git/forks/emsdk/upstream/bin/clang++" -target wasm32-unknown-emscripten -fignore-exceptions -mllvm -combiner-global-alias-analysis=false -mllvm -enable-emscripten-sjlj -mllvm -disable-lsr --sysroot=/Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot -DEMSCRIPTEN -Xclang -iwithsysroot/include/fakesdl -Xclang -iwithsysroot/include/compat -v check_log1p.cpp -c -o /var/folders/6f/wyccw81j5kj1c7v6p47zr9fr0000gn/T/emscripten_temp_3rlg3nz6/check_log1p_0.o
clang version 19.0.0git (https:/github.com/llvm/llvm-project 7cfffe74eeb68fbb3fb9706ac7071f8caeeb6520)
Target: wasm32-unknown-emscripten
Thread model: posix
InstalledDir: /Users/warren/repos/git/forks/emsdk/upstream/bin
 (in-process)
 "/Users/warren/repos/git/forks/emsdk/upstream/bin/clang-19" -cc1 -triple wasm32-unknown-emscripten -emit-obj -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name check_log1p.cpp -mrelocation-model static -mframe-pointer=none -ffp-contract=on -fno-rounding-math -mconstructor-aliases -target-cpu generic -fvisibility=hidden -debugger-tuning=gdb -fdebug-compilation-dir=/Users/warren/code-snippets/c++/test-emscripten-stuff -v -fcoverage-compilation-dir=/Users/warren/code-snippets/c++/test-emscripten-stuff -resource-dir /Users/warren/repos/git/forks/emsdk/upstream/lib/clang/19 -D EMSCRIPTEN -isysroot /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot -internal-isystem /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/wasm32-emscripten/c++/v1 -internal-isystem /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/c++/v1 -internal-isystem /Users/warren/repos/git/forks/emsdk/upstream/lib/clang/19/include -internal-isystem /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/wasm32-emscripten -internal-isystem /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include -fdeprecated-macro -ferror-limit 19 -fgnuc-version=4.2.1 -fskip-odr-check-in-gmf -fcxx-exceptions -fignore-exceptions -fexceptions -fcolor-diagnostics -iwithsysroot/include/fakesdl -iwithsysroot/include/compat -mllvm -combiner-global-alias-analysis=false -mllvm -enable-emscripten-sjlj -mllvm -disable-lsr -o /var/folders/6f/wyccw81j5kj1c7v6p47zr9fr0000gn/T/emscripten_temp_3rlg3nz6/check_log1p_0.o -x c++ check_log1p.cpp
clang -cc1 version 19.0.0git based upon LLVM 19.0.0git default target x86_64-apple-darwin22.6.0
ignoring nonexistent directory "/Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/wasm32-emscripten/c++/v1"
ignoring nonexistent directory "/Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/wasm32-emscripten"
#include "..." search starts here:
#include <...> search starts here:
 /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/fakesdl
 /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/compat
 /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include/c++/v1
 /Users/warren/repos/git/forks/emsdk/upstream/lib/clang/19/include
 /Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/include
End of search list.
 /Users/warren/repos/git/forks/emsdk/upstream/bin/clang --version
 /Users/warren/repos/git/forks/emsdk/upstream/bin/wasm-ld -o a.out.wasm -L/Users/warren/repos/git/forks/emsdk/upstream/emscripten/cache/sysroot/lib/wasm32-emscripten /var/folders/6f/wyccw81j5kj1c7v6p47zr9fr0000gn/T/emscripten_temp_3rlg3nz6/check_log1p_0.o -lGL-getprocaddr -lal -lhtml5 -lprintf_long_double-debug -lstubs-debug -lnoexit -lc-debug -ldlmalloc -lcompiler_rt -lc++-noexcept -lc++abi-debug-noexcept -lsockets -mllvm -combiner-global-alias-analysis=false -mllvm -enable-emscripten-sjlj -mllvm -disable-lsr /var/folders/6f/wyccw81j5kj1c7v6p47zr9fr0000gn/T/tmp3en6quuulibemscripten_js_symbols.so --strip-debug --export=emscripten_stack_get_end --export=emscripten_stack_get_free --export=emscripten_stack_get_base --export=emscripten_stack_get_current --export=emscripten_stack_init --export=_emscripten_stack_alloc --export=__get_temp_ret --export=__set_temp_ret --export=__wasm_call_ctors --export=_emscripten_stack_restore --export-if-defined=__start_em_asm --export-if-defined=__stop_em_asm --export-if-defined=__start_em_lib_deps --export-if-defined=__stop_em_lib_deps --export-if-defined=__start_em_js --export-if-defined=__stop_em_js --export-if-defined=main --export-if-defined=__main_argc_argv --export-if-defined=fflush --export-table -z stack-size=65536 --no-growable-memory --initial-heap=16777216 --no-entry --stack-first --table-base=1
 /Users/warren/repos/git/forks/emsdk/upstream/bin/llvm-objcopy a.out.wasm a.out.wasm --remove-section=.debug* --remove-section=producers
 /Users/warren/repos/git/forks/emsdk/upstream/bin/wasm-emscripten-finalize --dyncalls-i64 --pass-arg=legalize-js-interface-exported-helpers a.out.wasm -o a.out.wasm --detect-features
 /Users/warren/repos/git/forks/emsdk/node/18.20.3_64bit/bin/node /Users/warren/repos/git/forks/emsdk/upstream/emscripten/src/compiler.mjs /var/folders/6f/wyccw81j5kj1c7v6p47zr9fr0000gn/T/tmp4ac3r4z7.json

[sbc100]

The log1p function looks like it takes a double. Did you mean to write log1pl? Do you get the same results with log1pl?

https://github.com/emscripten-core/emscripten/blob/82b0d740924f589b320c624191497e1ab7857edc/system/lib/libc/musl/include/math.h#L269

BTW, the code for this function is part of the upstream musl project and lives in system/lib/libc/musl/src/math/log1pl.c.

[WarrenWeckesser]

If I am reading https://en.cppreference.com/w/cpp/numeric/math/log1p correctly, log1p should be overloaded to accept all three floating point types. The code works as expected with other C++ compilers.

I tried with log1pl, and still got the result corresponding to the double calculation.

[sbc100]

If I am reading https://en.cppreference.com/w/cpp/numeric/math/log1p correctly, log1p should be overloaded to accept all three floating point types. The code works as expected with other C++ compilers.

I tried with log1pl, and still got the result corresponding to the double calculation.

Sounds like the problem is perhaps with the implementation of log1pl. Perhaps you could re-write the example in pure C to show that the issue with musl (the underlying C library) rather than with libc++. Assuming that the issue can be reproduced in pure C then I think the next step would be to take a look at the implementation of log1pl to see if there is anything obvious stopping it from working with 128-bit float.

BTW, Wasm doesn't actually support 128-bit floats at the instruction level so any time you use long doulbe you will end up pulling in software fallbacks which are very slow and requires a bunch of extra code. If you can instead avoid long double in your code that would probably give you better results under wasm. is that possible in your case?

[WarrenWeckesser]

is that possible in your case?

I'm actually working on NumPy, and one of the platforms that is included in NumPy's test suite is pyodide, and pyodide is built on emscripten. I had a unit test failing on the pyodide platform, so I've had to dig into how emscripten handles floating point. The calculation that I'm working on loses precision for certain regions of the input values, and one way to avoid that loss is to using higher precision variables when computing in that region.

[WarrenWeckesser]

It looks like musl doesn't implement full long double precision (i.e. 113 bit mantissa). From log1pl.c:

#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
// TODO: broken implementation to make things compile
long double log1pl(long double x)
{
    return log1p(x);
}
#endif

Other functions (e.g. logl, expl, etc) have the same behavior.

[sbc100]

Do you have some way to disable long double support in NumPy at compile time? Since it doesn't exist at the Wasm level you might as well fall back to using whatever emulation support NumPy has internally for greater precision floating point.

[WarrenWeckesser]

NumPy doesn't have any emulation for long double. Whatever the underlying compiler provides is what is used. So the underlying implementation of long double might be: the same as double, 80 bit extended precision, IEEE float128, or IBM double-double. NumPy can and does override implementations of specific functions in the standard library if it is found that the standard library version has deficiencies. So in theory, NumPy could replace all these bad musl functions with different implementations. But that means developing (or finding license-compatible versions of) float128 versions of a lot of functions, and that should really be done in musl, not NumPy.