slow nextafter

[mglisse]

nextafter is used in FPU.h to work around the lack of directed rounding, but it is very slow. Here is a dump with some ideas for replacing it, and some tests to check that it is not too broken. Others are mentioned in #5510, in particular things like multiplying with 1±eps as a double.
Note that the optimization target should be wasm, the performance of this code compiled with gcc on x86_64 may not be relevant. Browsers are much worse at generating efficient code than the C++ compilers we are used to. And I don't know if, for standard functions like nextafter, the browser has a way to call the native function (in which case it may be hard to beat) or if it has to call a wasm version of nextafter, in which case the following code should help.
In some cases in Interval_nt.h, we may have extra information (already know the sign), but taking advantage of it does not seem worth the trouble.

#include <math.h>
#include <assert.h>
#include <string.h>
#include <limits>
#include <vector>
// special case of nextafter when y=+inf and x cannot be NaN
double nextdouble(double x)
{
  //CGAL_assertion(!std::isnan(x));
  assert(!std::isnan(x));
  int64_t X;
  // C++20: use std::bit_cast instead of memcpy
  memcpy(&X, &x, 8);
#if 0
  if(X==0x7ff0000000000000) return x; // x = +inf
  if(X==0x8000000000000000) return std::numeric_limits<double>::denorm_min(); // x = -0
  // if(X>=0) ++X; else --X;
  X += (X<0) ? -1 : 1; // (X >> 63) | 1
#else
  if(X>=0) {
    if(X==0x7ff0000000000000) return x; // x = +inf
    ++X;
    // How about computing ++X, read as double, then minsd(x,+inf), which if x is nan returns inf?
  } else {
    if(X==0x8000000000000000) return std::numeric_limits<double>::denorm_min(); // x = -0
    --X;
    // We could check the borrow to detect the special value.
    // It is probably necessary to handle negative 0, but it is worth making sure.
  }
#endif
  memcpy(&x, &X, 8);
  return x;
}

static void check(double d){
  double a = nextafter(d, std::numeric_limits<double>::infinity());
  double b = nextdouble(d);
  assert(memcmp(&a,&b,8)==0);
}
int main(){
  for(int s=-1;s<2;s+=2){
    for(double D : { 0., std::numeric_limits<double>::denorm_min(), std::numeric_limits<double>::min(), 1., std::numeric_limits<double>::max(), std::numeric_limits<double>::infinity() }){
      double d = s * D;
      check(d);
      check(nextafter(d, -std::numeric_limits<double>::infinity()));
      check(nextafter(d, std::numeric_limits<double>::infinity()));
    }
  }
}