Skip to content

Commit

Permalink
Refactored gmp library, cleaned up implementation of constants
Browse files Browse the repository at this point in the history
  • Loading branch information
@jimlambert
jimlambert committed Sep 12, 2023
1 parent 9b1c98d commit 0f18d6b
Show file tree
Hide file tree
Showing 15 changed files with 298 additions and 501 deletions.
7 changes: 7 additions & 0 deletions include/grid_synth/complex.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include <iomanip>
#include <iostream>

#include "gmp_functions.hpp"

namespace staq {
namespace grid_synth {

Expand Down Expand Up @@ -142,6 +144,11 @@ inline complex<mpf_class> operator*(complex<double> x, complex<mpf_class> z) {
z.a() * x.b() + z.b() * x.a());
}


inline mpf_class abs(const complex<mpf_class>& z) {
return gmpf::sqrt(z.real() * z.real() + z.imag() * z.imag());
}

} // namespace grid_synth
} // namespace staq

Expand Down
187 changes: 53 additions & 134 deletions include/grid_synth/constants.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,27 +2,67 @@
#define CONSTANTS_HPP

#include <cmath>
#include <complex>
#include <iostream>

#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;
Expand All @@ -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

Expand Down
6 changes: 3 additions & 3 deletions include/grid_synth/diophantine_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)) {
Expand All @@ -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;
Expand Down
Loading

0 comments on commit 0f18d6b

Please sign in to comment.