Fix #379: <complex>: asinh(-1e+307 + 2e+307i) is wrong

[statementreply]

Description

Fix a few wrong formula (possibly typo) in complex asinh, acosh and acos. (#379)

Checklist

Be sure you've read README.md and understand the scope of this repo.

If you're unsure about a box, leave it unchecked. A maintainer will help you.

• Identifiers in product code changes are properly _Ugly as per
  https://eel.is/c++draft/lex.name#3.1 or there are no product code changes.
• The STL builds successfully and all tests have passed (must be manually
  verified by an STL maintainer before automated testing is enabled on GitHub,
  leave this unchecked for initial submission).
• These changes introduce no known ABI breaks (adding members, renaming
  members, adding virtual functions, changing whether a type is an aggregate
  or trivially copyable, etc.).
• These changes were written from scratch using only this repository,
  the C++ Working Draft (including any cited standards), other WG21 papers
  (excluding reference implementations outside of proposed standard wording),
  and LWG issues as reference material. If they were derived from a project
  that's already listed in NOTICE.txt, that's fine, but please mention it.
  If they were derived from any other project (including Boost and libc++,
  which are not yet listed in NOTICE.txt), you must mention it here,
  so we can determine whether the license is compatible and what else needs
  to be done.

[Billy O'Neal (BillyONeal)]

Amy Wishnousky (@amyw-msft) Do we need to make corresponding changes to the UCRT?

[statementreply]

Amy Wishnousky (@amyw-msft) Do we need to make corresponding changes to the UCRT?

I found another correctness bug in UCRT complex (but not in STL) during testing. I didn't report that bug in #379 because I believed this is not the place to report UCRT bugs.

G:\Temp>type ahyp2.cpp
#include <complex>
#include <complex.h>
#include <cstdio>

int main() {
    using namespace std;
    const auto e_stl = asinh(2.0i);
    const auto e_ucrt = casinh(_Cbuild(0.0, 2.0));

    puts("asinh(2.0i)");
    puts("  Expected: ( 1.3169578969248167086..., 1.5707963267948966192...)");
    printf("  STL     : (% -25.17g,% -25.17g)\n", e_stl.real(), e_stl.imag());
    printf("  UCRT    : (% -25.17g,% -25.17g)\n", creal(e_ucrt), cimag(e_ucrt));
}

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

ahyp2.cpp
Microsoft (R) Incremental Linker Version 14.24.28314.0
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:ahyp2.exe
ahyp2.obj

G:\Temp>.\ahyp2
asinh(2.0i)
  Expected: ( 1.3169578969248167086..., 1.5707963267948966192...)
  STL     : ( 1.3169578969248166      , 1.5707963267948966      )
  UCRT    : ( 0                       , 1.5707963267948966      )

[Billy O'Neal (BillyONeal)]

You are correct that this isn't the place for UCRT bugs, but the implementations are similar -- I believe they both came from Dinkumware once upon a time, so when one is broken, usually they both are.

[Stephan T. Lavavej (StephanTLavavej)]

Can you provide test values that exercise each of the altered lines?

[statementreply]

Can you provide test values that exercise each of the altered lines?

Assuming 64-bit long double

L1234

assert(nearly_equal_partwise(acos(-2e+307 + 1e+307i), 2.677945044588987 - 708.3914896859491i));
assert(nearly_equal_partwise(acos(-4e+37f + 2e+37if), 2.6779451f - 87.386663if));
assert(nearly_equal_partwise(acos(-2e+307L + 1e+307il), 2.677945044588987L - 708.3914896859491il));

L1242

assert(nearly_equal_partwise(acos(2e+307 + 1e+307i), 0.4636476090008061 - 708.3914896859491i));
assert(nearly_equal_partwise(acos(4e+37f + 2e+37if), 0.4636476f - 87.386663if));
assert(nearly_equal_partwise(acos(2e+307L + 1e+307il), 0.4636476090008061L - 708.3914896859491il));

assert(nearly_equal_partwise(acos(2e+307 - 1e+307i), 0.4636476090008061 + 708.3914896859491i));
assert(nearly_equal_partwise(acos(4e+37f - 2e+37if), 0.4636476f + 87.386663if));
assert(nearly_equal_partwise(acos(2e+307L - 1e+307il), 0.4636476090008061L + 708.3914896859491il));

L1305

assert(nearly_equal_partwise(acosh(-1e+307 + 2e+307i), 708.3914896859491 + 2.0344439357957027i));
assert(nearly_equal_partwise(acosh(-2e+37f + 4e+37if), 87.386663f + 2.0344439if));
assert(nearly_equal_partwise(acosh(-1e+307L + 2e+307il), 708.3914896859491L + 2.0344439357957027il));

L1380

assert(nearly_equal_partwise(asinh(-1e+307 + 2e+307i), -708.3914896859491 + 1.1071487177940904i));
assert(nearly_equal_partwise(asinh(-2e+37f + 4e+37if), -87.386663f + 1.1071488if));
assert(nearly_equal_partwise(asinh(-1e+307L + 2e+307il), -708.3914896859491L + 1.1071487177940904il));

[Billy O'Neal (BillyONeal)]

Hello statementreply!

I'm seeing that the new implementation doesn't pass the test cases you enumerated above assuming nearly_equal_partwise tests that the results are within +/- 1 ULP. e.g. nearly_equal_partwise(acos(4e+37f + 2e+37if), 0.4636476f - 87.386663if) is off by 2 ULP. I observe that the complex implementation calls several base transcendentals that each individually have a +/- 1 ULP policies which can lead to this outcome; is that expected?

Thanks again for your contribution!

[Stephan T. Lavavej (StephanTLavavej)]

That seems expected to me; if you're concerned about the 2 ULP differences, I recommend filing an issue to investigate later.

[Billy O'Neal (BillyONeal)]

I still would like to give statementreply a couple days to comment in case they want to before merging just in case the results I observe are unexpected.

[statementreply]

Billy O'Neal (@BillyONeal) That's expected. This PR fixes some coding errors in the original implementation (which doesn't guarantee +/- 1 ULP accuracy) of complex asinh, acosh and acos. I don't intend to improve their accuracy in this PR.

[Billy O'Neal (BillyONeal)]

Sounds good, thanks again for your contribution!

[Stephan T. Lavavej (StephanTLavavej)]

This is a great example of a bug that we would have had a much harder time fixing before open-source! We really appreciate your analysis, test cases, and fix. 😸