Issues with the sqrt() function in C++

Question

So, I was writing code in C++, which required an intermediate step to check whether a number was a perfect square. I wrote the following code.

int sqrt_of_t = (int)sqrt(t);
if (sqrt_of_t*sqrt_of_t != t)
{
    cout << "NO" << endl;
}

This code gives the correct results in my system, but it fails when passing it through an online judge in Codeforces. The case where it fails doesn't have any overflow associated with it or anything (really small test cases). So, can anyone explain where it went wrong and suggest some alternate method to check if a number is a perfect square or not, which will work on all systems and not show behaviors like this. Here t is an int too.

Welcome to the impossibility of perfectly representing all floating point numbers. — user4581301, Dec 25 '22 at 06:49
[Related issue](https://stackoverflow.com/questions/25678481/why-does-pown-2-return-24-when-n-5-with-my-compiler-and-os/25678721#25678721). Do not use floating point functions for integer-based problems. Your cast to `int` makes the potential problem worse. — PaulMcKenzie, Dec 25 '22 at 06:49
General case: [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) — user4581301, Dec 25 '22 at 06:50
@PaulMcKenzie Is there anyway to solve this? Instead of say do a binary search to find the square root (works only for ints, which will satisfy my use case) — Vedanta Mohapatra, Dec 25 '22 at 06:52
@TomSolid Are you sure about that? `sqrt(1e301)` is an integer, but I'm pretty sure `1e301` is not a square number. — john, Dec 25 '22 at 06:52
@VedantaMohapatra -- The website "codeforces" is one that asks random puzzle questions. This particular question was more than likely designed so that code such as yours will fail. The next thing for you to do is to do your own research into figuring out how to get exact square roots without using floating point. — PaulMcKenzie, Dec 25 '22 at 06:54
@kotatsuyaki Yes it is, it's roughly equal to `3e150`. At that magnitude all floating point numbers are integers. — john, Dec 25 '22 at 07:25
Does this answer your question? [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) — Fedor, Dec 25 '22 at 07:36
Personally, if the requirement is to check if an integral value was the square of another integral value, I wouldn't use floating point at all. It is easy enough to find algorithms that use only integral operations to compute an "integer square root", defined as something like "the integral square root of a positive integral value `x` is the largest positive integral value that, when multiplied by itself, gives a result less than or equal to `x`". For example, the integral square root of `26` is `5` (not `6`) but `5*5` is not equal to `26`, so `26` is not a perfect square. — Peter, Dec 25 '22 at 08:46
"The case where it fails doesn't have any overflow associated with it or anything". How do you know? — n. m. could be an AI, Dec 27 '22 at 19:46
@n.m. The test cases were pretty small. There didn't seem to be any calculations which might have overflowed an int. — Vedanta Mohapatra, Dec 27 '22 at 19:51

John Zwinck · Answer 1 · 2022-12-25T07:40:33.250

5

sqrt() returns a floating point number which you cast to int, which truncates any fractional part. The problem is that floating point cannot represent all integers exactly, so you may end up with something like 19.99999999999999 which you expect to be 20 but is actually 19 when cast to an integer.

To fix it, use rounding instead:

long sqrt_of_t = lrint(sqrt(t));

edited Dec 25 '22 at 07:40

answered Dec 25 '22 at 06:55

John Zwinck

239,568
38
324
436

Well, floating point can represent integers correctly as long as as distance of highest and lowest one-bits is not larger than the mantissa – and I'd assume this is the case for all input relevant for the test. However `sqrt` itself might lose precision during calculations or on intermediate results. – Aconcagua Dec 25 '22 at 07:35
@Aconcagua: I tested it and I believe it works correctly. Do you have an example `t` where it fails, where `t` is representable in a `long long`? – John Zwinck Dec 25 '22 at 07:37
@Aconcagua: If I understand your line of thought properly, the reason you can't find an error is because of the second step `sqrt_of_t*sqrt_of_t != t`. In the first step we find an integer which is near to the square root of the input, and even if the input is not perfectly representable in a `double`, the maximum error is not so much as to produce a wrong guess for the integer square root with my rounding solution. The fact that we check the guess with `sqrt_of_t*sqrt_of_t != t` then removes all doubt: the guess is validated in integer space with no inaccuracy. – John Zwinck Dec 26 '22 at 07:45
1

Oh, just noticing (thought error of mine before): False positives we indeed cannot get; if we discover `root*root == value` then we have found the root already, of course... Remains the case of false negatives. Obviously, if some of the up to 63 one-bits do not fit into the mantissa then the least significant bits get cut away. These are at most 10 digits, though, in contrast to 53 yet remaining – just as obviously a ratio far to small to provoke a *significant* rounding error for the root... So now seeing, this occurs, the result still will indeed be correct. – Aconcagua Dec 27 '22 at 10:48
2

`std::sqrt` is required by the IEEE standard to be correctly rounded from the infinitely precise result. In particular, the exact result is produced if it can be represented in the floating-point type. If `sqrt(400)` returns something like 19.99999999999999 then it is broken. – n. m. could be an AI Dec 27 '22 at 19:50

selbie · Answer 2 · 2022-12-26T02:30:37.950

4

sqrt, on many systems returns an approximation.

For example, sqrt(25) might return something like 4.99999999.

Hence, 4.99999999 * 4.99999999 is slightly less than 25.

My advice would be to do a binary search across the number space to see if the number is a perfect square. Avoid floating point whenever you need precise results.

bool isPerfectSquare(long long t)
{
    bool result = false;
    if ((t == 0) || (t == 1)) {
        return true;
    }
    if (t < 0) {
        return false;
    }

    long long low = 1;
    long long high = t / 2;

    while (low < high)
    {
        auto mid = (high + low) / 2;
        auto sq = mid * mid;
        if (sq == t) {
            result = true;
            break;
        }
        if (sq < t) {
            low = mid + 1;
        }
        else {
            high = mid - 1;
        }
    }
    return result;
}

edited Dec 26 '22 at 02:30

answered Dec 25 '22 at 06:55

selbie

100,020
15
103
173

I was also thinking of doing this similarly. Thanks – Vedanta Mohapatra Dec 25 '22 at 06:56
Not an "estimate," but an approximation. – DYZ Dec 25 '22 at 07:00
1

This is very inefficient, see https://godbolt.org/z/7Ph4qsxnq – John Zwinck Dec 25 '22 at 07:01
1

@JohnZwinck - an "inefficient" solution is better than a solution that doesn't work. I'm sure there's something in my college calculus book for computing the square root in faster time, but I'll stand by my answer. – selbie Dec 26 '22 at 00:39
2

std::sqrt is required by the IEEE standard to be correctly rounded from the infinitely precise result. In particular, the exact result is produced if it can be represented in the floating-point type. If sqrt(25) returns anything but exactly 5, it is broken. – n. m. could be an AI Dec 27 '22 at 19:53

Gene · Answer 3 · 2022-12-27T19:13:42.510

Here is Knuth's very interesting algorithm for computing integer square roots with only shift and add. It rounds down for non-square inputs.

uint32_t isqrt1(uint32_t x) {
  uint32_t r = 0, r2 = 0; 
  for (int p = 15; p >= 0; --p) {
    uint32_t tr2 = r2 + (r << (p + 1)) + (1u << (p << 1));
    if (tr2 <= x) {
      r2 = tr2;
      r |= (1u << p);
    }
  }
  return r;
}

This works by trying to set each bit to 1, high to low, maintaining the square of the prospective root computed so far. Each bit is "or"ed into the result if doing so produces a square no greater than the input value. It can be modified to detect the case where the prospect is an exact square.

bool is_exact_square(uint32_t x) {
  if (x == 0) return true;
  uint32_t r = 0, r2 = 0; 
  for (int p = 15; p >= 0; --p) {
    uint32_t tr2 = r2 + (r << (p + 1)) + (1u << (p << 1));
    if (tr2 == x) return true;
    if (tr2 < x) {
      r2 = tr2;
      r |= (1u << p);
    }
  }
  return false;
}

I'm adding the for general interest. The binary search suggestion is good. Maybe better unless you're working on a machine without fast multiply.

Issues with the sqrt() function in C++

3 Answers3