Radius grows exponentially for long-running calculations

[cgd8d]

Consider the program:

#include "arb.h"
#include <iostream>
int main()
{
    arb_t x;
    arb_init(x);
    arb_add_error_2exp_si(x, -30); // Add 2^(-30) to radius of x.
    arb_t y;
    arb_init(y); // y = 0, exactly.
    slong prec = 64;

    for(size_t i = 0; i < 16; i++) {
        for(size_t j = 0; j < (1ULL << 27); j++) {
            arb_add(y, y, x, prec);
        }
        std::cout << "After "<< i+1 << " * 2^27 iterations, y = " << arb_get_str(y, 10, 0) << std::endl;
    }
}

When run, you will see that the radius of y grows in a piecewise-linear fashion, but that the slope doubles every 2^29 calculations. The implication is that for calculations running longer than 2^29 calculations, the radius grows exponentially rather than linearly.

This magic number 2^29 comes, I think, from the choice of MAG_BITS = 30 in https://github.com/fredrik-johansson/arb/blob/master/mag.h#L126. An ideal solution would be to increase MAG_BITS for calculations which are expected to be long-running, but that doesn't work out-of-the-box.

Minimal fix: This behavior could simply be documented, since long-running calculations won't be wrong, just expensive due to precision required.

Preferred: Enable increasing MAG_BITS, either by deciding what the longest calculation you wish to support should be (MAG_BITS = 62?), or by making the code work for a range of values for MAG_BITS depending on user needs. I don't have any understanding of the ripple effects involved in increasing MAG_BITS, though.

Thanks for maintaining a very nice piece of code!

[fredrik-johansson]

This is a rather extreme case. In the rare application where this is an issue, could you not use divide and conquer summation?

Alternatively, instead of adding into an arb_t you could manually add the midpoints into an arb_t and the radii into another arb_t where you control the precision. That is at least assuming that you are just doing one big summation and not using the partial sums as part of the computation of new terms.

I agree that MAG_BITS = 62 on 64-bit systems would be nicer, but increasing MAG_BITS would have performance repercussions.

[fredrik-johansson]

As you said, it could be worth documenting this somewhere.

In the rare application where this is an issue, could you not use divide and conquer summation?

Recursively or also just one layer: instead of adding to the final variable y all the time, use a different arb_t variable z for the inner loop with < 2^29 terms, and add this partial sum to y at the each step of the outer loop. This way I got:

1: [+/- 0.137]
2: [+/- 0.274]
3: [+/- 0.411]
4: [+/- 0.548]
5: [+/- 0.684]
6: [+/- 0.821]
7: [+/- 0.958]
8: [+/- 1.10]
9: [+/- 1.24]
10: [+/- 1.37]
11: [+/- 1.51]
12: [+/- 1.65]
13: [+/- 1.78]
14: [+/- 1.92]
15: [+/- 2.06]
16: [+/- 2.19]

[cgd8d]

Yes, this test case is just meant to illustrate the phenomenon quickly. My own use case involves looking at long-running partial sums (looking at each of the partial sums in turn) and also at partial products. You've given some good strategies for reducing the impact of this issue, which I appreciate. (For your other suggestion to re-group the operations, I think it requires two adds per iteration -- one to update z, and one to compute temp = y + z at each iteration so further calculations can be done -- so I will profile my application to see if that is prohibitive.)

My one minor point of push-back is that you say this is a rather extreme case. The test program is, but the need to understand loss of accuracy over the course of a long-running calculations seems like a core application of something like arb. Examples of when re-grouping won't work include any system with feedback, and I can imagine that arb would be a great tool for modeling a complex system with feedback because it lets you track when floating-point errors have accumulated enough to affect the accuracy of your model. For those uses, 2^29 will come quickly, so setting the MAG_BITS=62 flag might be desirable.

But those aren't my use cases; for me, re-grouping or doing manual summation/multiplication work fine.

[fredrik-johansson]

What I meant is not that it's an extreme case to do a long-running calculation; I think it's an extreme case that additions alone are the problem. Especially if you have any kind of feedback, it's likely that you will need several bits of higher precision anyway.