14
\$\begingroup\$

When sieving primes up to some 64-bit number n I need to determine an upper bound for the greatest potential factor of n, which is the greatest integer that is not greater than the (true) square root of n.

Computing this as std::sqrt(double(n)) is fraught with pitfalls. The standard IEEE double has only 53 bits in its significand, which means that some rounding is often inevitable. One consequence is that the result cannot safely be converted to uint32_t even though the true result could never exceed 32 bits. For example, sqrt(double(uint64_t(0) - 1)) usually rounds up to \$2^{32}\$ and if that is cast to uint32_t then the result becomes 0.

Currently I'm using a function like this:

uint32_t max_factor32 (uint64_t n)
{
   double r = std::sqrt(double(n));

   if (r < UINT32_MAX)
   {
      uint32_t r32 = uint32_t(r);

      return r32 - (uint64_t(r32) * r32 > n);
   }

   return UINT32_MAX;
}

Obviously, this works only if the value returned by sqrt() is reasonably close to the true result, so that the cast to integer yields either the desired result or a number that is exactly one greater.

I think this code should be reasonably safe, even in the face of floating-point hardware/software operating with unknown rounding modes and at an unknown level of strictness. But is that truly so?

Different rounding modes should not affect the correctness of the function, but what about different compiler settings for optimisation and floating point strictness? It seems to me that normally a compiler should only err on the side of excess precision by not rounding intermediate results to IEEE double, in which case my code should work fine as it is. Can different compiler switches throw a monkey wrench into my scheme or is the code safely portable to worlds and compilers unknown?

P.S.: what I want to avoid is the 'tail wagging the dog' style caveats like "this code works only with a strictly-conforming compiler and with strict floating-point precision enabled". On the other hand, if the function max_factor32() were to mess up quietly then the result of the containing program would likely be wrong, and even crashes are possible in cases like the above example where 0 is returned when the true result would be UINT32_MAX.

\$\endgroup\$
1
  • 1
    \$\begingroup\$ "When sieving primes up to some 64-bit number n I need to determine ..." --> I suggest good alternative exists that do not do this computation. Yet since that original problem is not posted, just have to wait until it is. \$\endgroup\$ Commented Dec 8, 2023 at 2:47

6 Answers 6

6
\$\begingroup\$

Your code is OK for prime sieving.

If the rounding mode is nearest-ties-to-even and std::sqrt is correctly rounded (as required by IEEE 754), your code returns

  • \$\lfloor\sqrt{n}\rfloor + 1\$ for inputs in the range 0xfffffffdfffffc00 .. 0xfffffffe00000000, and
  • \$\lfloor\sqrt{n}\rfloor\$ for all other inputs.

If you make the if-condition r <= UINT32_MAX instead of <, the code will always return \$\lfloor\sqrt{n}\rfloor\$.
Your version doesn't consider the few cases where the correct result is UINT32_MAX - 1 but sqrt(double(n)) rounds to UINT32_MAX from below.

Here's the proof:

  • We want to compute \$r = \lfloor\sqrt{n}\rfloor\$, which we can rewrite to r*r <= n < (r+1)*(r+1).
  • Casting uint64_t → double and std::sqrt are lossy but non-decreasing operations, so we have
    sqrt(r*r) <= sqrt(n) <= sqrt((r+1)*(r+1)). Both inequalities are now <=.
    We have to special-case UINT32_MAX because (UINT32_MAX+1)*(UINT32_MAX+1) doesn't fit in uint64_t.
  • If the rounding mode is nearest-ties-to-even and std::sqrt is correctly rounded, then
    x = sqrt((uint64_t)x * x) for any uint32_t x. (You can test it exhaustively.)
  • So we can rewrite the inequalities to r <= sqrt(n) <= r+1 and check only these two cases.

Some considerations:

  • Directed rounding modes are messier.
    The proof is the same for rounding mode toward +∞. But when rounding toward −∞ or toward 0, we only get sqrt(double((uint64_t)x * x)) ∈ {x-1, x} for x >= 94906267. The final inequality r-1 <= sqrt(double(n)) <= r+1 means that you have to look for errors in both directions.
  • If std::sqrt isn't correctly rounded, the final step breaks down and you're better off with an exact integer square root.

Here's my version (only for the nearest-ties-to-even rounding mode). I like to handle the special case first.

uint32_t max_factor32(uint64_t n)
{
    // handle special case before doing the sqrt
    if (n >= uint64_t(UINT32_MAX) * UINT32_MAX) return UINT32_MAX;

    // now it fits in 32 bits
    uint32_t r = std::sqrt(n);
    return r - (uint64_t(r) * r > n);
}
\$\endgroup\$
4
  • \$\begingroup\$ This works OK as long as double has at least half the precision of uint64_t n. I suspect it fails for uint128_t n \$\endgroup\$ Commented Dec 8, 2023 at 2:44
  • \$\begingroup\$ It seems to work up to 2¹⁰⁶. After that, x = sqrt((uint128_t)x * x) breaks down because the doubles start skipping odd integers after 2⁵³. \$\endgroup\$
    – Řrřola
    Commented Dec 9, 2023 at 11:11
  • 1
    \$\begingroup\$ @chux-ReinstateMonica I've thought about it some more. It definitely doesn't work for n = (2⁵⁴+1)² because we're choosing between r and r-1 but the nearest possible doubles (outputs of sqrt) are 2⁵⁴ and 2⁵⁴+4, so there's no way we can hit the correct result. I couldn't find a smaller counterexample. \$\endgroup\$
    – Řrřola
    Commented Dec 9, 2023 at 21:00
  • \$\begingroup\$ Good thought experiment. Using FP for integer problems often runs into such portability issues. Hence, for integer tasks, I use FP with caution and prefer integer based solutions. \$\endgroup\$ Commented Dec 9, 2023 at 21:10
5
\$\begingroup\$

Usually you can just test n < r*r instead which also has the benefit of being much cheaper than a sqrt. You do need to handle potential integer overflow for n > 264 - 2*232 + 1.

You can't hoist r*r out of the loop the way you can hoist sqrt(n) out but integer multiplication is cheap on any modern processor.

\$\endgroup\$
4
  • \$\begingroup\$ I think that's a really good idea, at least if the tests have reasonably low frequency. Trade a sloooow sqrt for a bunch of faster MULs. However, in the current case (sieve) the outer loop accounts for a significant part of the runtime, since most primes step less than a single time on average in a typical sieve window, and the number of prime factors is the same order of magnitude as the sieve window size (203,280,221 primes vs. sieve window between 1e6 and 1e9). In this case hoisting + sqrt is easier on impatient nerves... ;-) \$\endgroup\$
    – DarthGizka
    Commented Nov 6, 2014 at 17:25
  • \$\begingroup\$ However, if the loop contains the product anyway then the test is essentially free (since the compiler computes it only once). This doesn't answer the floaty question but it's a valuable optimisation tip that can often avoid the floaty question, up to a certain point. Unfortunately I can give you only one uptick... It's worth at least three. \$\endgroup\$
    – DarthGizka
    Commented Nov 6, 2014 at 17:33
  • \$\begingroup\$ Also, there may be a way to make your idea work even for high-frequency loops, by using strength reduction (n += (r << 1) + 1), especially since the term ((r << 1) + 1) occurs naturally in a sieve loop. New ideas, new avenues... Thanks! \$\endgroup\$
    – DarthGizka
    Commented Nov 6, 2014 at 17:41
  • \$\begingroup\$ The strength reduction idea is neat! On modern processors a branch is generally more expensive than a multiply, even a highly predictable branch like this, so I would expect the multiply to be effectively free anyway. \$\endgroup\$ Commented Nov 6, 2014 at 17:46
4
\$\begingroup\$

In order to get a fix on the issue I dug into scripture (C/C++ standards, math textbook) and ran extensive experiments.

The standards mandate that the result of sqrt(n) be within half an ulp of the true result. Since the subsequent rounding to uint32_t sheds 'excess' precision, the final result must be within one ulp of the true value (depending on rounding mode).

However, the square root is exact in this sense only for the double \$n' = n + h\$, which is the result of the initial conversion of the uint64_t n to a double n' with only 53 bits of precision.

The delta h is bounded by the ulp of n, which corresponds to \$(63 - 53 + 1) = 11\$ bits for values between \$2^{63}\$ and \$2^{64}-1\$. The estimate for the absolute error of \$\sqrt{n + h}\$ is \$\frac h {2 \sqrt n}\$; the ulp is constant over that range and equal to \$2^{11}\$. This maximises the absolute error for values close to the beginning of the range, e.g. \$n = 2^{63}\$. The value there is \$\frac{2^{11}}{(2 \sqrt{2^{63}})}\$ or \$\frac{2^{-21}}{\sqrt 2}\$. With other words, the error is negligible.

Since there are two precision-change roundings (casts) and one rounding computation involved, the absolute error for the final result can be greater than one even for a standards-conforming compiler operating in a fully standards-conforming mode.

When testing uncovered the existence of cases where my original max_factor32() diverged from the ideal result, I switched to the following 'luxury' version that should be exact under any and all circumstances, regardless of floating-strictness, optimisation settings or whatnot:

template<typename word_t>
word_t square (word_t n)  {  return word_t(n * n);  }

uint32_t max_factor32_b (uint64_t n)
{
   uint64_t r = uint64_t(std::sqrt(double(n)));

   while (r < UINT32_MAX && square(r) < n)  
      ++r;            
   while (r > UINT32_MAX || square(r) > n)  
      --r;

   return uint32_t(r);
}

That code can be optimised by hoisting a test out of a loop and similar operations, but I left it as is, to preserve the strange antisymmetry between the two cases.

Tests with trillions of computations in near the upper end of the range uncovered only 1025 cases where the two versions of the function returned different results, and in all these cases the simpler code reported a value that was higher by one than the ideal result (\$2^{32}-1\$ instead of \$2^{31}-2\$). For the Sieve of Eratosthenes that makes virtually no difference but other code might not be so forgiving, and other compilers might not be as well-behaved as gcc 4.8.1 and VC++ 2013...

The outlier results were for values of n between 0xFFFFFFFDFFFFFC00 and 0xFFFFFFFE00000000.

\$\endgroup\$
2
  • \$\begingroup\$ "The standards mandate that the result of sqrt(n) be within half an ulp of the true result." --> Please cite the standard and its requirement. I did not find such a requirement in the C nor C++ spec. \$\endgroup\$ Commented Dec 8, 2023 at 3:23
  • \$\begingroup\$ @chux-ReinstateMonica: I do not remember where I got that from and a cursory search of my standard stash did not turn up anything useful. The only thing I do remember is that the wording was on the lines 'the representable number closest to the true result'. Rounding modes can increase that to almost a full ULP, obviously. Pretty much all other sources - like gcc or MSVC docs - document 'within 1 ULP', probably in order to avoid talking about rounding and to leave room for optimisation. \$\endgroup\$
    – DarthGizka
    Commented Dec 22, 2023 at 14:41
3
\$\begingroup\$

I like DarthGizka's integer square root function with error correction, so I made a generic version of it. The code requires a C++14 compiler so It currently only works with clang.

template <typename T>
constexpr T sqrt_helper(T x, T lo, T hi)
{
  if (lo == hi)
    return lo;

  const T mid = (lo + hi + 1) / 2;
  if (x / mid < mid)
    return sqrt_helper<T>(x, lo, mid - 1);
  else
    return sqrt_helper(x, mid, hi);
}

template <typename T>
constexpr T ct_sqrt(T x)
{
  return sqrt_helper<T>(x, 0, x / 2 + 1);
}

template <typename T>
T isqrt(T x)
{
  T r = (T) std::sqrt((double) x);
  T sqrt_max = ct_sqrt(std::numeric_limits<T>::max());

  while (r < sqrt_max && r * r < x)
    r++;
  while (r > sqrt_max || r * r > x)
    r--;

  return r;
}

The above version uses two correction iterations for most numbers. As DarthGizka mentioned it is possible to optimize this, below is a version that uses no correction iteration for most numbers. Depending on your compiler and CPU this version can be faster.

template <typename T>
T isqrt(T x)
{
  T r = (T) std::sqrt((double) x);
  r = std::min(r, ct_sqrt(std::numeric_limits<T>::max()));

  while (r * r > x)
    r--;
  while (x - r * r > r * 2)
    r++;

  return r;
}
\$\endgroup\$
2
  • \$\begingroup\$ Very nice indeed! I took a long hard look at (lo + hi + 1) since expressions like that are very delicate, but it is perfectly all right for the parameters it gets from ct_sqrt(). I wish I had a compiler that accepts -std=c++14... \$\endgroup\$
    – DarthGizka
    Commented Dec 31, 2014 at 14:21
  • \$\begingroup\$ Have a look at my post on stackoverflow, it contains a C++11 version: stackoverflow.com/a/27709195/363778 \$\endgroup\$
    – user62153
    Commented Dec 31, 2014 at 14:34
1
\$\begingroup\$

If you have C99, I believe you can use sqrtl(), which will return a long double with the answer you need.

However, given your case, where all you really want is the integer you should run your sieve till, I'd recommend taking a look at the solutions provided here.

\$\endgroup\$
2
  • 2
    \$\begingroup\$ With many compilers, long double is just the same 64-bit float as double; the standard doesn't mandate that it actually be bigger (i.e. Intel's 80-bit float as in BC++, or a 128-bit float as on some workstation compilers). The topic you mentioned didn't have any solutions that weren't either 'use the sqrt(), Luke' or some slow iterative computation... I'd really like to crack the correctness thing for this particular example, so that I can use what I learnt in other, similar cases. \$\endgroup\$
    – DarthGizka
    Commented Nov 6, 2014 at 11:56
  • \$\begingroup\$ Even if a long double has sufficient bits to represent a 64-bit integer exactly, the rounding issue could still occur. Unless I'm mistaken, even a hypothetical sqrtlllll(UINT64_MAX) would still return 2^32 instead of 2^32-1. \$\endgroup\$
    – DarthGizka
    Commented Nov 6, 2014 at 14:32
1
\$\begingroup\$

An alterative to risking floating point issues is to perform an integer based square root:

C based code

uint32_t usqrt64(uint64_t n) {
  if (n < 2) {
    return (uint32_t)n;
  }
  uint32_t s = usqrt64(n >> 2) << 1;
  uint32_t l = s + 1;
  if ((uint64_t)l * l > n) {
    return s;
  }
  return l;
}

Or perhaps a non-recursive version

uint32_t usqrt64(uint64_t t) {
  uint64_t s, b;
  for (b = 0, s = t; s;  b++, s >>= 1) {
    ;
  }

  s = (uint64_t)1 << (b >> 1);
  if (b & 1) {
    s += s >> 1;
  }

  do {
    b = t / s;
    s = (s + b) >> 1;
  } while (b < s);

  return (uint32_t)s;
}

or (adapted from Integer_square_root)

// Binary search
uint32_t usqrt(uint64_t y) {
    if (y == UINT64_MAX) {
      return UINT32_MAX;
    }
    uint64_t L = 0;
    uint64_t R = y + 1;

    while (L + 1 != R) {
      uint64_t M = (L + R) / 2;
      if (M * M <= y) {
        L = M;
      } else {
        R = M;
      }
    }

    return (uint32_t) L;
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.