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.

107 lines
2.9 KiB

#pragma once
#include <cmath>
#include "double_extension.hpp"
#include "rational.hpp"
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(i64_rational const& n, unsigned /* precision */) { return rational_cast<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(i64_rational const& n, unsigned /* precision */)
{
return dbl_complex(rational_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);
}