Polyhedral Homotopy Continuation Method for solving sparse polynomial system, optimized by only tracing real zeros
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

141 lines
4.0 KiB

4 days ago
// 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