Fast integer solution of x(x-1)/2 = c

Question

Given a non-negative integer c, I need an efficient algorithm to find the largest integer x such that

x*(x-1)/2 <= c

Equivalently, I need an efficient and reliably accurate algorithm to compute:

x = floor((1 + sqrt(1 + 8*c))/2)        (1)

For the sake of defineteness I tagged this question C++, so the answer should be a function written in that language. You can assume that c is an unsigned 32 bit int.

Also, if you can prove that (1) (or an equivalent expression involving floating-point arithmetic) always gives the right result, that's a valid answer too, since floating-point on modern processors can be faster than integer algorithms.

Answer 1

If you're willing to assume IEEE doubles with correct rounding for all operations including square root, then the expression that you wrote (plus a cast to double) gives the right answer on all inputs.

Here's an informal proof. Since c is a 32-bit unsigned integer being converted to a floating-point type with a 53-bit significand, 1 + 8*(double)c is exact, and sqrt(1 + 8*(double)c) is correctly rounded. 1 + sqrt(1 + 8*(double)c) is accurate to within one ulp, since the last term being less than 2**((32 + 3)/2) = 2**17.5 implies that the unit in the last place of the latter term is less than 1, and thus (1 + sqrt(1 + 8*(double)c))/2 is accurate to within one ulp, since division by 2 is exact.

The last piece of business is the floor. The problem cases here are when (1 + sqrt(1 + 8*(double)c))/2 is rounded up to an integer. This happens if and only if sqrt(...) rounds up to an odd integer. Since the argument of sqrt is an integer, the worst cases look like sqrt(z**2 - 1) for positive odd integers z, and we bound

z - sqrt(z**2 - 1) = z * (1 - sqrt(1 - 1/z**2)) >= 1/(2*z)

by Taylor expansion. Since z is less than 2**17.5, the gap to the nearest integer is at least 1/2**18.5 on a result of magnitude less than 2**17.5, which means that this error cannot result from a correctly rounded sqrt.

Adopting Yakk's simplification, we can write

(uint32_t)(0.5 + sqrt(0.25 + 2.0*c))

without further checking.

Answer 2

If we start with the quadratic formula, we quickly reach sqrt(1/4 + 2c), round up at 1/2 or higher.

Now, if you do that calculation in floating point, there can be inaccuracies.

There are two approaches to deal with these inaccuracies. The first would be to carefully determine how big they are, determine if the calculated value is close enough to a half for them to be important. If they aren't important, simply return the value. If they are, we can still bound the answer to being one of two values. Test those two values in integer math, and return.

However, we can do away with that careful bit, and note that sqrt(1/4 + 2c) is going to have an error less than 0.5 if the values are 32 bits, and we use doubles. (We cannot make this guarantee with floats, as by 2^31 the float cannot handle +0.5 without rounding).

In essense, we use the quadratic formula to reduce it to two possibilities, and then test those two.

uint64_t eval(uint64_t x) {
  return x*(x-1)/2;
}
unsigned solve(unsigned c) {
  double test = sqrt( 0.25 + 2.*c );
  if ( eval(test+1.) <= c )
    return test+1.
  ASSERT( eval(test) <= c );
  return test;
}

Note that converting a positive double to an integral type rounds towards 0. You can insert floors if you want.

Answer 3

This may be a bit tangential to your question. But, what caught my attention is the specific formula. You are trying to find the triangular root of T_{n - 1} (where T_n is the nth triangular number).

I.e.:

T_n = n * (n + 1) / 2

and

T_n - n = T_{n - 1} = n * (n - 1) / 2

From the nifty trick described here, for T_n we have:

n = int(sqrt(2 * c))

Looking for n such that T_{n - 1} ≤ c in this case doesn't change the definition of n, for the same reason as in the original question.

Computationally, this saves a few operations, so it's theoretically faster than the exact solution (1). In reality, it's probably about the same.

Neither this solution or the one presented by David are as "exact" as your (1) though.

floor((1 + sqrt(1 + 8*c))/2) (blue) vs int(sqrt(2 * c)) (red) vs Exact (white line)

---

floor((1 + sqrt(1 + 8*c))/2) (blue) vs int(sqrt(0.25 + 2 * c) + 0.5 (red) vs Exact (white line)

My real point is that triangular numbers are a fun set of numbers that are connected to squares, pascal's triangle, Fibonacci numbers, et. al.

As such there are loads of identities around them which might be used to rearrange the problem in a way that didn't require a square root.

Of particular interest may be that T_n + T_{n - 1} = n²

I'm assuming you know that you're working with a triangular number, but if you didn't realize that, searching for triangular roots yields a few questions such as this one which are along the same topic.