real_polyhedral_homotopy_continuation/reference/include/num_traits.hpp

// This file is part of Bertini 2.
//
// num_traits.hpp is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
// num_traits.hpp is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with num_traits.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
//  Copyright(C) 2015 - 2021 by Bertini2 Development Team
//
//  See <http://www.gnu.org/licenses/> for a copy of the license,
//  as well as COPYING.  Bertini2 is provided with permitted
//  additional terms in the b2/licenses/ directory.

// individual authors of this file include:
// silviana amethyst, university of wisconsin eau claire

/**
\file num_traits.hpp

\brief Provides an Eigen-like NumTraits struct for querying traits of a number
type.

The bertini::NumTraits struct provides NumDigits and NumFuzzyDigits functions.
*/

#pragma once

#include <cmath>

#include "double_extensions.hpp"

namespace bertini {
template <typename T>
T random_unit();

template <typename T>
struct NumTraits {};

template <>
struct NumTraits<double> {
  inline static unsigned num_digits() { return 16; }

  inline static unsigned num_fuzzy_digits() { return 14; }

  inline static unsigned tolerance_to_digits(double tol) {
    return ceil(-log10(tol));
  }

  inline static double from_rational(mpq_rational const& n,
                                     unsigned /* precision */) {
    return double(n);
  }

  using Real = double;
  using Complex = dbl_complex;
};

template <>
struct NumTraits<dbl_complex> {
  inline static unsigned num_digits() { return 16; }

  inline static unsigned num_fuzzy_digits() { return 14; }

  inline static dbl_complex from_rational(mpq_rational const& n,
                                          unsigned /* precision */) {
    return dbl_complex(static_cast<double>(n), 0);
  }

  using Real = double;
  using Complex = dbl_complex;
};

inline unsigned precision_increment() { return 10; }

inline unsigned double_precision() { return 16; }

inline unsigned lowest_multiple_precision() { return 20; }

inline unsigned max_precision_allowed() { return 1000; }

/**
\brief Get the precision of a number.

For doubles, this is trivially 16.
*/
inline unsigned precision(double) { return double_precision(); }

/**
For complex doubles, throw if the requested precision is not double_precision.
*/
inline void precision(double, unsigned prec) {
  if (prec != double_precision()) {
    std::stringstream err_msg;
    err_msg << "trying to change precision of a double to " << prec;
    throw std::runtime_error(err_msg.str());
  }
}

/**
\brief Get the precision of a number.

For complex doubles, this is trivially 16.
*/
inline unsigned precision(dbl_complex) { return double_precision(); }

/**
For complex doubles, throw if the requested precision is not double_precision.
*/
inline void precision(dbl_complex, unsigned prec) {
  if (prec != double_precision()) {
    std::stringstream err_msg;
    err_msg << "trying to change precision of a double to " << prec;
    throw std::runtime_error(err_msg.str());
  }
}

inline dbl_complex rand_complex() {
  using std::abs;
  using std::sqrt;
  static std::default_random_engine generator;
  static std::uniform_real_distribution<double> distribution(-1.0, 1.0);
  dbl_complex returnme(distribution(generator), distribution(generator));
  return returnme / sqrt(abs(returnme));
}

template <>
inline dbl_complex random_unit<dbl_complex>() {
  static std::default_random_engine generator;
  static std::uniform_real_distribution<double> distribution(-1.0, 1.0);
  dbl_complex returnme(distribution(generator), distribution(generator));
  return returnme / abs(returnme);
}
}  // namespace bertini