I have a power series with all terms non-negative which I want to evaluate to some arbitrarily set precision p (the length in binary digits of a MPFR floating-point mantissa). The result should be faithfully rounded. The issue is that I don't know when should I stop adding terms to the result variable, that is, how do I know when do I already have p + 32 accurate summed bits of the series? 32 is just an arbitrarily chosen small natural number meant to facilitate more accurate rounding to p binary digits.

This is my original series

0 <= h <= 1
series_orig(h) := sum(n = 0, +inf, a(n) * h^n)

But I actually need to calculate an arbitrary derivative of the above series (m is the order of the derivative):

series(h, m) := sum(n = m, +inf, a(n) * (n - m + 1) * ... * n * h^(n - m))

The rational number sequence a is defined like so:

a(n) := binomial(1/2, n)^2
      = (((2*n)!/(n!)) / (n! * 4^n * (2*n - 1)))^2

So how do I know when to stop summing up terms of series?

Is the following maybe a good strategy?

  1. compute in p * 4 (which is assumed to be greater than p + 32).
  2. at each point be able to recall the current partial sum and the previous one.
  3. stop looping when the previous and current partial sums are equal if rounded to precision p + 32.
  4. round to precision p and return.


I'm doing this with MPFI, an interval arithmetic addon to MPFR. Thus the [mpfi] tag.

Attempts to get relevant formulas and equations

Guided by Eric in the comments, I have managed to derive a formula for the required working precision and an equation for the required number of terms of the series in the sum.

A problem, however, is that a nice formula for the required number of terms is not possible.

Someone more mathematically capable might instead be able to achieve a formula for a useful upper bound, but that seems quite difficult to do for all possible requested result precisions and for all possible values of m (the order of the derivative). Note that the formulas need to be easily computable so they're ready before I start computing the series.

Another problem is that it seems necessary to assume the worst case for h (h = 1) for there to be any chance of a nice formula, but this is wasteful if h is far from the worst case, that is if h is close to zero.

  • 1
    Define what you mean by “accurate rounding.” In general, what is called “correctly rounded” requires being able to calculate the function to arbitrary precision regardless of the final number of digits desired. This is because, if the ideal result is near a transition point for rounding, you have to calculate the function precisely enough that the error is less than the distance between the transition point and the calculated value. And that distance has no lower limit greater than zero. Jul 29, 2021 at 11:19
  • 1
    For example, if you want two digits, and the ideal result is near 10.1 (using binary), it should be rounded to 10 if it is 10.1 or less and to 11 if it is greater than 10.1. If the ideal result is 10.1000…0001, where that last 1 appears after a thousand bits, then you need to calculate the ideal result to a thousand bits to make the rounding decision. Jul 29, 2021 at 11:20
  • 1
    On the other hand, if you will settle for “faithful rounding,” which means, for ideal results that are not exactly representable, producing either of the two nearest representable results (the one just above the ideal result or the one just below it), then it suffices to determine the error between the calculated result and the ideal result is less than 1 ULP. Jul 29, 2021 at 11:21
  • 1
    Okay, so now that the criterion is faithful rounding, what you need to show in order to conclude the calculated result is faithfully rounded is that the total error in the calculated result is less than 1 ULP. The total error includes all the rounding errors made in the additions, multiplications, and other operations used to calculate the result plus the sum of the unevaluated terms of the series. Jul 29, 2021 at 11:26
  • 1
    If MPFR has a routine to calculate n!, I suspect it may be correctly rounded, but I cannot point to documentation for that. a(n) looks well-behaved numerically. So I expect the complete rounding error in calculating the terms to be readily controllable by evaluating them with enough precision. You need information on the domain for h and how big n can get, and you need a bound on the sum of the unevaluated terms of the series as a function of n (the index of the last evaluated term). Jul 29, 2021 at 12:49


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.