GitLab

Update: I see now why we indeed need this function. Let me go through my thought process anyway.

The issue mentioned in bisection.hh regarding unpredictable rounding from 80 to 64 bits is worrisome. It's also mentioned here:

https://gcc.gnu.org/wiki/x87note

However, from my understanding (and please correct me if I'm wrong), this only happens if instructions for x87 FPUs are generated, and in particular not if SSE is enabled. With SSE, rounding is reasonable and predictable.

So now there are two questions:

1. when do we generate code for x87
2. can we really trigger an issue such as this when we do?

Question no. 1 appears to depend considerably on the compiler. It's actually hard to do this with clang. I need to use the command line

clang++ -m32 -mno-sse

The flat -mfpmath=387 cannot be used on its own and something like -m32 -mtune=generic did not suffice: SSE was still used.

With gcc, this is quite different. As soon as gcc is in 32bit mode, it will use x87 FP arithmetic. To keep it from doing so, we need to use

g++ -m32 -msse2 -mfpmath=sse

The former flag is implied by march for architectures that support sse2. But the latter is not.

To test the implications of this issue, I wrote two test programs. The first one bisects the interval [-M_PI_2, M_PI_2] in search of 0, i.e. the zero of std::sin:

#include <cassert>
#include <iostream>
#include <random>

int main() {
  size_t const num_samples = 100;
  std::default_random_engine generator;
  std::uniform_real_distribution<double> left_distribution(-M_PI_2, 0);
  std::uniform_real_distribution<double> right_distribution(0, M_PI_2);
  size_t conv_counter = 0;
  for (size_t i = 0; i < num_samples; ++i) {
    double left = left_distribution(generator);
    double right = right_distribution(generator);
    // we need approximately 1000 steps to go down to 1e-300
    for (size_t j = 0; j < 1500; ++j) {
      double const middle = 0.5 * left + 0.5 * right;
      if (middle == left or middle == right) {
        assert(std::abs(middle) < 1e-14);
        conv_counter++;
        break;
      }
      double const mvalue = std::sin(middle);
      if (mvalue > 0)
        right = middle;
      else
        left = middle;
      // std::cout << "[" << left << ", " << right << "]" << std::endl;
    }
  }
  std::cout << conv_counter << std::endl;
  return conv_counter == num_samples ? 0 : 1;
}

The second program checks if the the midpoint between two consecutive floating point numbers really needs to be equal to at least one of them.

#include <cmath>
#include <iostream>

int
main()
{
  double x = 1.0;
  double y = std::nextafter(x, 2.0);
  double middle = (x+y)/2.0;
  //std::cout << middle << std::endl;
  if (x == middle or y == middle)
    std::cout << "Behaving as expected" << std::endl;
  else
    std::cout << "Help, where are we?" << std::endl;
}

As I've mentioned before: To see proper behaviour, you can use something like

clang++ -std=c++11 -m32 foo.cc

or

g++ -std=c++11 -m32 -msse2 -mfpmath=sse foo.cc

For the opposite, we need

clang++ -std=c++11 -m32 -mno-sse foo.cc

or

g++ -std=c++11 -m32 -O1 foo.cc

(note the need for optimisation, at least -O1).