Atanh returning values slightly outside the specified error bound

[MattyMuir]

Describe the bug
The Sleef_atanhd1_u10purec function is returning values with an error of slightly more than 1ULP for arguments around 10^-17.

I've also experienced similar issues with other atanh functions, such as Sleef_atanhd4_u10avx2. This may affect others also.

Command lines and logs
I built SLEEF using

cmake -G Ninja -D CMAKE_C_COMPILER="C:/Program Files/LLVM/bin/clang.exe" -D SLEEF_BUILD_GNUABI_LIBS=OFF -D SLEEF_BUILD_TESTS=OFF -B ./build
cd build
cmake --build .

Configure log: CMakeConfigureLog.txt
Build log: cmake_build_log.txt

To Reproduce
Here's a minimal reproducer

#include <iostream>

#include <sleef-config.h>
#include <sleef.h>

int main()
{
    double x = 8.0e-17;
    double y = Sleef_atanhd1_u10purec(x);
    if (y < x)
        std::cout << "Error!\n";
}

The CMake used to build this example is

cmake_minimum_required(VERSION 3.13)
project(sleefAtanhMinimal)
add_executable(sleefAtanhMinimal main.cpp)

target_link_libraries(sleefAtanhMinimal PUBLIC "C:/Users/matty/Desktop/sleef/build/lib/sleef.lib")
target_include_directories(sleefAtanhMinimal PUBLIC "C:/Users/matty/Desktop/sleef/build/include")

(I just hardcoded the paths to SLEEF for simplicity)
And I built using the commands

cmake -G Ninja -D CMAKE_C_COMPILER="C:/Program Files/LLVM/bin/clang.exe" -D CMAKE_CXX_COMPILER="C:/Program Files/LLVM/bin/clang.exe" -B ./build
cd build
cmake --build .

Expected behavior
I expected that "Error!" would not be printed to the console. Since atanh(x) > x for all positive x, by returning a value strictly less than x, it's produced a result more than 1ULP away from the exact value.

Environment

• Building natively
• Architecture: x86_64
• OS: Windows 10 Home version 22H2 build 19045.4780
• CMake Version: 3.26.3
• Makefile Generator: Ninja
• Compiler: Clang version 18.1.6, target x86_64-pc-windows-msvc

I'm more than happy to provide any other information you need.

[joanaxcruz]

Hey, thank you so much for the thorough issue report, it helps us a lot. We will try reproducing this inaccuracy and come back to you with more questions if we need.

[joanaxcruz]

I built sleef with the configuration you posted here.

When i link with an extended version of your program:

int main(int argc, char **argv) {
  double in = 8.0e-17;
  double result = Sleef_atanhd1_u10purec(in);
  printf("Result was %a \n", result);
  if (result < in)
     printf("Error!\n");

  double reference_result = atanh(in);
  printf("Reference Result was %a \n", reference_result);
}

I obtain the following output:

Result was 0x1.70ef54646d496p-54 
Error!
Reference Result was 0x1.70ef54646d497p-54 

In the documentation (https://sleef.org/2-references/libm/#sleef_atanh_u10), we advertise that the error bound for double-precision atanh functions is 1.0 ULP is what we see in this case.
Could you check on your side if this is the output that you're observing? Is there any case where you see a difference that is bigger than 1ULP between the results and the reference.

Note: I am using linux on x86 instead of windows, compiling with gcc-11 - there is a difference in our environment, so differences are possible.

[MattyMuir]

I observe the same output on my machine.

In your program you compare against a reference value computed at double precision, but if you compare against a reference value computed at higher precision, the error is slightly more than 1.0ULP.

I thought that ULP error was measured relative to the exact result, and I've been running my tests under that assumption and all the other functions have passed.

I think a relevant portion of the docs is:
"All the functions in the library are thoroughly tested and confirmed that the evaluation error is within the designed limit by comparing the returned values against high-precision evaluation"

Let me know if I've misunderstood the way that error is measured, and I'll close the issue.

[joanaxcruz]

Hi, so I made an adaptation to the previous program I sent:

int main(int argc, char **argv) {
  double in = 8.0e-17;
  double result = Sleef_atanhd1_u10purec(in);
  printf("Result was %a \n", result);
  if (result < in)
     printf("Error!\n");

  double reference_dp = atanh(in);
  printf("Reference double precision result was %a \n", reference_dp);

  double reference_ld = atanhl(in);
  printf("Reference long double precision result was %a \n", reference_ld);

  mpfr_t in_mpfr, out_mpfr;
  mpfr_init2 (in_mpfr, 200);
  mpfr_init2 (out_mpfr, 200);
  mpfr_set_d (in_mpfr, in, GMP_RNDN);
  mpfr_atanh(out_mpfr, in_mpfr, GMP_RNDN);
  double reference_mpfr = mpfr_get_ld(out_mpfr, GMP_RNDN);
  printf("Reference mpfr result was %a \n", reference_mpfr);
  mpfr_clear (in_mpfr);
  mpfr_clear (out_mpfr);
}

And this is the output I obtain:

Result was 0x1.70ef54646d496p-54 
Error!
Reference double precision result was 0x1.70ef54646d497p-54 
Reference long double precision result was 0x1.70ef54646d497p-54 
Reference mpfr result was 0x1.70ef54646d497p-54 

In this example I use mpfr_atanh (rounded to long double) and GLIBC’s atanhl as my high accuracy references, and still don't see more than 1.0ULP error. Could you let me know which high precision references you are using?

[MattyMuir]

Hi I've written an extended reproducer using MPFR as my reference.

mpfr_prec_t GlobalPrecision = 200;

void PrintULPError(double approx, mpfr_srcptr reference)
{
	// If you're worried about precision being lost here,
	// convert approx and approxNudged to MPFR and perform the subtraction after,
	// you'll get the same result
	double approxNudged = nextafter(approx, INFINITY);
	double ulp = approxNudged - approx;
	
	// Compute error as (reference - approx) / ulp
	mpfr_t err;
	mpfr_init2(err, GlobalPrecision);
	mpfr_sub_d(err, reference, approx, MPFR_RNDN);
	mpfr_div_d(err, err, ulp, MPFR_RNDN);

	mpfr_printf("ULP Error: %.20RDf", err);

	mpfr_clear(err);
}

int main()
{
	double x = 8e-17;

	// Compute approx atanh(x)
	double approx = Sleef_atanhd1_u10purec(x);

	// Convert x to mpfr
	mpfr_t xMpfr;
	mpfr_init2(xMpfr, GlobalPrecision);
	mpfr_set_d(xMpfr, x, MPFR_RNDN);

	// Compute reference atanh(x)
	mpfr_t reference;
	mpfr_init2(reference, GlobalPrecision);
	mpfr_atanh(reference, xMpfr, MPFR_RNDN);

	// Display error
	PrintULPError(approx, reference);

	mpfr_clear(xMpfr);
	mpfr_clear(reference);
}

This is the output I observe

ULP Error: 1.00000000000000001384

[joanaxcruz]

Hi, just letting you know we have been able to reproduce the issue. Currently deciding how to proceed about this, will follow up on this soon.

[joanaxcruz]

After having a look at tester2dp.c where we test Sleef_atanhd1_u10purec, I noticed that even though we use high precision for intermediate error calculations, we truncate the error at double precision, and hence we do not consider errors that cannot be captured by this.