8

I stumbled onto this code here.

Generators doubleSquares(int value)
{
    Generators result;
    for (int i = 0; i <= std::sqrt(value / 2); i++) // 1
    {
        double j = std::sqrt(value - std::pow(i, 2)); // 2
        if (std::floor(j) == j) // 3
            result.push_back( { i, static_cast<int>(j) } ); // 4
    }
    return result;
}

Am I wrong to think that //3 is dangerous ?

ApproachingDarknessFish
  • 14,133
  • 7
  • 40
  • 79
NoSenseEtAl
  • 28,205
  • 28
  • 128
  • 277
  • 2
    I think it's safe -- `std::floor` returns an integer value (though in a `float` or `double`) and all integers can be represented without the slight inaccuracies that can affect fractional values. – Joe Jan 21 '14 at 01:57
  • 1
    @Joe: At some point (values very far from zero), floating point representation will not necessarily represent exact integers, either. That's unlikely to be seen in practice here, and even so I'm not sure (anyone? anyone?) whether the logic is still sound because the sqrt would demonstrate the same variance? – Ben Zotto Jan 21 '14 at 02:01
  • @Joe What if `j` itself has a fractional value? EG, if `sqrt(16)` returns 4.000000000000001 and the comparison returns false? – ApproachingDarknessFish Jan 21 '14 at 02:01
  • "Any integer with absolute value less than 2^53 can be exactly represented in the double precision format" (http://en.wikipedia.org/wiki/Floating_point) – Ben Zotto Jan 21 '14 at 02:03

3 Answers3

10

This code is not guaranteed by the C++ standard to work as desired.

Some low-quality math libraries do not return correctly rounded values for pow, even when the inputs have integer values and the mathematical result can be exactly represented. sqrt may also return an inaccurate value, although this function is easier to implement and so less commonly suffers from defects.

Thus, it is not guaranteed that j is exactly an integer when you might expect it to be.

In a good-quality math library, pow and sqrt will always return correct results (zero error) when the mathematical result is exactly representable. If you have a good-quality C++ implementation, this code should work as desired, up to the limits of the integer and floating-point types used.

Improving the Code

This code has no reason to use pow; std::pow(i, 2) should be i*i. This results in exact arithmetic (up to the point of integer overflow) and completely avoids the question of whether pow is correct.

Eliminating pow leaves just sqrt. If we know the implementation returns correct values, we can accept the use of sqrt. If not, we can use this instead:

for (int i = 0; i*i <= value/2; ++i)
{
    int j = std::round(std::sqrt(value - i*i));
    if (j*j + i*i == value)
        result.push_back( { i, j } );
}

This code only relies on sqrt to return a result accurate within .5, which even a low-quality sqrt implementation should provide for reasonable input values.

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
  • I think you can rely on `sqrt` returning the proper value for IEEE floats, but this doesn't apply to `pow`. – user541686 Jan 21 '14 at 02:06
  • 1
    @Mehrdad: Why do you think that? – Eric Postpischil Jan 21 '14 at 02:07
  • Because [*"The IEEE standard goes further than just requiring the use of a guard digit. It gives an algorithm for addition, subtraction, multiplication, division and **square root**, and requires that implementations produce the same result as that algorithm."*](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) – user541686 Jan 21 '14 at 02:08
  • 4
    @Mehrdad: Yes, if a square root implementation conforms to IEEE 754, then it returns a correctly rounded result. But that is not the question. The C++ `std::sqrt` function **is not bound to IEEE 754**. A C++ implementation might choose to conform to IEEE 754 or it might choose not to conform. So, why do you think `std::sqrt` is C++ returns a correctly rounded value? – Eric Postpischil Jan 21 '14 at 02:10
  • I *never said* `sqrt` returns a correctly rounded value. I said it returns a correctly rounded value **for IEEE floats**. You asked me why, and I linked you to the reason. Now you're incorrectly generalizing my statement. – user541686 Jan 21 '14 at 02:11
  • 3
    @Mehrdad: A C++ implementation may use IEEE-754 floating-point values and formats but not use IEEE-754 operations and requirements. So saying “for IEEE floats” is insufficient to say that `std::sqrt` is bounded to IEEE 754. In fact, many C++ implementations do use IEEE-754 values and formats but do not conform to IEEE-754 requirements for correctness. – Eric Postpischil Jan 21 '14 at 02:12
  • 1
    The goal of my comment was to get a point across (i.e. that `sqrt` and `pow` are fundamentally very different in terms of required accuracy), not to pedantically nitpick over the word choice. Seems like you understood the point so I'm done. – user541686 Jan 21 '14 at 02:18
1

There are two different, but related, questions:

  1. Is j an integer?
  2. Is j likely to be the result of a double calculation whose exact result would be an integer?

The quoted code asks the first question. It is not correct for asking the second question. More context would be needed to be certain which question should be being asked.

If the second question should be being asked, you cannot depend only on floor. Consider a double that is greater than 2.99999999999 but less than 3. It could be the result of a calculation whose exact value would be 3. Its floor is 2, and it is greater than its floor by almost 1. You would need to compare for being close to the result of std:round instead.

Patricia Shanahan
  • 25,849
  • 4
  • 38
  • 75
0

I would say it is dangerous. One should always test for "equality" of floating point numbers by comparing the difference between the two numbers with an acceptably small number, e.g.:

#include <math.h>
...
if (fabs(std::floor(j) - j) < eps) {
...

... where eps is a number that is acceptably small for one's purpose. This approach is essential unless one can guarantee that the operations return exact results, which may be true for some cases (e.g. IEEE-754-compliant systems) but the C++ standard does not require that this be true. See, for instance Cross-Platform Issues With Floating-Point Arithmetics in C++.

Simon
  • 10,679
  • 1
  • 30
  • 44
  • In that case, better replace `std::floor` with something else. – Markus Mayr Jan 21 '14 at 12:02
  • @MarkusMayr: Good catch. The inexact result could be slightly less than the integer, which would lead to `floor()` giving the next-lower integer instead. It would be better to use a `round()` function, either from C++11 or one of the options in the answers to [round() for float in C++](http://stackoverflow.com/questions/485525/round-for-float-in-c). – Simon Jan 21 '14 at 21:18