<complex>: complex exp incorrect results

[amyw-msft]

Bug Description When using std::complex, results may be off by 1 ULP.

Command-line test case Note: This test case also compares against what the user reported the Fortran result was.

amd64:

C:\Users\amyw\source\exp_fortran>type test.cpp
#include <cmath>
#include <complex>
#include <stdio.h>
#include <stdint.h>
#include <string.h>

struct test_case
{
    double initial_real;
    double initial_imag;

    const char * expected_real;
    const char * expected_imag;

    double wolfram_real;
    double wolfram_imag;
};

uint64_t dtoui64(const double d)
{
    uint64_t ui64;
    memcpy(&ui64, &d, sizeof(d));
    return ui64;
}

const char * dbl_cmp(const double real_lhs, const double imag_lhs, const double real_rhs, const double imag_rhs)
{
    if (   dtoui64(real_lhs) == dtoui64(real_rhs)
        && dtoui64(imag_lhs) == dtoui64(imag_rhs)
    ) {
        return "PASS";
    } else {
        return "FAIL";
    }
}

void test(const test_case t)
{
    printf("exp( (%.30f, %.30f) ):\n", t.initial_real, t.initial_imag);

    std::complex<double> stl_temp(t.initial_real, t.initial_imag);
    std::complex<double> stl_result;

    stl_result = exp(stl_temp);

    const char * const expected_real = t.expected_real;
    const char * const expected_imag = t.expected_imag;

    char wolfram_real[64];
    char wolfram_imag[64];

    char stl_actual_real[64];
    char stl_actual_imag[64];

    sprintf_s(wolfram_real, "%.30f", t.wolfram_real);
    sprintf_s(wolfram_imag, "%.30f", t.wolfram_imag);

    sprintf_s(stl_actual_real, "%.30f", stl_result.real());
    sprintf_s(stl_actual_imag, "%.30f", stl_result.imag());

    printf("\tWolfram Real: %s\n", wolfram_real);
    printf("\tFortran Real: %s - %s\n",    expected_real, strcmp(   expected_real, wolfram_real) ? "FAIL" : "PASS");
    printf("\tSTL     Real: %s - %s\n",  stl_actual_real, strcmp( stl_actual_real, wolfram_real) ? "FAIL" : "PASS");

    printf("\tWolfram Imag: %s\n", wolfram_imag);
    printf("\tFortran Imag: %s - %s\n",    expected_imag, strcmp(   expected_imag, wolfram_imag) ? "FAIL" : "PASS");
    printf("\tSTL     Imag: %s - %s\n",  stl_actual_imag, strcmp( stl_actual_imag, wolfram_imag) ? "FAIL" : "PASS");

    printf("\tWolfram Binary: 0x%llx + 0x%llx\n", dtoui64(t.wolfram_real), dtoui64(t.wolfram_imag));
    printf("\tSTL     Binary: 0x%llx + 0x%llx - %s\n", dtoui64(stl_result.real()), dtoui64(stl_result.imag()), dbl_cmp(t.wolfram_real, t.wolfram_imag, stl_result.real(), stl_result.imag()));
}

int main()
{
    static constexpr test_case tests[] =
    {
        {   0, 1.12813019434902717108570868731,
            "0.428350136129215020019955773023", "0.903612837933416512825601785153",
             0.4283501361292150480680021695357,  0.9036128379334165575193607327889},
        {   0, -1.95503828587548733608514339721,
            "-0.374856477775750818182132206857", "-0.927082855557990526129685804335",
             -0.374856477775750824325912673562,   -0.927082855557990581658174511326},
    };

    for (auto&& t : tests) {
        test(t);
    }

    return 0;
}

C:\Users\amyw\source\exp_fortran>cl /EHsc /W4 /WX test.cpp
Microsoft (R) C/C++ Optimizing Compiler Version 19.28.29617 for x64
Copyright (C) Microsoft Corporation.  All rights reserved.

test.cpp
Microsoft (R) Incremental Linker Version 14.28.29617.0
Copyright (C) Microsoft Corporation.  All rights reserved.

/out:test.exe
test.obj

C:\Users\amyw\source\exp_fortran>test.exe
exp( (0.000000000000000000000000000000, 1.128130194349027171085708687315) ):
        Wolfram Real: 0.428350136129215075531107004281
        Fortran Real: 0.428350136129215020019955773023 - FAIL
        STL     Real: 0.428350136129215020019955773023 - FAIL
        Wolfram Imag: 0.903612837933416512825601785153
        Fortran Imag: 0.903612837933416512825601785153 - PASS
        STL     Imag: 0.903612837933416512825601785153 - PASS
        Wolfram Binary: 0x3fdb6a16b07a6049 + 0x3fecea6578656ec6
        STL     Binary: 0x3fdb6a16b07a6048 + 0x3fecea6578656ec6 - FAIL
exp( (0.000000000000000000000000000000, -1.955038285875487336085143397213) ):
        Wolfram Real: -0.374856477775750818182132206857
        Fortran Real: -0.374856477775750818182132206857 - PASS
        STL     Real: -0.374856477775750818182132206857 - PASS
        Wolfram Imag: -0.927082855557990637151988266851
        Fortran Imag: -0.927082855557990526129685804335 - FAIL
        STL     Imag: -0.927082855557990637151988266851 - PASS
        Wolfram Binary: 0xbfd7fda6062f6600 + 0xbfedaaa9aa29b93c
        STL     Binary: 0xbfd7fda6062f6600 + 0xbfedaaa9aa29b93c - PASS

x86:

C:\Users\amyw\source\exp_fortran>test.exe
exp( (0.000000000000000000000000000000, 1.128130194349027171085708687315) ):
        Wolfram Real: 0.428350136129215075531107004281
        Fortran Real: 0.428350136129215020019955773023 - FAIL
        STL     Real: 0.428350136129215020019955773023 - FAIL
        Wolfram Imag: 0.903612837933416512825601785153
        Fortran Imag: 0.903612837933416512825601785153 - PASS
        STL     Imag: 0.903612837933416512825601785153 - PASS
        Wolfram Binary: 0x3fdb6a16b07a6049 + 0x3fecea6578656ec6
        STL     Binary: 0x3fdb6a16b07a6048 + 0x3fecea6578656ec6 - FAIL
exp( (0.000000000000000000000000000000, -1.955038285875487336085143397213) ):
        Wolfram Real: -0.374856477775750818182132206857
        Fortran Real: -0.374856477775750818182132206857 - PASS
        STL     Real: -0.374856477775750818182132206857 - PASS
        Wolfram Imag: -0.927082855557990637151988266851
        Fortran Imag: -0.927082855557990526129685804335 - FAIL
        STL     Imag: -0.927082855557990526129685804335 - FAIL
        Wolfram Binary: 0xbfd7fda6062f6600 + 0xbfedaaa9aa29b93c
        STL     Binary: 0xbfd7fda6062f6600 + 0xbfedaaa9aa29b93b - FAIL

Expected behavior Our results should agree with results from Wolfram Alpha (which is also what is returned by libstdc++ and libc++)

STL version

Microsoft Visual Studio Community 2019 Version 16.8.2

Additional context This affects the UCRT (MSFT-31014699) as well, which has a different implementation that agrees with MSVC STL. Also tracked by DevCom-1271319 and Microsoft-internal VSO-1255553.

[leuri397]

I've investigated this issue and it is not related to complex numbers. It is related to double precision and parsing a double number.

On hex representation of values you can see, that the difference between wolfram value and STL value is just one bit. Actual value for real part in first test is 0.4283501361292150480680021695357 (which is returned by wolfram), but when it was converted to double from string it became 0.428350136129215075531107004281, but it is not closest approximation to real value (probably because it just run out of precise significant digits and gave up parsing on ...4806...), but when it was calculated as cosine (when real part is 0 it simplifies to cos(imag)) of double value (which was precise double representation) it gets closer approximation of 0.428350136129215020019955773023.

All in all it is not a bug of complex, it is just question of double precision (which is related to atof() function implementation in UCRT library)