`pow(std::complex<double>(randdouble,0), 2)`

\(\notin \mathbb{R}\)

### August 25, 2017

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);
}
}
```