Incorrect const float values in xvalues.cpp files

[maksqwe]

Incorerct const float values in x_values.cpp files?

There are several constants from xfvalues.cpp, xlvalues.cpp, xvalues.cpp:

extern const float _FXbig = (NBITS + 1) * 347L / 1000;
extern const double _Xbig = (NBITS + 1) * 347L / 1000;
extern const long double _LXbig = (NBITS + 1) * 347L / 1000;

A diffrerence only in NBITS value and result type.
But all this constants (_FXbig, _Xbig, _LXbig) ignores fractional part.
for example:
extern const float _FXbig = (NBITS + 1) * 347L / 1000; // _FXbig == 8.0f;

In fact _FXbig should equel to 8.32800... if use cast to float:
extern const float _FXbig = static_cast<float>(NBITS + 1) * 347L / 1000;

This values are used to calculate sin(complex)/cos(complex): xlsinh.cpp

_CRTIMP2_PURE long double __CLRCALL_PURE_OR_CDECL _LSinh(long double x, long double y) {
    // compute y * sinh(x), |y| <= 1
	...
    default: // finite
		...
        if (x < _LRteps._Long_double) {
            x *= y; // x tiny
        } else if (x < 1.0L) {
            long double w = x * x;
            x += x * w * _LPoly(w, p, NP - 1);
            x *= y;
        } else if (x < _LXbig) { // worth adding in exp(-x)
            _LExp(&x, 1.0L, -1);
            x = y * (x - 0.25L / x);
        } else {
		...
    }
	...
}

Because I have not test cases to check, I just reported this issue as a question.
Does the loss of accuracy affect the result?
@statementreply

[statementreply]

It looks like a bug.

These values are used as a threshold to decide whether exp(-x) is so tiny compared to exp(x) that it could be discarded when computing sinh(x) and cosh(x). Using a value too small would affect the accuracy of the result.

Test case:

G:\Temp>type xvalues_xbig.cpp
#include <cmath>
#include <complex>
#include <cstdio>
using namespace std;

int main() {
    const float real_sinh_f = sinh(8.0f);
    const float complex_sinh_f = real(sinh(8.0f + 0if));
    printf("     expected          = 1490.47882578... (0x1.749ea514e...p+10)\n");
    printf("     sinh(8.0f)        = %-16.9g (%.6a)\n", real_sinh_f, real_sinh_f);
    printf("real(sinh(8.0f + 0if)) = %-16.9g (%.6a)\n", complex_sinh_f, complex_sinh_f);
    printf("\n");

    const float real_cosh_f = cosh(8.0f);
    const float complex_cosh_f = real(cosh(8.0f + 0if));
    printf("     expected          = 1490.47916125... (0x1.749eaa93f...p+10)\n");
    printf("     cosh(8.0f)        = %-16.9g (%.6a)\n", real_cosh_f, real_cosh_f);
    printf("real(cosh(8.0f + 0if)) = %-16.9g (%.6a)\n", complex_cosh_f, complex_cosh_f);
    printf("\n");

    const double real_sinh_d = sinh(18.0);
    const double complex_sinh_d = real(sinh(18.0 + 0i));
    printf("     expected          = 32829984.568665247954... (0x1.f4f22091940bb255...p+24)\n");
    printf("     sinh(18.0)        = %-24.17g (%a)\n", real_sinh_d, real_sinh_d);
    printf("real(sinh(18.0 + 0i))  = %-24.17g (%a)\n", complex_sinh_d, complex_sinh_d);
    printf("\n");

    const double real_cosh_d = cosh(18.0);
    const double complex_cosh_d = real(cosh(18.0 + 0i));
    printf("     expected          = 32829984.568665263184... (0x1.f4f22091940bf3bf...p+24)\n");
    printf("     cosh(18.0)        = %-24.17g (%a)\n", real_cosh_d, real_cosh_d);
    printf("real(cosh(18.0 + 0i))  = %-24.17g (%a)\n", complex_cosh_d, complex_cosh_d);
}

G:\Temp>cl /EHsc /W4 /WX xvalues_xbig.cpp
用于 x64 的 Microsoft (R) C/C++ 优化编译器 19.27.29016 版
版权所有(C) Microsoft Corporation。保留所有权利。

xvalues_xbig.cpp
Microsoft (R) Incremental Linker Version 14.27.29016.0
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:xvalues_xbig.exe
xvalues_xbig.obj

G:\Temp>.\xvalues_xbig
     expected          = 1490.47882578... (0x1.749ea514e...p+10)
     sinh(8.0f)        = 1490.47888       (0x1.749ea6p+10)
real(sinh(8.0f + 0if)) = 1490.47913       (0x1.749eaap+10)

     expected          = 1490.47916125... (0x1.749eaa93f...p+10)
     cosh(8.0f)        = 1490.47913       (0x1.749eaap+10)
real(cosh(8.0f + 0if)) = 1490.47913       (0x1.749eaap+10)

     expected          = 32829984.568665247954... (0x1.f4f22091940bb255...p+24)
     sinh(18.0)        = 32829984.568665247       (0x1.f4f22091940bbp+24)
real(sinh(18.0 + 0i))  = 32829984.568665259       (0x1.f4f22091940bep+24)

     expected          = 32829984.568665263184... (0x1.f4f22091940bf3bf...p+24)
     cosh(18.0)        = 32829984.568665262       (0x1.f4f22091940bfp+24)
real(cosh(18.0 + 0i))  = 32829984.568665259       (0x1.f4f22091940bep+24)

[StephanTLavavej]

Ordinarily I would like to call the CRT for math machinery when possible - but its C99 complex is powered by code from the same source, so I expect it'll be equally affected.

[MattStephanson]

@StephanTLavavej
I'm interested in working on this issue, but are you saying that it would need a separate fix in the CRT? Or will patching x?values.cpp fix it in both places?

Also, in addition to the integer truncation issue, I think the calculation (NBITS + 1) has an off-by-one error. Either the formula was derived to make the exp(-x) smaller than one ULP, instead of half an ULP, or the implementer didn't realize that NBITS is the physical, not effective, size of the mantissa. I tested every value near the crossover, and it doesn't matter for float but does make a difference for double. The current formula gives a crossover of x=18.391. The exp(-x)/2 term is 5.15e-9, which is more than half the gap of 7.45e-9. Pushing the crossover to 18.738 gives no difference compared to always doing the subtraction.

[StephanTLavavej]

@MattStephanson The CRT and the STL have separate codebases, so fixing the STL's xvalues.cpp files will fix the bug for std::complex without needing the CRT to be fixed simultaneously. However, the CRT's <complex.h> will be affected until it receives a separate fix. (The CRT's affected file appears to be src\appcrt\tran\_values.c, not shipped publicly AFAICT. The history here is that MSVC's STL was originally licensed from Dinkumware in the 1990s, while the code that became the UCRT was originally written by MS back in the 1980s. However, when we added support for C99 complex circa 2010, we licensed that code from Dinkumware - hence the shared origin.)

Apologies for the delayed reply - I fell behind on my GitHub notifications. We look forward to a PR! 😸