CGAL ApproxDouble Kernel: pushing naive inexact operations to the limits

fuzzy_double.h

// https://gist.github.com/ochafik/84e4e2455cd8a37b3cf10e0ebfaf7d2
// Copyright 2021 Google LLC.
// SPDX-License-Identifier: Apache-2.0
//
// Context: "Can Nef_polyhedron_3 (and minkowski_sum_3 ) be used with inexact kernels?"
// https://github.com/CGAL/cgal/issues/5490c
#pragma once

#include <cmath>
#include <type_traits>

class FuzzyDouble {
public:
  static constexpr double epsilon = 0.00000001;
  double value;

  template <typename T, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
  FuzzyDouble(T value) : value(static_cast<double>(value)) {}
  FuzzyDouble(const FuzzyDouble& other) { *this = other; }
  FuzzyDouble() : value(0.0) {}

  operator double() const { return value; }
  FuzzyDouble &operator=(const FuzzyDouble& other) { value = other.value; return *this; }
  FuzzyDouble operator-() const { return -value; }

#define OP_TEMPLATE(T) \
  template <typename T, std::enable_if_t<(std::is_arithmetic<T>::value || std::is_assignable<FuzzyDouble, T>::value), bool> = true>

  OP_TEMPLATE(T) bool operator==(const T& x) const { return std::abs(value - static_cast<double>(x)) <= FuzzyDouble::epsilon;}
  OP_TEMPLATE(T) bool operator!=(const T& x) const { return !(*this == x); }
  OP_TEMPLATE(T) bool operator<(const T& x) const { return value + FuzzyDouble::epsilon < static_cast<double>(x); }
  OP_TEMPLATE(T) bool operator<=(const T& x) const { return value <= static_cast<double>(x) + FuzzyDouble::epsilon; }
  OP_TEMPLATE(T) bool operator>(const T& x) const { return value > static_cast<double>(x) + FuzzyDouble::epsilon; }
  OP_TEMPLATE(T) bool operator>=(const T& x) const { return value + FuzzyDouble::epsilon >= static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator*(const T& x) const { return value * static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator/(const T& x) const { return value / static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator+(const T& x) const { return value + static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble operator-(const T& x) const { return value - static_cast<double>(x); }
  OP_TEMPLATE(T) FuzzyDouble& operator*=(const T& x) { value *= static_cast<double>(x); return *this; }
  OP_TEMPLATE(T) FuzzyDouble& operator/=(const T& x) { value /= static_cast<double>(x); return *this; }
  OP_TEMPLATE(T) FuzzyDouble& operator+=(const T& x) { value += static_cast<double>(x); return *this; }
  OP_TEMPLATE(T) FuzzyDouble& operator-=(const T& x) { value -= static_cast<double>(x); return *this; }
};

std::ostream& operator<<(std::ostream& out, const FuzzyDouble& n) { out << n.value; return out; }
std::istream& operator>>(std::istream& in, FuzzyDouble& n) { in >> n.value; return in; }

namespace std { FuzzyDouble sqrt(const FuzzyDouble& x) { return std::sqrt(x.value); } }

namespace CGAL {
  inline FuzzyDouble min BOOST_PREVENT_MACRO_SUBSTITUTION(const FuzzyDouble& x,const FuzzyDouble& y){ return std::min(x.value, y.value); }
  inline FuzzyDouble max BOOST_PREVENT_MACRO_SUBSTITUTION(const FuzzyDouble& x,const FuzzyDouble& y){ return std::max(x.value, y.value); }

  CGAL_DEFINE_COERCION_TRAITS_FOR_SELF(FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(short    ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(int      ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(long     ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(float    ,FuzzyDouble)
  CGAL_DEFINE_COERCION_TRAITS_FROM_TO(double   ,FuzzyDouble)

  template <>
  class Algebraic_structure_traits<FuzzyDouble> : public Algebraic_structure_traits_base<FuzzyDouble,Field_with_sqrt_tag> {
  public:
    typedef Tag_false       Is_exact;
    typedef Tag_true        Is_numerical_sensitive;
  };

  template <>
  class Real_embeddable_traits<FuzzyDouble>: public INTERN_RET::Real_embeddable_traits_base<FuzzyDouble,CGAL::Tag_true>
  {
    typedef Algebraic_structure_traits<Type>        AST;
  public:

    typedef Tag_true                                Is_real_embeddable;
    typedef Uncertain<bool>                         Boolean;
    typedef Uncertain<CGAL::Comparison_result>      Comparison_result;
    typedef Uncertain<CGAL::Sign>                   Sign;

    typedef AST::Is_zero    Is_zero;

    struct Is_finite: public CGAL::cpp98::unary_function<Type,Boolean>{
      inline Boolean operator()(const Type &x) const { return !std::isinf(x.value) && !std::isnan(x.value); };
    };

    struct Abs: public CGAL::cpp98::unary_function<Type,Type>{
      inline Type operator()(const Type &x) const { return std::abs(x.value); };
    };

    struct To_interval: public CGAL::cpp98::unary_function<Type,std::pair<double,double> >{
      inline std::pair<double,double> operator()(const Type &x)const{
        return std::make_pair(x.value - FuzzyDouble::epsilon, x.value + FuzzyDouble::epsilon);
      };
    };
  };

} //namespace CGAL