Does the average of floats satisfy conditions for machine-precision bisection?

user_11130854 0
Viewed 59

Suppose you want to run a bisection algorithm/binary search down to machine precision, which will always terminate in a few hundred steps due to exponential halving of the range of floats. For this to work, the following conditions should be satisfied: Let floats A<B and M = (A+B)/2. Then

  1. A < M < B if and only if A and B are not neighboring floats.

  2. A=M or M=B if and only if A and B are neighboring floats.

Is this always guaranteed in floating point arithmetic?

If not, is there any reasonable definition of a midpoint M for which these conditions hold? (Obviously, defining M as the upper neighbor of A would work, but that would not be a reasonable definition.)

Edit 1: As pointed out in the comments, the sum of A+B may overflow. I am not necessarily asking about this specific sequence of operations, something like M = A/2 + B/2 would also be valid, as well as other midpoint methods.

Edit 2: https://scicomp.stackexchange.com/questions/20369/robust-computation-of-the-mean-of-two-numbers-in-floating-point/20379 is related, for the weaker condition that min(A,B) <= M <= max(A,B), which works for M=A+(B/2−A/2) and also does not overflow. Another interesting approach for bisection is based on reinterpreting floats as integers: https://www.juliabloggers.com/bisecting-floating-point-numbers/

floating-accuracy numeric precision floating-point
1 Answers

This is an implementation of a midpoint method which actually does not use the average, but rather splits the set of floats between A and B at their median position. This is based on https://www.juliabloggers.com/bisecting-floating-point-numbers/, and I believe it satisfies the requirements. Code review is appreciated.

#include <cstdint> 
#include <math.h> 
#include <cstring>
#include <cassert>


double bisect_float(double L, double R) {
        
    // Use the usual float rules for combining non-finite numbers
    if(!isfinite(L) || !isfinite(R)){
        return(L+R);
    }
    
    if (L==R){return(L);}
    

    // Always return 0.0 when inputs have opposite sign and neither is zero
    if ( (L>0.0 && R<0.0) || (L<0.0 && R>0.0) )   {
        return(0.0);
    }
    
    // check whether we need to change the sign
    bool negate = (L < 0.0 || R < 0.0);
        
    // reinterpret floats as integers
    // this gives their ordering index
    uint64_t Li;
    uint64_t Ri;
    assert(sizeof(uint64_t) == sizeof L);
    assert(sizeof(uint64_t) == sizeof R);
    L = abs(L); // remove sign bit
    R = abs(R); // remove sign bit
    memcpy(&Li, &L, sizeof Li);
    memcpy(&Ri, &R, sizeof Ri);
    
    // find the median float representation
    uint64_t res_i;
    if (Li>Ri){
        std::swap(Li, Ri);
    }
    res_i = Li + ((Ri-Li)>>1);
    //  NOTE: (Li>>1) + (Ri>>1) would be wrong for (3,5), as it would be 1+2=3 and hence not the midpoint
    
    // reinterpret median index as float
    double res;
    memcpy(&res, &res_i, sizeof res);
    
    if (negate){
        return(-res);
    } else {
        return(res);
    }   
}
user_111308540
Related