1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137

#include <iostream>
#include <type_traits>
#include <tuple>
#include <cmath>
template <typename Float>
struct float_info{};
template <>
struct float_info<float>{
typedef unsigned integer_type;
static const unsigned bits = 32;
static const unsigned exponent_size = 8;
static const unsigned exponent_mask = (1 << exponent_size)  1;
static const unsigned mantissa_size = 23;
static const integer_type mantissa_mask = ((integer_type)1 << mantissa_size)  1;
static const unsigned exponent_baseline = (1 << (exponent_size  1))  2;
};
template <>
struct float_info<double>{
typedef unsigned long long integer_type;
static const unsigned bits = 64;
static const unsigned exponent_size = 11;
static const unsigned exponent_mask = (1 << exponent_size)  1;
static const unsigned mantissa_size = 52;
static const integer_type mantissa_mask = ((integer_type)1 << mantissa_size)  1;
static const unsigned exponent_baseline = (1 << (exponent_size  1))  2;
};
template <typename T>
typename std::enable_if<std::is_floating_point<T>::value, bool>::type
about_equal(T a, T b, int digits){
return false;
}
template <typename T>
class FloatOperations{
typedef float_info<T> fi;
typedef typename fi::integer_type it;
static void normalize(T &l, T &r){
//If l and r have different signs, make l the one with the
//greater absolute value and flip both their signs until l is
//also positive, then make r equal to l + abs(l  r),
//maintaining the delta but ensuring that r is positive.
auto ls = l < 0;
auto rs = r < 0;
if (ls != rs){
if (ls){
l = l;
r = r;
}
if (l < r)
std::swap(l, r);
r = l + l  r;
}
}
public:
static_assert(std::numeric_limits<T>::is_iec559, "We need standard floats.");
static it to_int(T x){
it ret;
memcpy(&ret, &x, sizeof(x));
return ret;
}
static T to_float(it x){
T ret;
memcpy(&ret, &x, sizeof(x));
return ret;
}
static T zero_exponent(T x){
auto [sign, exponent, mantissa] = decompose(x);
return recompose(sign, fi::exponent_baseline, mantissa);
}
static auto decompose(T x){
auto y = to_int(x);
unsigned sign = y >> (fi::bits  1);
int exponent = (y >> (fi::bits  1  fi::exponent_size)) & fi::exponent_mask;
it mantissa = y & fi::mantissa_mask;
return std::make_tuple(sign, exponent, mantissa);
}
static auto recompose(unsigned sign, int exponent, it mantissa){
it ret = sign;
ret <<= fi::exponent_size;
ret = exponent & fi::exponent_mask;
ret <<= fi::mantissa_size;
ret = mantissa & fi::mantissa_mask;
return to_float(ret);
}
static bool approx_equals(T l, T r, unsigned digits){
//log2(10) = 3.3219... vvvvvvvvv to get ceil()
digits = (digits * 33219 + (10000  1)) / 10000;
normalize(l, r);
auto [l_sign, l_exponent, l_mantissa] = decompose(l);
auto [r_sign, r_exponent, r_mantissa] = decompose(r);
if (l_exponent != r_exponent)
return !digits;
l_mantissa >>= fi::mantissa_size  digits;
r_mantissa >>= fi::mantissa_size  digits;
return l_mantissa == r_mantissa;
}
static bool approx_equals_no_black_magic(T l, T r, unsigned digits){
//log2(10) = 3.3219... vvvvvvvvv to get ceil()
digits = (digits * 33219 + (10000  1)) / 10000;
normalize(l, r);
if (l < 0){
l = l;
r = r;
}
const auto log2 = 1.0 / log(2);
l = l * pow(2, floor((digits  log(l)) * log2));
r = r * pow(2, floor((digits  log(r)) * log2));
return floor(l) == floor(r);
}
};
const int digits_of_precision = 5;
const double epsilon = pow(10.0, (double)digits_of_precision);
bool usual_approx_equals(double a, double b){
return abs(a  b) < epsilon;
}
int main(){
typedef FloatOperations<double> fo;
double big_divisor = pow(2.0, 32.0);
double l1 = 3.1415926535897932384626433832795 / big_divisor;
double r1 = 2.7182818284590452353602874713527 / big_divisor;
double l2 = 1125899906842624.0;
double r2 = 1125899906842625.0;
std::cout << "False positive: usual_approx_equals(" << l1 << ", " << r1 << ") = " << usual_approx_equals(l1, r1) << std::endl;
std::cout << "False negative: usual_approx_equals(" << l2 << ", " << r2 << ") = " << usual_approx_equals(l2, r2) << std::endl;
std::cout << "True negative approx_equals(" << l1 << ", " << r1 << ") = " << fo::approx_equals(l1, r1, digits_of_precision) << std::endl;
std::cout << "True positive approx_equals(" << l2 << ", " << r2 << ") = " << fo::approx_equals(l2, r2, digits_of_precision) << std::endl;
std::cout << "True negative approx_equals_no_black_magic(" << l1 << ", " << r1 << ") = " << fo::approx_equals_no_black_magic(l1, r1, digits_of_precision) << std::endl;
std::cout << "True positive approx_equals_no_black_magic(" << l2 << ", " << r2 << ") = " << fo::approx_equals_no_black_magic(l2, r2, digits_of_precision) << std::endl;
return 0;
}
