Yep, that's right, if you use std::pow to square a real-valued complex number in C++11, you don't get a real number

dear lord, save me from egregious numerical error

So, there I was simplifying my syntax tree in Bertini2. I wrote some code, which eliminates 0's and 1's, and reduces the depth of the tree, by flattening it. And, I discovered that some systems don't have a nice property -- that the function values at random points before and after simplification differ. By a lot. Why?

Well, consider this SO post from /u/Schaki. They ask why, when using C++11, pow(i,2) has non-zero but small imaginary part. Like, 1e-16 or so. The answer is that there is a new implementation, and set of templates, for C++11 complex arithmetic. And it's scary.

Why is this so scary? Because pure-real complex numbers don't stay pure-real. And that causes enormous amounts of drift through a long sequence of complex operations.

\(\mathbb{R}\) must remain real under those operations for which it should. This is one of the greatest lessons I have learned from implementing Bertini_real, my surface and curve decomposer. If a number is real, theoretically, then make it be as real as the metal will let it.

My solution to this problem in Bertini2 was to add an overload for pow(std::complex<double>, int) which uses successive squaring. It restores sanity, in the sense that the following two tests pass, where they wouldn't with just std::pow.

BOOST_AUTO_TEST_CASE(complex_pow)
{
    auto x = bertini::rand_complex();

    using bertini::pow;
    BOOST_CHECK_EQUAL(pow(x, 2), x*x);
}


BOOST_AUTO_TEST_CASE(complex_pow_on_real_stays_real)
{
    using bertini::pow;

    for (int ii=0; ii<100; ++ii)
    {
        auto x = dbl(bertini::RandReal());
        auto result = pow(x, 2);
        BOOST_CHECK_EQUAL(result, x*x);
        BOOST_CHECK_EQUAL(imag(result), 0);
    }
}