diff --git a/include/grid_synth/complex.hpp b/include/grid_synth/complex.hpp index e8c267a8..3ff10f2a 100644 --- a/include/grid_synth/complex.hpp +++ b/include/grid_synth/complex.hpp @@ -5,6 +5,8 @@ #include #include +#include "gmp_functions.hpp" + namespace staq { namespace grid_synth { @@ -142,6 +144,11 @@ inline complex operator*(complex x, complex z) { z.a() * x.b() + z.b() * x.a()); } + +inline mpf_class abs(const complex& z) { + return gmpf::sqrt(z.real() * z.real() + z.imag() * z.imag()); +} + } // namespace grid_synth } // namespace staq diff --git a/include/grid_synth/constants.hpp b/include/grid_synth/constants.hpp index 80ec5532..4e53f136 100644 --- a/include/grid_synth/constants.hpp +++ b/include/grid_synth/constants.hpp @@ -2,27 +2,67 @@ #define CONSTANTS_HPP #include -#include #include #include "types.hpp" +#include "gmp_functions.hpp" +#include "complex.hpp" namespace staq { namespace grid_synth { -inline real_t TOL = 1e-14; -inline real_t PI = 3.14159265358979323846264338327950288; // M_PI not standard -inline long int DEFAULT_GMP_PREC = 17; -inline real_t SQRT2 = std::sqrt(2); -inline real_t INV_SQRT2 = 1 / SQRT2; -inline real_t HALF_INV_SQRT2 = 1 / (2 * SQRT2); -inline cplx_t OMEGA(INV_SQRT2, INV_SQRT2); -inline cplx_t OMEGA_CONJ(INV_SQRT2, -INV_SQRT2); -inline cplx_t Im(0, 1); + +struct MultiPrecisionConstants { + real_t tol; + real_t pi; + long int default_gmp_prec; + real_t sqrt2; + real_t inv_sqrt2; + real_t half_inv_sqrt2; + cplx_t omega; + cplx_t omega_conj; + cplx_t im; + real_t log_lambda; + real_t sqrt_lambda; + real_t sqrt_lambda_inv; +}; + +inline MultiPrecisionConstants initialize_constants(long int prec) { + long int default_gmp_prec = 4 * prec + 19; + mpf_set_default_prec(log2(10) * default_gmp_prec); + real_t tol = gmpf::pow(real_t(10), -default_gmp_prec + 2); + real_t pi = gmpf::gmp_pi(); + real_t sqrt2 = gmpf::sqrt(real_t(2)); + real_t inv_sqrt2 = real_t(real_t(1) / sqrt2); + real_t half_inv_sqrt2 = real_t(real_t(1) / (real_t(2) * sqrt2)); + cplx_t omega = cplx_t(inv_sqrt2, inv_sqrt2); + cplx_t omega_conj = cplx_t(inv_sqrt2, -inv_sqrt2); + real_t log_lambda = gmpf::log10(real_t(1)+sqrt2); + real_t sqrt_lambda = gmpf::sqrt(real_t(1)+sqrt2); + real_t sqrt_lambda_inv = sqrt(-real_t(1) + sqrt2); + cplx_t im = cplx_t(real_t(0), real_t(1)); + + return MultiPrecisionConstants{ + tol,pi,default_gmp_prec,sqrt2,inv_sqrt2,half_inv_sqrt2,omega,omega_conj,im, + log_lambda, sqrt_lambda, sqrt_lambda_inv}; +} + +inline MultiPrecisionConstants MP_CONSTS = initialize_constants(10); + +#define TOL MP_CONSTS.tol +#define PI MP_CONSTS.pi +#define DEFAULT_GMP_PREC MP_CONSTS.default_gmp_prec +#define SQRT2 MP_CONSTS.sqrt2 +#define INV_SQRT2 MP_CONSTS.inv_sqrt2 +#define HALF_INV_SQRT2 MP_CONSTS.half_inv_sqrt2 +#define OMEGA MP_CONSTS.omega +#define OMEGA_CONJ MP_CONSTS.omega_conj +#define Im MP_CONSTS.im +#define LOG_LAMBDA MP_CONSTS.log_lambda +#define SQRT_LAMBDA MP_CONSTS.sqrt_lambda +#define SQRT_LAMBDA_INV MP_CONSTS.sqrt_lambda_inv + inline int MAX_ATTEMPTS_POLLARD_RHO = 200; -inline real_t LOG_LAMBDA = std::log10(1 + std::sqrt(2)); -inline real_t SQRT_LAMBDA = std::sqrt(1 + std::sqrt(2)); -inline real_t SQRT_LAMBDA_INV = std::sqrt(-1 + std::sqrt(2)); const int KMIN = 0; const int KMAX = 10000000; @@ -38,127 +78,6 @@ const str_t DEFAULT_TABLE_FILE = "./.s3_table_file.csv"; // on average we only need 2 attempts so 5 is playing it safe const int MAX_ATTEMPTS_SQRT_NEG_ONE = 100; -inline void initialize_constants(int prec, int max_attempts) { - DEFAULT_GMP_PREC = 4 * prec + 19; - mpf_set_default_prec(log2(10) * DEFAULT_GMP_PREC); - // TOL = pow(real_t(10),-DEFAULT_GMP_PREC+2); - TOL = real_t("10e" + std::to_string(-DEFAULT_GMP_PREC + 2)); - SQRT2 = sqrt(real_t(2)); - INV_SQRT2 = real_t(real_t(1) / SQRT2); - HALF_INV_SQRT2 = real_t(real_t(1) / (real_t(2) * SQRT2)); - OMEGA = cplx_t(INV_SQRT2, INV_SQRT2); - OMEGA_CONJ = cplx_t(INV_SQRT2, -INV_SQRT2); - Im = cplx_t(real_t(0), real_t(1)); - MAX_ATTEMPTS_POLLARD_RHO = max_attempts; -} - -/* -// Tolerance for equality when comparing floats. Default is set to guarantee -// known edge cases. -inline real_t TOL = 1e-17; -inline real_t PI = 3.14159265358979323846264338327950288 ; // M_PI not standard -inline long int DEFAULT_GMP_PREC = 500; -inline real_t SQRT2 = sqrt(real_t(2)); -inline real_t INV_SQRT2 = real_t(1) / SQRT2; -inline real_t HALF_INV_SQRT2 = real_t(1) / (real_t(2) * SQRT2); -inline cplx_t OMEGA(INV_SQRT2, INV_SQRT2); -inline cplx_t OMEGA_CONJ(INV_SQRT2, -INV_SQRT2); -inline cplx_t Im(real_t(0), real_t(1)); - -const double LOW_PREC_TOL = 1e-17; // tolerance for float equality as low - // precisions - -const int KMIN = 0; // Default minimum k value -const int KMAX = 10000000; // Default maximum k value -const int COLW = 10; // Default width of columns for data output -const int PREC = 5; // Default precision for output -int MAX_ATTEMPTS_POLLARD_RHO = 200; -const int POLLARD_RHO_INITIAL_ADDEND = 1; -const int POLLARD_RHO_START = 2; -const int MOD_SQRT_MAX_DEPTH = 20; - -const int MAX_ITERATIONS_FERMAT_TEST = 5; -const str_t DEFAULT_TABLE_FILE = "./.s3_table_file.csv"; - -// on average we only need 2 attempts so 5 is playing it safe -const int MAX_ATTEMPTS_SQRT_NEG_ONE = 10; -*/ -/* - * These need to be initialized once the required precision is known. - */ -// class MultiPrecisionConstants { -// private: -// inline static long int gmp_prec_; // precision of gmp ints -// inline static real_t tol_; // tolerance for float equaliy and zero -// checking inline static real_t pi_; inline static real_t sqrt2_; -// inline static real_t inv_sqrt2_; -// inline static real_t half_inv_sqrt2_; -// inline static cplx_t omega_; -// inline static cplx_t omega_conj_; -// inline static cplx_t im_; -// -// inline static mpf_class gmp_pi_(const mpf_class& tol) { -// using namespace std; -// real_t three("3"); -// real_t lasts("0"); -// real_t t = three; -// real_t s("3"); -// real_t n("1"); -// real_t na("0"); -// real_t d("0"); -// real_t da("24"); -// while(abs(s-lasts)>tol) { -// lasts = s; -// n = n+na; -// na = na + mpf_class("8"); -// d = d+da; -// da = da+mpf_class("32"); -// t = (t*n) / d; -// s += t; -// } -// return s; -// } -// -// public: -// /* -// * Accepts the precision required for the solution of the -// approximate -// * synthesis problem -// */ -// static void initialize(long int prec) { -// gmp_prec_ = -4*prec-19; -// mpf_set_default_prec(gmp_prec_); -// tol_ = real_t(str_t(-4*prec-16,'0')+"1"); -// pi_ = gmp_pi_(tol_); -// sqrt2_ = sqrt(real_t(2)); -// inv_sqrt2_ = real_t(1) / sqrt2_; -// half_inv_sqrt2_ = real_t(1) / (real_t(2)*sqrt2_); -// omega_ = cplx_t(inv_sqrt2_,inv_sqrt2_); -// omega_conj_ = cplx_t(inv_sqrt2_, -inv_sqrt2_); -// im_ = cplx_t(real_t(0),real_t(1)); -// } -// -// static long int gmp_prec() noexcept { return gmp_prec_; } -// static real_t tol() noexcept { return tol_; } -// static real_t pi() noexcept { return pi_; } -// static real_t sqrt2() noexcept { return sqrt2_; } -// static real_t inv_sqrt2() noexcept { return inv_sqrt2_; } -// static real_t half_inv_sqrt2() noexcept { return half_inv_sqrt2_; } -// static cplx_t omega() noexcept { return omega_; } -// static cplx_t omega_conj() noexcept { return omega_conj_; } -// static cplx_t im() noexcept { return im_; } -// }; - -// #define PI MultiPrecisionConstants::pi() -// #define TOL MultiPrecisionConstants::tol() -// #define DEFAULT_GMP_PRECISION MultiPrecisionConstants::gmp_prec() -// #define SQRT2 MultiPrecisionConstants::sqrt2() -// #define INV_SQRT2 MultiPrecisionConstants::inv_sqrt2() -// #define HALF_INV_SQRT2 MultiPrecisionConstants::half_inv_sqrt2() -// #define OMEGA MultiPrecisionConstants::omega() -// #define OMEGA_CONJ MultiPrecisionConstants::omega_conj() -// #define Im MultiPrecisionConstants::im() - } // namespace grid_synth } // namespace staq diff --git a/include/grid_synth/diophantine_solver.hpp b/include/grid_synth/diophantine_solver.hpp index a1a45e5d..d670cfe1 100644 --- a/include/grid_synth/diophantine_solver.hpp +++ b/include/grid_synth/diophantine_solver.hpp @@ -67,7 +67,7 @@ inline bool modular_sqrt(int_t& answer, const int_t& a, const int_t& p) { s += 1; } - int_t num_tries = 2 * int_t(log10(p) * log10(p)); + int_t num_tries = 2 * int_t(gmpf::log10(p) * gmpf::log10(p)); int_t j = 0; while ((mod_pow(z, int_t((p - 1) / 2), p) == 1) && (j < num_tries)) { @@ -88,14 +88,14 @@ inline bool modular_sqrt(int_t& answer, const int_t& a, const int_t& p) { int_t depth = 1; while (t != 0 && t != 1) { int_t i = 0; - while ((int_t(mod_pow(t, int_t(pow(2, i)), p))) != 1) { + while ((int_t(mod_pow(t, int_t(gmpf::pow(2, i)), p))) != 1) { i += 1; if (i == m) { return false; } } - int_t b = int_t(mod_pow(c, int_t(pow(2, m - i - 1)), p)); + int_t b = int_t(mod_pow(c, int_t(gmpf::pow(2, m - i - 1)), p)); m = i; c = int_t(mod_pow(b, 2, p)); t = (t * b * b) % p; diff --git a/include/grid_synth/gmp_functions.hpp b/include/grid_synth/gmp_functions.hpp index 3a65b007..06025e7e 100644 --- a/include/grid_synth/gmp_functions.hpp +++ b/include/grid_synth/gmp_functions.hpp @@ -3,27 +3,28 @@ #include #include +#include -#include "complex.hpp" -#include "constants.hpp" #include "utils.hpp" namespace staq { -namespace grid_synth { - -// TODO add high precision sine and cosine functions - -inline mpf_class gmp_pi(const mpf_class& tol) { - using namespace std; - real_t three("3"); - real_t lasts("0"); - real_t t = three; - real_t s("3"); - real_t n("1"); - real_t na("0"); - real_t d("0"); - real_t da("24"); - while (abs(s - lasts) > tol) { +namespace gmpf { + +inline mpf_class gmp_abs(const mpf_class& x) { return sgn(x) * x; } + +// computes +inline mpf_class gmp_pi() { + long int tol_exp = std::log10(2)*mpf_get_default_prec(); + mpf_class tol("1e-" + std::to_string(tol_exp)); + mpf_class three("3"); + mpf_class lasts("0"); + mpf_class t = three; + mpf_class s("3"); + mpf_class n("1"); + mpf_class na("0"); + mpf_class d("0"); + mpf_class da("24"); + while (gmp_abs(s - lasts) > tol) { lasts = s; n = n + na; na = na + mpf_class("8"); @@ -35,37 +36,37 @@ inline mpf_class gmp_pi(const mpf_class& tol) { return s; } -inline mpz_class min(const mpz_class& x, const mpz_class& y) noexcept { +inline mpz_class gmp_min(const mpz_class& x, const mpz_class& y) noexcept { return mpz_class(x * (x < y) + (1 - (x < y)) * y); } -inline mpz_class max(const mpz_class& x, const mpz_class& y) noexcept { +inline mpz_class gmp_max(const mpz_class& x, const mpz_class& y) noexcept { return mpz_class(x * (x > y) + (1 - (x > y)) * y); } -inline mpf_class min(const mpf_class& x, const mpf_class& y) noexcept { +inline mpf_class gmp_min(const mpf_class& x, const mpf_class& y) noexcept { return mpf_class(x * (x < y) + (1 - (x < y)) * y); } -inline mpf_class max(const mpf_class& x, const mpf_class& y) noexcept { +inline mpf_class gmp_max(const mpf_class& x, const mpf_class& y) noexcept { return mpf_class(x * (x > y) + (1 - (x > y)) * y); } -inline mpz_class floor(const mpf_class& x) { +inline mpz_class gmp_floor(const mpf_class& x) { mpf_class f; mpf_floor(f.get_mpf_t(), x.get_mpf_t()); return mpz_class(f); } -inline mpz_class ceil(const mpf_class& x) { +inline mpz_class gmp_ceil(const mpf_class& x) { mpf_class f; mpf_ceil(f.get_mpf_t(), x.get_mpf_t()); return mpz_class(f); } -inline mpz_class round(const mpf_class& x) { +inline mpz_class gmp_round(const mpf_class& x) { mpf_class f, c; mpf_floor(f.get_mpf_t(), x.get_mpf_t()); @@ -96,14 +97,16 @@ inline mpf_class pow(const mpf_class& base, signed long int exponent) { return output; } -inline mpf_class fleq(const mpf_class& lhs, const mpf_class& rhs, - const mpf_class& tol = TOL) { - return (lhs < rhs) || (abs(lhs - rhs) < TOL); +inline mpf_class gmp_leq(const mpf_class& lhs, const mpf_class& rhs) { + long int tol_exp = std::log10(2)*lhs.get_prec(); + mpf_class tol = mpf_class("1e-" + std::to_string(tol_exp)); + return (lhs < rhs) || (abs(lhs - rhs) < tol); } -inline mpf_class fgeq(const mpf_class& lhs, const mpf_class& rhs, - const mpf_class& tol = TOL) { - return (lhs > rhs) || (abs(lhs - rhs) < TOL); +inline mpf_class gmp_geq(const mpf_class& lhs, const mpf_class& rhs) { + long int tol_exp = std::log10(2)*lhs.get_prec(); + mpf_class tol = mpf_class("1e-" + std::to_string(tol_exp)); + return (lhs > rhs) || (abs(lhs - rhs) < tol); } /* @@ -166,56 +169,59 @@ inline mpf_class log2(const mpf_class& x) { */ inline mpf_class reduce_angle(const mpf_class& phi) { mpf_class result = phi; - while (result > PI) - result -= 2 * PI; - while (result < -PI) - result += 2 * PI; + mpf_class pi = gmp_pi(); + while (result > pi) + result -= mpf_class("2") * pi; + while (result < pi) + result += mpf_class("2") * pi; return result; } -inline mpf_class sin(const mpf_class& theta, const mpf_class& tol = TOL) { - real_t phi = reduce_angle(theta); - mpz_class i = 1; - mpf_class lasts = 0; +//TODO improve accuracy of sin and cos by expanding about pi/2, pi/3, pi/4, and pi/6 +//depending on the input angle. +inline mpf_class sin(const mpf_class& theta) { + long int initial_prec = theta.get_prec(); + long int tol_exp = std::log10(2)*initial_prec; + mpf_class phi = reduce_angle(theta); + mpz_class i(1); + mpf_class lasts(0); mpf_class s = phi; - mpf_class fact = 1; - mpf_class num = phi; - mpf_class sign = 1; - while (abs(s - lasts) > tol) { + mpf_class fact(1); + mpf_class num(phi); + mpf_class sign(1); + while (gmp_abs(s - lasts) > mpf_class("1e-" + std::to_string(tol_exp))) { lasts = s; - i += 2; - fact *= i * (i - 1); + i += mpf_class("2"); + fact *= i * (i - mpf_class("1")); num *= phi * phi; - sign *= -1; + sign *= mpf_class("-1"); s += sign * (num / fact); } return s; } -inline mpf_class cos(const mpf_class& theta, const mpf_class& tol = TOL) { - real_t phi = reduce_angle(theta); - mpz_class i = 0; - mpf_class lasts = 0; - mpf_class s = 1; - mpf_class fact = 1; - mpf_class num = 1; - mpf_class sign = 1; - while (abs(s - lasts) > tol) { + +inline mpf_class cos(const mpf_class& theta) { + long int initial_prec = theta.get_prec(); + long int tol_exp = std::log10(2)*theta.get_prec(); + mpf_class phi = reduce_angle(theta); + mpz_class i(0); + mpf_class lasts(0); + mpf_class s(1); + mpf_class fact(1); + mpf_class num(1); + mpf_class sign(1); + while (gmp_abs(s - lasts) > mpf_class("1e-" + std::to_string(tol_exp))) { lasts = s; - i += 2; - fact *= i * (i - 1); + i += mpf_class("2"); + fact *= i * (i - mpf_class("1")); num *= phi * phi; - sign *= -1; + sign *= mpf_class("-1"); s += sign * (num / fact); } return s; } -inline mpf_class abs(const complex& z) { - return sqrt(z.real() * z.real() + z.imag() * z.imag()); -} - -inline mpf_class abs(const mpf_class& x) { return sgn(x) * x; } inline mpf_class sqrt(const mpf_class& x) { mpf_class output; @@ -223,7 +229,7 @@ inline mpf_class sqrt(const mpf_class& x) { return output; } -} // namespace grid_synth +} // namespace gmpf } // namespace staq #endif // GMP_FUNCTIONS_HPP diff --git a/include/grid_synth/grid_solvers.hpp b/include/grid_synth/grid_solvers.hpp index 5038dfcb..b70b659e 100644 --- a/include/grid_synth/grid_solvers.hpp +++ b/include/grid_synth/grid_solvers.hpp @@ -27,17 +27,16 @@ namespace grid_synth { template inline int_t lower_bound_a(const real_t& xlo, const int_t& b, const real_t& tol) { - using namespace std; real_t lowera_double = xlo - (b * SQRT2); real_t decimal; int_t intpart; - decimal = abs(decimal_part(lowera_double, intpart)); + decimal = gmpf::gmp_abs(gmpf::decimal_part(lowera_double, intpart)); if ((lowera_double < 0) && ((1 - decimal) < tol)) - return floor(lowera_double); + return gmpf::gmp_floor(lowera_double); if ((lowera_double > 0) && (decimal < tol)) - return floor(lowera_double); + return gmpf::gmp_floor(lowera_double); - return ceil(lowera_double); + return gmpf::gmp_ceil(lowera_double); } /** @@ -47,15 +46,14 @@ inline int_t lower_bound_a(const real_t& xlo, const int_t& b, template inline int_t upper_bound_a(const bound_t& xhi, const int_t& b, const real_t& tol) { - using namespace std; real_t uppera_double = xhi - (b * SQRT2); real_t decimal; int_t intpart; - decimal = abs(decimal_part(uppera_double, intpart)); + decimal = gmpf::gmp_abs(gmpf::decimal_part(uppera_double, intpart)); if ((uppera_double > 0) && ((1 - decimal) < tol)) - return ceil(uppera_double); + return gmpf::gmp_ceil(uppera_double); - return floor(uppera_double); + return gmpf::gmp_floor(uppera_double); } /** @@ -65,17 +63,16 @@ inline int_t upper_bound_a(const bound_t& xhi, const int_t& b, template inline int_t lower_bound_b(const bound_t& xlo, const bound_t& yhi, const bound_t& tol) { - using namespace std; real_t lowerb_double = (xlo - yhi) * HALF_INV_SQRT2; real_t decimal; int_t intpart; - decimal = abs(decimal_part(lowerb_double, intpart)); + decimal = gmpf::gmp_abs(gmpf::decimal_part(lowerb_double, intpart)); if ((lowerb_double) < 0 && ((1 - decimal) < tol)) - return floor(lowerb_double); + return gmpf::gmp_floor(lowerb_double); if ((lowerb_double) > 0 && (decimal < tol)) - return floor(lowerb_double); + return gmpf::gmp_floor(lowerb_double); - return ceil(lowerb_double); + return gmpf::gmp_ceil(lowerb_double); } /** @@ -85,15 +82,14 @@ inline int_t lower_bound_b(const bound_t& xlo, const bound_t& yhi, template inline int_t upper_bound_b(const bound_t& xhi, const bound_t& ylo, const real_t& tol) { - using namespace std; real_t upperb_double = (xhi - ylo) * HALF_INV_SQRT2; real_t decimal; int_t intpart; - decimal = abs(decimal_part(upperb_double, intpart)); + decimal = gmpf::gmp_abs(gmpf::decimal_part(upperb_double, intpart)); if ((upperb_double > 0) && (1 - decimal < tol)) - return ceil(upperb_double); + return gmpf::gmp_ceil(upperb_double); - return floor(upperb_double); + return gmpf::gmp_floor(upperb_double); } /** @@ -106,7 +102,7 @@ inline int_t upper_bound_b(const bound_t& xhi, const bound_t& ylo, */ template inline int_t find_scale_exponent(const Interval& interval) { - real_t ratio = log10(interval.width()) / LOG_LAMBDA; + real_t ratio = gmpf::log10(interval.width()) / LOG_LAMBDA; return int_t(ratio) + 1; } @@ -159,8 +155,6 @@ template inline zsqrt2_vec_t oneD_scaled_grid_solver(const Interval& A, const Interval& B, const real_t tol = TOL) { - using namespace std; - zsqrt2_vec_t solns; int_t k = find_scale_exponent(A); Interval scaled_A = A; diff --git a/include/grid_synth/regions.hpp b/include/grid_synth/regions.hpp index 560fcf51..80b10ca3 100644 --- a/include/grid_synth/regions.hpp +++ b/include/grid_synth/regions.hpp @@ -75,7 +75,7 @@ class Interval { bool contains(const bound_t& x, const real_t tol = TOL) const { return ((hi_ - x) * (x - lo_) > 0) || - (abs((hi_ - x) * (x - lo_)) < tol); + (gmpf::gmp_abs((hi_ - x) * (x - lo_)) < tol); } Interval operator+(const bound_t& shift_factor) const { @@ -225,16 +225,16 @@ class Ellipse { real_t e_; void get_z_and_e_() { - z_ = real_t(real_t(real_t("0.5") * log10(real_t(D_(1, 1) / D_(0, 0)))) / + z_ = real_t(real_t(real_t("0.5") * gmpf::log10(real_t(D_(1, 1) / D_(0, 0)))) / LOG_LAMBDA); - e_ = sqrt(D_(1, 1) * D_(0, 0)); + e_ = gmpf::sqrt(D_(1, 1) * D_(0, 0)); } mat_t get_mat_from_axes_(const real_t& semi_major_axis, const real_t& semi_minor_axis, const real_t& angle) { - real_t ct = cos(angle); - real_t st = sin(angle); + real_t ct = gmpf::cos(angle); + real_t st = gmpf::sin(angle); real_t inva = real_t("1") / semi_minor_axis; real_t invb = real_t("1") / semi_major_axis; @@ -268,7 +268,7 @@ class Ellipse { // <<(T * msq - sqrt(T * T * msq * msq - 4 * msq))<< endl; real_t a1 = 0; real_t a2 = 0; - if (abs(T * T * msq * msq - real_t("4") * msq) < TOL) { + if (gmpf::gmp_abs(T * T * msq * msq - real_t("4") * msq) < TOL) { a1 = sqrt(T * msq) / real_t("2"); a2 = sqrt(T * msq) / real_t("2"); } else { @@ -331,7 +331,7 @@ class Ellipse { */ Ellipse(const real_t& angle, const real_t& eps) : angle_(angle) { real_t r0 = (real_t("3") - eps * eps) / real_t("3"); - center_ = vec_t{r0 * cos(angle), r0 * sin(angle)}; + center_ = vec_t{r0 * gmpf::cos(angle), r0 * gmpf::sin(angle)}; semi_major_axis_ = (real_t("2") / sqrt(real_t("3")) * eps * sqrt(real_t("1") - (eps * eps / real_t("4")))); @@ -365,8 +365,8 @@ class Ellipse { void rescale(const real_t scale) { D_ = (real_t("1") / (scale * scale)) * D_; - semi_minor_axis_ = semi_minor_axis_ * abs(scale); - semi_major_axis_ = semi_major_axis_ * abs(scale); + semi_minor_axis_ = semi_minor_axis_ * gmpf::gmp_abs(scale); + semi_major_axis_ = semi_major_axis_ * gmpf::gmp_abs(scale); center_ = scale * center_; get_z_and_e_(); @@ -385,7 +385,7 @@ class Ellipse { bool contains(const vec_t& point, const real_t& tol = TOL) const { using namespace std; real_t x = (point - center_).transpose() * D_ * (point - center_); - return (x < real_t("1")) || (abs(x - real_t("1")) < tol); + return (x < real_t("1")) || (gmpf::gmp_abs(x - real_t("1")) < tol); } bool contains(const real_t& x, const real_t& y, diff --git a/include/grid_synth/rings.hpp b/include/grid_synth/rings.hpp index fc8da73b..e207992b 100644 --- a/include/grid_synth/rings.hpp +++ b/include/grid_synth/rings.hpp @@ -52,14 +52,14 @@ class ZSqrt2 { } ZSqrt2 self_sqrt() const { - using namespace std; + using namespace gmpf; real_t a, b; - int_t b_squared_plus = (a_ + sqrt((*this).norm())) / 4; - int_t b_squared_minus = (a_ - sqrt((*this).norm())) / 4; + int_t b_squared_plus = int_t((a_ + sqrt((*this).norm())) / real_t(4)); + int_t b_squared_minus = int_t((a_ - sqrt((*this).norm())) / real_t(4)); - if (pow(round(sqrt(b_squared_plus)), 2) == b_squared_plus) + if (pow(gmp_round(sqrt(b_squared_plus)), 2) == b_squared_plus) b = sqrt(b_squared_plus); - else if (pow(round(sqrt(b_squared_minus)), 2) == b_squared_minus) + else if (pow(gmp_round(sqrt(b_squared_minus)), 2) == b_squared_minus) b = sqrt(b_squared_minus); else { std::cout << "(" << a_ << "," << b_ << ")" @@ -76,7 +76,7 @@ class ZSqrt2 { } else a = b_ / (2 * b); - return ZSqrt2(round(a), round(b)); + return ZSqrt2(gmpf::gmp_round(a), gmpf::gmp_round(b)); } // Arithmetic operators @@ -99,8 +99,8 @@ class ZSqrt2 { ZSqrt2 operator/(const ZSqrt2& Z) const { using namespace std; real_t mag = Z.norm(); - int_t a = round(real_t(a_ * Z.a() - 2 * b_ * Z.b()) / mag); - int_t b = round(real_t(b_ * Z.a() - a_ * Z.b()) / mag); + int_t a = gmpf::gmp_round(real_t(a_ * Z.a() - 2 * b_ * Z.b()) / mag); + int_t b = gmpf::gmp_round(real_t(b_ * Z.a() - a_ * Z.b()) / mag); return ZSqrt2(a, b); } @@ -305,7 +305,7 @@ class ZOmega { (c_ - a_) / 2); } - bool w() const { return round(w_.real()) != 0; } + bool w() const { return gmpf::gmp_round(w_.real()) != 0; } ZOmega dot() const { return ZOmega(-a_, b_, -c_, d_); } ZOmega conj() const { return ZOmega(-c_, -b_, -a_, d_); } @@ -372,7 +372,7 @@ class ZOmega { c_ += Z.c(); d_ += Z.d(); - ZSqrt2 shift(0, round(w_.real()) != 0 && Z.w()); + ZSqrt2 shift(0, gmpf::gmp_round(w_.real()) != 0 && Z.w()); alpha_ += Z.alpha() + shift; beta_ += Z.beta() + shift; @@ -387,7 +387,7 @@ class ZOmega { c_ -= Z.c(); d_ -= Z.d(); - ZSqrt2 shift(0, round(w_.real()) != 0 && Z.w()); + ZSqrt2 shift(0, gmpf::gmp_round(w_.real()) != 0 && Z.w()); alpha_ -= Z.alpha() + shift; beta_ -= Z.beta() + shift; @@ -406,7 +406,7 @@ class ZOmega { // c_ = -a_*Z.b() - b_*Z.a() + c_*Z.d() + d_*Z.c(); // d_ = -a_*Z.c() - b_*Z.b() - c_*Z.a() + d_*Z.d(); // - // ZSqrt2 shift(0, round(w_.real()) != 0 && Z.w()); + // ZSqrt2 shift(0, gmpf::gmp_round(w_.real()) != 0 && Z.w()); // alpha_ *= Z.alpha() + shift; // beta_ *= Z.beta() + shift; diff --git a/include/grid_synth/rz_approximation.hpp b/include/grid_synth/rz_approximation.hpp index 708dea7a..3ff66e07 100644 --- a/include/grid_synth/rz_approximation.hpp +++ b/include/grid_synth/rz_approximation.hpp @@ -37,11 +37,11 @@ class RzApproximation { const int_t scale_exponent, const real_t theta, const real_t eps) : matrix_(u, t, scale_exponent, 0), eps_(eps), solution_found_(true) { - u_val_ = cplx_t(u.decimal().real() / pow(SQRT2, scale_exponent), - u.decimal().imag() / pow(SQRT2, scale_exponent)); - t_val_ = cplx_t(t.decimal().real() / pow(SQRT2, scale_exponent), - t.decimal().imag() / pow(SQRT2, scale_exponent)); - z_ = cplx_t(cos(theta), sin(theta)); + u_val_ = cplx_t(u.decimal().real() / gmpf::pow(SQRT2, scale_exponent), + u.decimal().imag() / gmpf::pow(SQRT2, scale_exponent)); + t_val_ = cplx_t(t.decimal().real() / gmpf::pow(SQRT2, scale_exponent), + t.decimal().imag() / gmpf::pow(SQRT2, scale_exponent)); + z_ = cplx_t(gmpf::cos(theta), gmpf::sin(theta)); } DOmegaMatrix matrix() const { return matrix_; } @@ -87,11 +87,11 @@ inline RzApproximation find_rz_approximation(const real_t& theta, real_t scaleB; while ((!solution_found) && k < max_k) { if (k % 2 == 0) { - scaleA = pow(2, k / 2); - scaleB = pow(2, k / 2); + scaleA = gmpf::pow(2, k / 2); + scaleB = gmpf::pow(2, k / 2); } else { - scaleA = pow(2, (k - 1) / 2) * SQRT2; - scaleB = -pow(2, (k - 1) / 2) * SQRT2; + scaleA = gmpf::pow(2, (k - 1) / 2) * SQRT2; + scaleB = -gmpf::pow(2, (k - 1) / 2) * SQRT2; } state[0].rescale(scaleA); state[1].rescale(scaleB); @@ -108,7 +108,7 @@ inline RzApproximation find_rz_approximation(const real_t& theta, candidate = candidate.reduce(); } - ZSqrt2 xi = ZSqrt2(int_t(pow(2, temp_k)), 0) - + ZSqrt2 xi = ZSqrt2(int_t(gmpf::pow(2, temp_k)), 0) - (candidate.conj() * candidate).to_zsqrt2(); ZOmega answer(0); solution_found = diophantine_solver(answer, xi); @@ -130,13 +130,12 @@ inline RzApproximation find_fast_rz_approximation(const real_t& theta, const int_t& kmin = KMIN, const int_t& kmax = KMAX, const real_t tol = TOL) { - using namespace std; // int_t k = 3 * int_t(log(real_t(1) / eps) / log(2)) / 2; int_t k = kmin; int_t max_k = kmax; bool solution_found = false; zomega_vec_t scaled_candidates; - vec_t z{cos(theta), sin(theta)}; + vec_t z{gmpf::cos(theta), gmpf::sin(theta)}; Ellipse eps_region(theta, eps); Ellipse disk(real_t("0"), real_t("0"), real_t("1"), real_t("1"), real_t("0")); @@ -149,11 +148,11 @@ inline RzApproximation find_fast_rz_approximation(const real_t& theta, while ((!solution_found) && (k < max_k)) { if (k % 2 == 0) { - scaleA = real_t(pow(2, k / 2)); - scaleB = real_t(pow(2, k / 2)); + scaleA = gmpf::pow(2, k / 2); + scaleB = gmpf::pow(2, k / 2); } else { - scaleA = real_t(pow(2, (k - 1) / 2)) * SQRT2; - scaleB = -real_t(pow(2, (k - 1) / 2)) * SQRT2; + scaleA = gmpf::pow(2, (k - 1) / 2) * SQRT2; + scaleB = -gmpf::pow(2, (k - 1) / 2) * SQRT2; } bboxA.rescale(scaleA, scaleA); @@ -182,7 +181,7 @@ inline RzApproximation find_fast_rz_approximation(const real_t& theta, candidate = candidate.reduce(); } - ZSqrt2 xi = ZSqrt2(int_t(pow(2, temp_k)), 0) - + ZSqrt2 xi = ZSqrt2(int_t(gmpf::pow(2, temp_k)), 0) - (candidate.conj() * candidate).to_zsqrt2(); ZOmega answer(0); solution_found = diophantine_solver(answer, xi); @@ -211,7 +210,7 @@ inline RzApproximation find_fast_rz_approximation(const real_t& theta, candidate = candidate.reduce(); } - ZSqrt2 xi = ZSqrt2(int_t(pow(2, temp_k)), 0) - + ZSqrt2 xi = ZSqrt2(int_t(gmpf::pow(2, temp_k)), 0) - (candidate.conj() * candidate).to_zsqrt2(); ZOmega answer(0); solution_found = diophantine_solver(answer, xi); @@ -230,122 +229,6 @@ inline RzApproximation find_fast_rz_approximation(const real_t& theta, return RzApproximation(); } -inline RzApproximation verbose_find_fast_rz_approximation( - const real_t& theta, const real_t& eps, const int_t& kmin = KMIN, - const int_t& kmax = KMAX, const real_t tol = TOL) { - using namespace std; - // int_t k = 3 * int_t(log(real_t(1) / eps) / log(2)) / 2; - // int_t k = kmin; - int_t k = 0; - int_t max_k = kmax; - bool solution_found = false; - zomega_vec_t scaled_candidates; - vec_t z{cos(theta), sin(theta)}; - Ellipse eps_region(theta, eps); - Ellipse disk(real_t("0"), real_t("0"), real_t("1"), real_t("1"), - real_t("0")); - state_t state{eps_region, disk}; - SpecialGridOperator G = optimize_skew(state); - real_t scaleA, scaleB; - - UprightRectangle bboxA = state[0].bounding_box(); - UprightRectangle bboxB = state[1].bounding_box(); - cout << "theta = " << theta << endl; - while ((!solution_found) && (k < max_k)) { - if (k % 2 == 0) { - scaleA = real_t(pow(2, k / 2)); - scaleB = real_t(pow(2, k / 2)); - } else { - scaleA = real_t(pow(2, (k - 1) / 2)) * SQRT2; - scaleB = -real_t(pow(2, (k - 1) / 2)) * SQRT2; - } - - cout << "k = " << k << endl; - bboxA.rescale(scaleA, scaleA); - bboxB.rescale(scaleB, scaleB); - - // Interval A_x = bboxA.x_interval().fatten(eps); - // Interval B_x = bboxB.x_interval().fatten(eps); - - // Interval A_y = bboxA.y_interval().fatten(eps); - // Interval B_y = bboxB.y_interval().fatten(eps); - Interval A_x = bboxA.x_interval(); - Interval B_x = bboxB.x_interval(); - - Interval A_y = bboxA.y_interval(); - Interval B_y = bboxB.y_interval(); - - // real_t min_width = (real_t(1) + SQRT2) * (real_t(1) + SQRT2); - - // if((A_x.width()*B_x.width() < min_width) || - // (A_y.width()*B_y.width() < min_width)) { - // k++; - // bboxA.rescale(1 / scaleA, 1 / scaleA); - // bboxB.rescale(1 / scaleB, 1 / scaleB); - // continue; - // } - - zsqrt2_vec_t alpha_solns = oneD_optimal_grid_solver(A_x, B_x, tol); - zsqrt2_vec_t beta_solns = oneD_optimal_grid_solver(A_y, B_y, tol); - - for (auto alpha_soln : alpha_solns) - for (auto beta_soln : beta_solns) { - ZOmega candidate = G * ZOmega(alpha_soln, beta_soln, 0); - if (((candidate.real() / scaleA) * z[0] + - (candidate.imag() / scaleA) * z[1]) > - real_t("1") - (eps * eps / real_t("2"))) { - int_t temp_k = k; - while (candidate.is_reducible()) { - temp_k -= 1; - candidate = candidate.reduce(); - } - - ZSqrt2 xi = ZSqrt2(int_t(pow(2, temp_k)), 0) - - (candidate.conj() * candidate).to_zsqrt2(); - ZOmega answer(0); - solution_found = diophantine_solver(answer, xi); - if (solution_found) { - return RzApproximation(candidate, answer, temp_k, theta, - eps); - } - } - } - - zsqrt2_vec_t shifted_alpha_solns = - oneD_optimal_grid_solver(A_x - INV_SQRT2, B_x + INV_SQRT2, tol); - zsqrt2_vec_t shifted_beta_solns = - oneD_optimal_grid_solver(A_y - INV_SQRT2, B_y + INV_SQRT2, tol); - - for (auto alpha_soln : shifted_alpha_solns) - for (auto beta_soln : shifted_beta_solns) { - ZOmega candidate = G * ZOmega(alpha_soln, beta_soln, 1); - if (((candidate.real() / scaleA) * z[0] + - (candidate.imag() / scaleA) * z[1]) > - real_t("1") - (eps * eps / real_t("2"))) { - int_t temp_k = k; - while (candidate.is_reducible()) { - temp_k -= 1; - candidate = candidate.reduce(); - } - - ZSqrt2 xi = ZSqrt2(int_t(pow(2, temp_k)), 0) - - (candidate.conj() * candidate).to_zsqrt2(); - ZOmega answer(0); - solution_found = diophantine_solver(answer, xi); - if (solution_found) { - return RzApproximation(candidate, answer, temp_k, theta, - eps); - } - } - } - - k++; - bboxA.rescale(real_t("1") / scaleA, real_t("1") / scaleA); - bboxB.rescale(real_t("1") / scaleB, real_t("1") / scaleB); - } - return RzApproximation(); -} - } // namespace grid_synth } // namespace staq diff --git a/include/grid_synth/states.hpp b/include/grid_synth/states.hpp index 5757d197..744dd1e8 100644 --- a/include/grid_synth/states.hpp +++ b/include/grid_synth/states.hpp @@ -32,18 +32,18 @@ inline int_t determine_shift(const state_t state) { inline mat_t sigma(int_t k) { if (k < 0) - return pow(SQRT_LAMBDA_INV, -k) * - mat_t{1, 0, 0, pow(LAMBDA.decimal(), -k)}; + return gmpf::pow(SQRT_LAMBDA_INV, -k) * + mat_t{1, 0, 0, gmpf::pow(LAMBDA.decimal(), -k)}; - return pow(SQRT_LAMBDA_INV, k) * mat_t{pow(LAMBDA.decimal(), k), 0, 0, 1}; + return gmpf::pow(SQRT_LAMBDA_INV, k) * mat_t{gmpf::pow(LAMBDA.decimal(), k), 0, 0, 1}; } inline mat_t tau(int_t k) { if (k < 0) - return pow(SQRT_LAMBDA_INV, -k) * - mat_t{pow(LAMBDA.decimal(), -k), 0, 0, pow(-1, -k)}; + return gmpf::pow(SQRT_LAMBDA_INV, -k) * + mat_t{gmpf::pow(LAMBDA.decimal(), -k), 0, 0, gmpf::pow(-1, -k)}; - return pow(SQRT_LAMBDA_INV, k) * mat_t{1, 0, 0, pow(-LAMBDA.decimal(), k)}; + return gmpf::pow(SQRT_LAMBDA_INV, k) * mat_t{1, 0, 0, gmpf::pow(-LAMBDA.decimal(), k)}; } /* @@ -58,7 +58,6 @@ inline state_t shift(const state_t& state, const int_t& k) { * Reduces skew(state) by 10% and returns the operator that did it. */ inline SpecialGridOperator reduce_skew(state_t& state) { - using namespace std; real_t initial_skew = skew(state); if (initial_skew < 15) return ID; @@ -85,24 +84,24 @@ inline SpecialGridOperator reduce_skew(state_t& state) { real_t z = state[0].z(); real_t zeta = state[1].z(); - if (static_cast(fgeq(state[0].D(0, 1), real_t("0")))) { - if (fgeq(z, real_t("-0.8")) && fleq(z, real_t("0.8")) && - fgeq(zeta, real_t("-0.8")) && fleq(zeta, real_t("0.8"))) { + if (static_cast(gmpf::gmp_geq(state[0].D(0, 1), real_t("0")))) { + if (gmpf::gmp_geq(z, real_t("-0.8")) && gmpf::gmp_leq(z, real_t("0.8")) && + gmpf::gmp_geq(zeta, real_t("-0.8")) && gmpf::gmp_leq(zeta, real_t("0.8"))) { G = G * R; state = R * state; // std::cout << "R" << std::endl; - } else if (fleq(z, real_t("0.3")) && fgeq(zeta, real_t("0.8"))) { + } else if (gmpf::gmp_leq(z, real_t("0.3")) && gmpf::gmp_geq(zeta, real_t("0.8"))) { G = G * K; state = K * state; // std::cout << "K" << std::endl; - } else if (fgeq(z, real_t("0.3")) && fgeq(zeta, real_t("0.3"))) { - real_t c = min(z, zeta); + } else if (gmpf::gmp_geq(z, real_t("0.3")) && gmpf::gmp_geq(zeta, real_t("0.3"))) { + real_t c = gmpf::gmp_min(z, zeta); int_t n = - max(1, int_t(floor(pow(LAMBDA.decimal(), c.get_ui()) / 2))); + gmpf::gmp_max(1, int_t(gmpf::gmp_floor(gmpf::pow(LAMBDA.decimal(), c.get_ui()) / 2))); G = G * A(n); state = A(n) * state; // std::cout << "A^" << n << std::endl; - } else if (fgeq(z, real_t("0.8")) && fleq(zeta, real_t("0.3"))) { + } else if (gmpf::gmp_geq(z, real_t("0.8")) && gmpf::gmp_leq(zeta, real_t("0.3"))) { G = G * K.dot(); state = K.dot() * state; // std::cout << "K.dot()" << std::endl; @@ -116,15 +115,15 @@ inline SpecialGridOperator reduce_skew(state_t& state) { exit(EXIT_FAILURE); } } else { - if (fgeq(z, real_t("-0.8")) && fleq(z, real_t("0.8")) && - fgeq(zeta, real_t("-0.8")) && fleq(zeta, real_t("0.8"))) { + if (gmpf::gmp_geq(z, real_t("-0.8")) && gmpf::gmp_leq(z, real_t("0.8")) && + gmpf::gmp_geq(zeta, real_t("-0.8")) && gmpf::gmp_leq(zeta, real_t("0.8"))) { G = G * R; state = R * state; // std::cout << "R" << std::endl; - } else if (fgeq(z, real_t("-0.2")) && fgeq(zeta, real_t("-0.2"))) { - real_t c = min(z, zeta); + } else if (gmpf::gmp_geq(z, real_t("-0.2")) && gmpf::gmp_geq(zeta, real_t("-0.2"))) { + real_t c = gmpf::gmp_min(z, zeta); int_t n = - max(1, int_t(floor(pow(LAMBDA.decimal(), c.get_ui()) / 2))); + gmpf::gmp_max(1, int_t(gmpf::gmp_floor(gmpf::pow(LAMBDA.decimal(), c.get_ui()) / 2))); G = G * B(n); state = B(n) * state; // std::cout << "B^" << n << std::endl; @@ -154,7 +153,6 @@ inline SpecialGridOperator reduce_skew(state_t& state) { * original normalization but with the skew reduced to its lowest possible value */ inline SpecialGridOperator optimize_skew(state_t& state) { - using namespace std; real_t scaleA = state[0].normalize(); real_t scaleB = state[1].normalize(); SpecialGridOperator G = ID; diff --git a/junk b/junk new file mode 100644 index 00000000..f26bbe46 --- /dev/null +++ b/junk @@ -0,0 +1,8 @@ +0.70710678118654752 440084436210484903928483593768847403658833986899536623923105 96259259108747347952536 +0.70710678118654752 440084436210484903928483593768847403658833986899536623923105 35194251937671638207864 + +0.70710678118654752 440084436210484903928483593768847403658833986899536623923105 35194251937671638207864 +0.70710678118654752 +440084436210484903928483593768847403658833986899536623923105i 35194251937671638207864 +0.70710678118654752 440084436210484903928483593768847403658833986899536623923105 35194251937671638207864 +0.70710678118654752 64172236478881629202753737022098913577754042624199671058910420861761390431913749669 diff --git a/tools/grid_synth.cpp b/tools/grid_synth.cpp index 2acd8d72..321abdb4 100644 --- a/tools/grid_synth.cpp +++ b/tools/grid_synth.cpp @@ -33,16 +33,17 @@ #include "grid_synth/exact_synthesis.hpp" #include "grid_synth/rz_approximation.hpp" #include "grid_synth/types.hpp" +#include "grid_synth/regions.hpp" int main(int argc, char** argv) { + using namespace staq; using namespace grid_synth; - using namespace std; - bool check, details, verbose, timer, more_verbose; + bool check=false, details=false, verbose=false, timer=false; real_t theta, eps; - vector thetas; - vector prec_lst; + std::vector thetas; + std::vector prec_lst; long int prec; int factor_effort; domega_matrix_table_t s3_table; @@ -79,41 +80,36 @@ int main(int argc, char** argv) { app.add_flag("-v, --verbose", verbose, "Include additional output during runtime such as runtime " "parameters and update on each step."); - app.add_flag("-V, --more-verbose", more_verbose, - "Output update on approximation function"); app.add_flag("--time", timer, "Time program"); CLI11_PARSE(app, argc, argv); if (verbose) { - cout << thetas.size() << " angles read." << '\n'; + std::cout << thetas.size() << " angles read." << '\n'; } - if (more_verbose) { - verbose = true; - } if (*read) { if (verbose) { - cout << "Reading s3_table from " << tablefile << '\n'; + std::cout << "Reading s3_table from " << tablefile << '\n'; } s3_table = read_s3_table(tablefile); } else if (*write) { if (verbose) { - cout << "Generating new table file and writing to " << tablefile + std::cout << "Generating new table file and writing to " << tablefile << '\n'; } s3_table = generate_s3_table(); write_s3_table(tablefile, s3_table); - } else if (ifstream(DEFAULT_TABLE_FILE)) { + } else if (std::ifstream(DEFAULT_TABLE_FILE)) { if (verbose) { - cout << "Table file found at default location " + std::cout << "Table file found at default location " << DEFAULT_TABLE_FILE << '\n'; } s3_table = read_s3_table(DEFAULT_TABLE_FILE); } else { if (verbose) { - cout << "Failed to find " << DEFAULT_TABLE_FILE + std::cout << "Failed to find " << DEFAULT_TABLE_FILE << ". Generating new table file and writing to " << DEFAULT_TABLE_FILE << '\n'; } @@ -121,147 +117,132 @@ int main(int argc, char** argv) { write_s3_table(DEFAULT_TABLE_FILE, s3_table); } - DEFAULT_GMP_PREC = 4 * prec + 19; - mpf_set_default_prec(log2(10) * DEFAULT_GMP_PREC); - TOL = pow(real_t(10), -DEFAULT_GMP_PREC + 2); - PI = gmp_pi(TOL); - SQRT2 = sqrt(real_t(2)); - INV_SQRT2 = real_t(real_t(1) / SQRT2); - HALF_INV_SQRT2 = real_t(real_t(1) / (real_t(2) * SQRT2)); - OMEGA = cplx_t(INV_SQRT2, INV_SQRT2); - OMEGA_CONJ = cplx_t(INV_SQRT2, -INV_SQRT2); - LOG_LAMBDA = log10(LAMBDA.decimal()); - SQRT_LAMBDA = sqrt(LAMBDA.decimal()); - SQRT_LAMBDA_INV = sqrt(LAMBDA_INV.decimal()); - Im = cplx_t(real_t(0), real_t(1)); - eps = pow(real_t(10), -prec); + MP_CONSTS = initialize_constants(prec); + eps = gmpf::pow(real_t(10), -prec); + + std::cout << std::fixed << std::setprecision(100) << SQRT2 << std::endl; if (*fact_eff) { MAX_ATTEMPTS_POLLARD_RHO = factor_effort; } if (verbose) { - cout << "Runtime Parameters" << '\n'; - cout << "------------------" << '\n'; - cout << setw(3 * COLW) << left << "TOL (Tolerance for float equality) " - << setw(1) << ": " << setw(3 * COLW) << left << scientific << TOL + std::cout << "Runtime Parameters" << '\n'; + std::cout << "------------------" << '\n'; + std::cout << std::setw(3 * COLW) << std::left << "TOL (Tolerance for float equality) " + << std::setw(1) << ": " << std::setw(3 * COLW) << std::left << std::scientific << TOL << '\n'; - cout << setw(3 * COLW) << left << "KMIN (Minimum scaling exponent) " - << setw(1) << ": " << setw(3 * COLW) << left << fixed << KMIN + std::cout << std::setw(3 * COLW) << std::left << "KMIN (Minimum scaling exponent) " + << std::setw(1) << ": " << std::setw(3 * COLW) << std::left << std::fixed << KMIN << '\n'; - cout << setw(2 * COLW) << left << "KMAX (Maximum scaling exponent) " - << setw(1) << ": " << setw(3 * COLW) << left << fixed << KMAX + std::cout << std::setw(2 * COLW) << std::left << "KMAX (Maximum scaling exponent) " + << std::setw(1) << ": " << std::setw(3 * COLW) << std::left << std::fixed << KMAX << '\n'; - cout << setw(3 * COLW) << left + std::cout << std::setw(3 * COLW) << std::left << "MAX_ATTEMPTS_POLLARD_RHO (How hard we try to factor) " - << setw(1) << ": " << setw(3 * COLW) << left + << std::setw(1) << ": " << std::setw(3 * COLW) << std::left << MAX_ATTEMPTS_POLLARD_RHO << '\n'; - cout << setw(3 * COLW) << left + std::cout << std::setw(3 * COLW) << std::left << "MAX_ITERATIONS_FERMAT_TEST (How hard we try to check " "primality) " - << setw(1) << ": " << setw(3 * COLW) << left + << std::setw(1) << ": " << std::setw(3 * COLW) << std::left << MAX_ITERATIONS_FERMAT_TEST << '\n'; } - cout << scientific; + std::cout << std::scientific; long long duration = 0; if (*prec_opt && *thetas_op) { - random_device rd; + std::random_device rd; random_numbers.seed(rd()); for (const auto& angle : thetas) { str_t common_case = check_common_cases(real_t(angle),eps); if(verbose) - cout << "Checking common cases..." << '\n'; + std::cout << "Checking common cases..." << '\n'; if(common_case != "") { if(details) - cout << "Angle is multiple of pi/4, answer is known exactly" << '\n'; + std::cout << "Angle is multiple of pi/4, answer is known exactly" << '\n'; if (check) { - cout << "Check flag = " << 1 << '\n'; + std::cout << "Check flag = " << 1 << '\n'; } - cout << common_case << endl; + std::cout << common_case << '\n'; return 0; } if(verbose) - cout << "No common cases found" << '\n'; + std::cout << "No common cases found" << '\n'; RzApproximation rz_approx; if (timer) { - auto start = chrono::steady_clock::now(); + auto start = std::chrono::steady_clock::now(); rz_approx = find_fast_rz_approximation( real_t(angle) * PI / real_t("-2"), eps); str_t op_str = synthesize(rz_approx.matrix(), s3_table); - auto end = chrono::steady_clock::now(); + auto end = std::chrono::steady_clock::now(); duration += - chrono::duration_cast(end - start) + std::chrono::duration_cast(end - start) .count(); } else { if (verbose) { - cout << "----\n"; - cout << "Finding approximation for angle = " << (angle) + std::cout << "----\n"; + std::cout << "Finding approximation for angle = " << (angle) << "..." << '\n'; } - if (more_verbose) { - rz_approx = verbose_find_fast_rz_approximation( - real_t(angle) * PI / real_t("-2"), eps); - } else { - rz_approx = find_fast_rz_approximation( - real_t(angle) * PI / real_t("-2"), eps); - } + rz_approx = find_fast_rz_approximation( + real_t(angle) * PI / real_t("-2"), eps); if (!rz_approx.solution_found()) { - cout << "No approximation found for RzApproximation. Try " + std::cout << "No approximation found for RzApproximation. Try " "changing factorization effort." << '\n'; return EXIT_FAILURE; } if (verbose) { - cout << "Approximation found. Synthesizing..." << '\n'; + std::cout << "Approximation found. Synthesizing..." << '\n'; } str_t op_str = synthesize(rz_approx.matrix(), s3_table); if (verbose) { - cout << "Synthesis complete." << '\n'; + std::cout << "Synthesis complete." << '\n'; } if (check) { - cout << "Check flag = " + std::cout << "Check flag = " << (rz_approx.matrix() == domega_matrix_from_str(full_simplify_str(op_str))) << '\n'; } if (details) { - real_t scale = pow(SQRT2, rz_approx.matrix().k()); - cout << "angle = " << angle << '\n'; - cout << rz_approx.matrix(); - cout << "u decimal value = " + real_t scale = gmpf::pow(SQRT2, rz_approx.matrix().k()); + std::cout << "angle = " << angle << '\n'; + std::cout << rz_approx.matrix(); + std::cout << "u decimal value = " << "(" << rz_approx.matrix().u().decimal().real() / scale << "," << rz_approx.matrix().u().decimal().imag() / scale << ")" << '\n'; - cout << "t decimal value = " + std::cout << "t decimal value = " << "(" << rz_approx.matrix().t().decimal().real() / scale << "," << rz_approx.matrix().t().decimal().imag() / scale << ")" << '\n'; - cout << "error = " << rz_approx.error() << '\n'; + std::cout << "error = " << rz_approx.error() << '\n'; str_t simplified = full_simplify_str(op_str); - string::difference_type n = + std::string::difference_type n = count(simplified.begin(), simplified.end(), 'T'); - cout << "T count = " << n << '\n'; - cout << "----" << '\n'; + std::cout << "T count = " << n << '\n'; + std::cout << "----" << '\n'; } for (auto& ch : full_simplify_str(op_str)) { - cout << ch << " "; + std::cout << ch << " "; } - cout << '\n'; + std::cout << '\n'; } } } if (timer) { - cout << "Duration = " << (static_cast(duration) / 1e6) << '\n'; + std::cout << "Duration = " << (static_cast(duration) / 1e6) << '\n'; } } diff --git a/unit_tests/tests/grid_synth/gmp_functions.cpp b/unit_tests/tests/grid_synth/gmp_functions.cpp index 66ffebd8..cc6e8854 100644 --- a/unit_tests/tests/grid_synth/gmp_functions.cpp +++ b/unit_tests/tests/grid_synth/gmp_functions.cpp @@ -3,7 +3,7 @@ #include "grid_synth/gmp_functions.hpp" using namespace staq; -using namespace grid_synth; +using namespace gmpf; TEST(GmpFunctions, min) { mpz_class x("10000012031201301030102301301023012301030103010301301023010101" @@ -11,24 +11,24 @@ TEST(GmpFunctions, min) { mpz_class y("12398198091845901809458349085918340581034850193845019839405810" "39485091834"); - EXPECT_TRUE(min(x, y) == x); - EXPECT_TRUE(min(y, x) == x); + EXPECT_TRUE(gmp_min(x, y) == x); + EXPECT_TRUE(gmp_min(y, x) == x); x = mpz_class("-13490123481023498091285094850984350980594810938450938409581" "8934514908249038290384"); y = mpz_class("-12910905890348590384509183045983019458091384509384509384509" "8130458103485013845094"); - EXPECT_TRUE(min(x, y) == x); - EXPECT_TRUE(min(y, x) == x); + EXPECT_TRUE(gmp_min(x, y) == x); + EXPECT_TRUE(gmp_min(y, x) == x); x = mpz_class("-51509183490598134905810934580913458091385039845093485093845" "0938450039485093850914"); y = mpz_class("394158094850934850913485903485019438590318501384509384509384" "5093840598103985091350"); - EXPECT_TRUE(min(x, y) == x); - EXPECT_TRUE(min(y, x) == x); + EXPECT_TRUE(gmp_min(x, y) == x); + EXPECT_TRUE(gmp_min(y, x) == x); } TEST(GmpFunctions, max) { @@ -37,50 +37,50 @@ TEST(GmpFunctions, max) { mpz_class y = mpz_class("12398198091845901809458349085918340581034850193845" "01983940581039485091834"); - EXPECT_TRUE(max(x, y) == y); - EXPECT_TRUE(max(y, x) == y); + EXPECT_TRUE(gmp_max(x, y) == y); + EXPECT_TRUE(gmp_max(y, x) == y); x = mpz_class("-13490123481023498091285094850984350980594810938450938409581" "8934514908249038290384"); y = mpz_class("-12910905890348590384509183045983019458091384509384509384509" "8130458103485013845094"); - EXPECT_TRUE(max(x, y) == y); - EXPECT_TRUE(max(y, x) == y); + EXPECT_TRUE(gmp_max(x, y) == y); + EXPECT_TRUE(gmp_max(y, x) == y); x = mpz_class("-51509183490598134905810934580913458091385039845093485093845" "0938450039485093850914"); y = mpz_class("394158094850934850913485903485019438590318501384509384509384" "5093840598103985091350"); - EXPECT_TRUE(max(x, y) == y); - EXPECT_TRUE(max(y, x) == y); + EXPECT_TRUE(gmp_max(x, y) == y); + EXPECT_TRUE(gmp_max(y, x) == y); } TEST(GmpFunctions, floor) { mpf_class x( "-1.123153451345634647367356735673567357635673567356735735673573573"); - EXPECT_TRUE(floor(x) == -2); + EXPECT_TRUE(gmp_floor(x) == -2); x = mpf_class("1." "213341545346456345647667356736573657356736573573567356735673" "5735735735"); - EXPECT_TRUE(floor(x) == 1); + EXPECT_TRUE(gmp_floor(x) == 1); } TEST(GmpFunctions, ceil) { mpf_class x( "-1.123153451345634647367356735673567357635673567356735735673573573"); - EXPECT_TRUE(ceil(x) == -1); + EXPECT_TRUE(gmp_ceil(x) == -1); x = mpf_class("1." "213341545346456345647667356736573657356736573573567356735673" "5735735735"); - EXPECT_TRUE(ceil(x) == 2); + EXPECT_TRUE(gmp_ceil(x) == 2); } TEST(GmpFunctions, round) { @@ -88,23 +88,23 @@ TEST(GmpFunctions, round) { mpf_class x( "-1.5452624566272736757567577777777777777777666666666666666666666"); - EXPECT_TRUE(round(x) == -2); + EXPECT_TRUE(gmp_round(x) == -2); x = mpf_class("-12001894518450983940850238599034852093502983049520923580385" "20385.04345234523523452345345"); - EXPECT_TRUE(round(x) == mpz_class("-120018945184509839408502385990348520935" + EXPECT_TRUE(gmp_round(x) == mpz_class("-120018945184509839408502385990348520935" "0298304952092358038520385")); x = mpf_class("1349218409." "231940295801948509438509834095810394859013845091384095810934" "580193845091385"); - EXPECT_TRUE(round(x) == mpz_class("1349218409")); + EXPECT_TRUE(gmp_round(x) == mpz_class("1349218409")); x = mpf_class( "4935810934580938." "66980981094819028419028409184309138240928409890902941390490284"); - EXPECT_TRUE(round(x) == mpz_class("4935810934580939")); + EXPECT_TRUE(gmp_round(x) == mpz_class("4935810934580939")); } diff --git a/unit_tests/tests/grid_synth/regions.cpp b/unit_tests/tests/grid_synth/regions.cpp index 7f0bc460..80f31e7b 100644 --- a/unit_tests/tests/grid_synth/regions.cpp +++ b/unit_tests/tests/grid_synth/regions.cpp @@ -126,10 +126,10 @@ TEST(EllipseTest, ZAndE) { real_t z = A.z(); real_t e = A.e(); - EXPECT_TRUE(abs(e * std::pow(LAMBDA.decimal().get_d(), -z.get_d()) - a) < - TOL); - EXPECT_TRUE(abs(e * std::pow(LAMBDA.decimal().get_d(), z.get_d()) - c) < - TOL); + EXPECT_TRUE(gmpf::gmp_abs(e * std::pow(LAMBDA.decimal().get_d(), -z.get_d()) - a) < + 1e-14); + EXPECT_TRUE(gmpf::gmp_abs(e * std::pow(LAMBDA.decimal().get_d(), z.get_d()) - c) < + 1e-14); } TEST(EllipseTest, Normalization) { diff --git a/unit_tests/tests/grid_synth/rings.cpp b/unit_tests/tests/grid_synth/rings.cpp index 33af539b..9c1817bf 100644 --- a/unit_tests/tests/grid_synth/rings.cpp +++ b/unit_tests/tests/grid_synth/rings.cpp @@ -76,7 +76,7 @@ TEST(ZSqrt2Inverse, LambdaInverse) { } TEST(ZOmegaConstructor, ConstructorEquality) { - double TOL = 1e-16; + MP_CONSTS.tol = 1e-16; ZOmega ztest(5, 6, -2, 1); diff --git a/unit_tests/tests/grid_synth/states.cpp b/unit_tests/tests/grid_synth/states.cpp index b0133e71..a49ebc64 100644 --- a/unit_tests/tests/grid_synth/states.cpp +++ b/unit_tests/tests/grid_synth/states.cpp @@ -1,13 +1,14 @@ #include "gtest/gtest.h" #include "grid_synth/states.hpp" +#include "grid_synth/gmp_functions.hpp" using namespace staq; using namespace grid_synth; TEST(ShiftState, ShiftToUnity) { Ellipse A(0, 0, 10, 5, 0.45 * PI); - Ellipse B(0, 0, pow(LAMBDA_INV.decimal(), 6), pow(LAMBDA_INV.decimal(), -6), + Ellipse B(0, 0, gmpf::pow(LAMBDA_INV.decimal(), 6), gmpf::pow(LAMBDA_INV.decimal(), -6), 0); real_t scaleA = A.normalize();