How do I properly set the precision of std::complex<mpf

I wanted to have an std::complex<mpf_class> with a specific precision (400), and read in the complex from std::cin. I didn't know how to set the precision in the declaration of the complex, so I tried this first:

std::complex<mpf_class > c;
c.real().set_prec(400);
c.imag().set_prec(400);
std::cin >> c;
std::cout << std::setprecision(400) << c;

But with the input (0.12345678123456781234567812345678,0), I got the output (0.123456781234567812346,0). This isn't right because 400 bits of precision should be able to hold all the digits of the input. I wanted to make sure that std::complex could hold higher precision variables so I tried this:

mpf_class a(0, 400);
mpf_class b(0, 400);
std::cin >> a >> b;
std::complex<mpf_class> c(a, b);
std::cout << std::setprecision(400) << c;

With the same input, 0.12345678123456781234567812345678 0, the output was (0.12345678123456781234567812345678,0). So it held all the digits, but I had to both declare and read the complex parts first; not ideal.

But I was still trying to figure out what was causing the first code segment to not work. One difference I noticed between the first two code segments was that the precision in the first was set after the declaration of the mpf_class, while the second segment sets the precision during the mpf_class declaration. I just wanted to check that you could update the precision of an mpf_class so I tried this:

mpf_class a(0, 32);
a.set_prec(400);
std::cin >> a;
std::cout << std::setprecision(400) << a;

With the input 0.12354678123456781234567812345678, the output was the same which means the precision can be updated.

Finally I thought there might be some issue with the complex stream extraction operator and high precision decimals. So I set the precision of the mpf_class before declaring the complex, then declared the complex, then read in the complex:

mpf_class a, b;
a.set_prec(400);
b.set_prec(400);
std::complex<mpf_class> c(a, b);
std::cin >> c;
std::cout << std::setprecision(400) << c;

But with this the weirdest thing happened. With input (0.12345678123456781234567812345678,0) the output was (0.1234567812345678123456781234567799999999999999999999999999273060010609330153935727917471428455421640301869890510557285753902098601227437,0)

TLDR

There are three steps to reading in a high precision complex: set precision, read values, declare complex. These steps can be permutated in a variety of orders. Obviously the setting of the precision has to be done before reading the values. But then there are still three different orderings with three different behaviors. Code segment 1 shows that the ordering of (declare complex, set precision, read values) doesn't work because the precision didn't get updated. Code segment 2 shows that the ordering of (set precision, read values, declare complex) works as expected, but is inconvenient because you have to declare the complex components separately. Finally code segment 4 shows that the ordering of (set precision, declare complex, read values) introduces weird rounding errors.

My questions are these: Can anyone explain these discrepancies? How can I declare an std::complex<mpf_class> with a given precision?

CodePudding user response：

The main point is that .real() returns a copy of the content of the complex, not a pointer/reference to it. Thus, statement:

c.real().set_prec(400);

is changing the precision of the copy, not the precision of the value inside the complex.

Example:

std::complex< std::string > a("foo","bar");
std::string s = a.real();
cout <<  "1-" << a << endl;
cout <<  "2-" << s << endl;

a.real().clear();
s.clear();
cout << "3-" << a << endl;
cout << "4-" << s << endl;

gives an output of:

1-(foo,bar)
2-foo
3-(foo,bar)
4-

being .clear() "ignored" after .real().

See definition of read:

// std::complex<T>::real
constexpr T real() const;

About the fourth code example in the OP question: operator >> of complex class usually does:

reads the real and imaginary parts as separate data
call the class constructor to create a new complex and assign it to the argument.

(see here as example, statement: __x = complex<_Tp>(__u, __v);)

thus, all previous value (including .set_prec effects) are lost.

Suggestion: using the mpf library, create your own C class of real numbers with the correct precision (assign it in the constructor) and use it as argument to complex template. The precision can be hard-coded or stored in a static variable or the class can be a template.

CodePudding user response：

Thanks to pasaba por aqui for explaining why the first code segment doesn't work. I have thought about it a little more and decided this might be an appropriate solution:

std::complex<mpf_class> c (mpf_class(0,400), mpf_class(0,400));
std::cin >> c;
std::cout << std::setprecision(400) << c;

This code segment properly sets the precision of the complex's individual parts without having to declare the individual parts before hand.

However, it still doesn't solve the weird rounding error problem that we saw in the OP code segment 4. If anyone has a better solution to the OP or can explain how the rounding error is introduced/how to avoid it, please comment or submit a new answer. I will keep this post open for a couple days and if I don't see anything else I will accept this answer as the answer to the OP.