From 8c10a34019bdc458ac547db2d6756d677dcef16e Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 18 Jul 2020 10:28:27 -0400 Subject: [PATCH 01/12] PSLQ [CI SKIP] --- doc/equations/algebraic.svg | 2 + doc/equations/basel_problem.svg | 2 + doc/equations/cfrac_accuracy.svg | 2 + doc/equations/khinchin_geometric.svg | 2 + doc/equations/khinchin_harmonic.svg | 2 + doc/equations/ln_basel_problem.svg | 2 + doc/equations/pslq_success.svg | 2 + doc/reference.qbk | 1 + doc/reference_pslq.qbk | 110 +++ example/pslq_demo.cpp | 86 +++ example/representations.cpp | 25 + example/representations/analyze.hpp | 116 +++ example/representations/cfrac_patterns.cpp | 167 +++++ example/representations/constants.cpp | 27 + example/representations/dottie.cpp | 25 + example/representations/golomb_dickman.cpp | 44 ++ example/representations/laplace_limit.cpp | 34 + example/representations/madelung.cpp | 26 + example/representations/nested_radical.cpp | 43 ++ .../representations/zeta_related_integral.cpp | 26 + include/boost/multiprecision/pslq.hpp | 665 ++++++++++++++++++ performance/pslq_performance.cpp | 69 ++ test/test_pslq.cpp | 70 ++ 23 files changed, 1548 insertions(+) create mode 100644 doc/equations/algebraic.svg create mode 100644 doc/equations/basel_problem.svg create mode 100644 doc/equations/cfrac_accuracy.svg create mode 100644 doc/equations/khinchin_geometric.svg create mode 100644 doc/equations/khinchin_harmonic.svg create mode 100644 doc/equations/ln_basel_problem.svg create mode 100644 doc/equations/pslq_success.svg create mode 100644 doc/reference_pslq.qbk create mode 100644 example/pslq_demo.cpp create mode 100644 example/representations.cpp create mode 100644 example/representations/analyze.hpp create mode 100644 example/representations/cfrac_patterns.cpp create mode 100644 example/representations/constants.cpp create mode 100644 example/representations/dottie.cpp create mode 100644 example/representations/golomb_dickman.cpp create mode 100644 example/representations/laplace_limit.cpp create mode 100644 example/representations/madelung.cpp create mode 100644 example/representations/nested_radical.cpp create mode 100644 example/representations/zeta_related_integral.cpp create mode 100644 include/boost/multiprecision/pslq.hpp create mode 100644 performance/pslq_performance.cpp create mode 100644 test/test_pslq.cpp diff --git a/doc/equations/algebraic.svg b/doc/equations/algebraic.svg new file mode 100644 index 000000000..cd7955b7a --- /dev/null +++ b/doc/equations/algebraic.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/basel_problem.svg b/doc/equations/basel_problem.svg new file mode 100644 index 000000000..f3d537cb3 --- /dev/null +++ b/doc/equations/basel_problem.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/cfrac_accuracy.svg b/doc/equations/cfrac_accuracy.svg new file mode 100644 index 000000000..6ffdf465d --- /dev/null +++ b/doc/equations/cfrac_accuracy.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/khinchin_geometric.svg b/doc/equations/khinchin_geometric.svg new file mode 100644 index 000000000..c8c033943 --- /dev/null +++ b/doc/equations/khinchin_geometric.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/khinchin_harmonic.svg b/doc/equations/khinchin_harmonic.svg new file mode 100644 index 000000000..be22f1fa1 --- /dev/null +++ b/doc/equations/khinchin_harmonic.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/ln_basel_problem.svg b/doc/equations/ln_basel_problem.svg new file mode 100644 index 000000000..fbf1ca5c3 --- /dev/null +++ b/doc/equations/ln_basel_problem.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/equations/pslq_success.svg b/doc/equations/pslq_success.svg new file mode 100644 index 000000000..241a7efa2 --- /dev/null +++ b/doc/equations/pslq_success.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/reference.qbk b/doc/reference.qbk index 3773479b9..4519091de 100644 --- a/doc/reference.qbk +++ b/doc/reference.qbk @@ -18,6 +18,7 @@ [include reference_mpfr_float.qbk] [include reference_cpp_bin_float.qbk] [include reference_cpp_dec_float.qbk] +[include reference_pslq.qbk] [include reference_internal_support.qbk] [include reference_backend_requirements.qbk] [include reference_header_structure.qbk] diff --git a/doc/reference_pslq.qbk b/doc/reference_pslq.qbk new file mode 100644 index 000000000..8711e14eb --- /dev/null +++ b/doc/reference_pslq.qbk @@ -0,0 +1,110 @@ +[/ + Copyright Nick Thompson, 2020 + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:pslq PSLQ] + + #include + namespace boost::multiprecision { + + template + std::string pslq(std::map const & dictionary, Real max_acceptable_solution_norm, std::ostream& = std::cout); + + template + bool is_algebraic(std::pair const & x, Real max_acceptable_solution_norm, std::ostream& = std::cout); + + template + std::string identify(std::pair const & x, Real max_acceptable_solution_norm, std::ostream& = std::cout); + } + +The PSLQ algorithm takes a list /x/ ∈ℝ[super n] of interesting constants, e.g., π, the Euler-Mascheroni constant γ, or the golden ration φ, and determines if there is an /integer relation/ between them. +By this we mean a list of integers /m/ ∈ℤ[super n], not all zero, such that /m⋅x/ = 0. +It is not obvious from this definition what the utility of the algorithm is, so let's examine an [@https://www.davidhbailey.com/dhbpapers/pslq-cse.pdf identity discovered by this method:] + +[$../equations/pslq_success.svg] + +Let's use the Boost implementation to rediscover this relation: + + using Real = cpp_bin_float_100; + std::map m; + m.emplace_back("0.0965512", "β"); + auto pi = boost::math::constants::pi(); + m.emplace_back(pow(pi,6), "π⁶"); + auto z3 = boost::math::constants::zeta_three(); + m.emplace_back(z3*z3, "ζ(3)²"); + std::string msg = pslq(m, 1e10); + + std::cout << msg << "\n"; + +We agonized over how to properly communicate this result, and in the end decided to just communicate the success or failure with a human-readable string. +Should you find a novel relation with `pslq`, it occupies the status of a major event-taking the time to read it should be perfectly acceptable. +Of course, if this is to be used in unanticipated ways, there is a lower-level API which gives machine readable status codes. + +If PSLQ fails to find a relation in a given dictionary, this is not proof of non-existence. +However, if a relation exists and PSLQ cannot find it, a lower bound on the 2-norm of /m/ is returned. +In this way, PSLQ can allow us to make statements like "π/e is probably not a rational number, because if π/e=p/q, then PSLQ shows that /p/[super 2] + /q/[super 2] > 10[super 200]." + +Conversely, if PSLQ finds an integer relation, that is not proof the relation is real. +In a large dictionary with each term computed to very low precision, it is easy to produce spurious relations. +This is why multiprecision is required: +If the dictionary has length /n/ and we want to find relations whose entries have height at most /d/, then we must compute the entries in our dictionary to /nd/ digits of precision. +/Try to make sure that every digit in the Real type you provide is correct!/ +If you only know the elements of the dictionary to double precision, do not upcast them to a multiprecision type. +This breaks various invariants used internally and makes it more likely that you will recover spurious relations. + +Even though PSLQ does not produce proofs, pessimism is unwarranted: PSLQ is a tool for discovery, not proof. +In addition, /verifying/ a suspected relation is much less computationally difficult than discovering it, so once a relation is found, we can compute the relevant terms of the dictionary to much higher precision to give additional evidence for the relation, and then search for a formal proof. + +Astute readers will notice that our example requires a very auspicious guess about the form the relation takes. +If we didn't know that β was somehow related to π⁶ and ζ(3)², we would've never found this relation. +This is a general property of the PSLQ algorithm: As its runtime scales quartically in the length of the dictionary, having a hunch about the form the relation will take is necessary. + +That said, simple modifications of the inputs allow us to recover different classes of relations. +Let's use the Basel problem as an example: + +[$../equations/basel_problem.svg] + +As stated, PSLQ cannot recover this relation unless π² is in the dictionary, which is kinda cheating. +However, we can take logarithms to obtain + +[$../equations/ln_basel_problem.svg] + +This shows how to recover multiplicative relations: Include the logarithms of small primes and logarithms of suspected terms in the dictionary. + +PSLQ can also determine if a number is algebraic. +A number is algebraic iff it is the root of a polynomial with integer coefficients. +If we suspect a number α is algebraic, we can apply PSLQ to the vector + +[$../equations/algebraic.svg] + +and if it is indeed algebraic with reasonable sized polynomial coefficients, then PSLQ will find it. +Communicating this idea has a somewhat different format than vanilla `pslq`, and hence we have provided a wrapper `is_algebraic` for precisely this purpose. + +Quartically scaling algorithms are rare, and users generally have little experience with them. +It is easy to construct a PSLQ dictionary that would only complete after a few thousand years. +Hence we have defaulted to printing a progress bar to the terminal. + + + [======> ] 10%, iteration 2000/20000, \u2016 m\u2016\u2082\u2265308, ETA:1.145 days + + +The progress bar is to be read as follows: The algorithm will complete in 20000 iterations. +At the time, it has proven that there are no integer relations with norm less that 308. +Under the assumption there are no integer relations, the algorithm will complete the demonstration that there are no integer solutions with norm less that `max_acceptable_solution_norm` in 1.145 days. +(If there are solutions, then they are found at essentially a random time.) +If the ETA is unacceptably high, note that you don't need to cancel the job. +You can stop it at any time and the 2-norm bound at the stop time still provides useful information. + +The progress bar is very useful, but in some cases it is not helpful (say, during unit tests). +We can silence the progress bar via by passing a `std::ostream` with the badbit set: + + std::ofstream ofs; + ofs.setstate(std::ios_base::badbit); + std::string s = boost::multiprecision::pslq(m, max_acceptable_norm_bound, ofs); + + + +[endsect] \ No newline at end of file diff --git a/example/pslq_demo.cpp b/example/pslq_demo.cpp new file mode 100644 index 000000000..4acdc6967 --- /dev/null +++ b/example/pslq_demo.cpp @@ -0,0 +1,86 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include +#include +#include +#include + +template +void ln2_plus_pi() +{ + using namespace boost::math::constants; + std::map m; + m.emplace(pi() + ln_two(), "(π+ln(2))"); + m.emplace(pi(), "π"); + m.emplace(ln_two(), "ln(2)"); + auto s = boost::multiprecision::pslq(m, 1e5); + if (!s.empty()) + { + std::cout << s << "\n"; + } + else + { + std::cout << "No relations found.\n"; + } +} + +template +void pi_root_two() +{ + using namespace boost::math::constants; + std::map m; + m.emplace(pi()/2, "π/2"); + m.emplace(root_two(), "√2"); + auto s = boost::multiprecision::pslq(m, 390); + if (!s.empty()) + { + std::cout << s << "\n"; + } + else + { + std::cout << "No relations found.\n"; + } +} + +template +void test_standard() +{ + std::map m = boost::multiprecision::small_pslq_dictionary(); + Real max_acceptable_norm_bound = 1e9; + std::ofstream ofs; + ofs.setstate(std::ios_base::badbit); + std::string s = boost::multiprecision::pslq(m, max_acceptable_norm_bound, ofs); + if (!s.empty()) { + std::cout << s << "\n"; + } else { + std::cout << "No relations found.\n"; + } +} + +int main() { + using Real = boost::multiprecision::number >; + test_standard(); + /*ln2_plus_pi(); + ln2_plus_pi(); + pi_root_two(); + pi_root_two(); + + std::map m = boost::multiprecision::standard_pslq_dictionary(); + Real max_acceptable_norm_bound = 1e9; + std::pair number{boost::math::lambert_w0(boost::math::constants::zeta_three()), "W(ζ(3))"}; + //std::ofstream ofs; + //ofs.setstate(std::ios_base::badbit); + std::string s = boost::multiprecision::identify(number, max_acceptable_norm_bound); + if (!s.empty()) + { + std::cout << s << "\n"; + } + else + { + std::cout << "No relations found.\n"; + }*/ +} \ No newline at end of file diff --git a/example/representations.cpp b/example/representations.cpp new file mode 100644 index 000000000..7d3e9aec6 --- /dev/null +++ b/example/representations.cpp @@ -0,0 +1,25 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2020 Nick Thompson. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include + +using boost::math::constants::gauss; + + +int main() { + using boost::multiprecision::mpfr_float; + using boost::multiprecision::representations; + mpfr_float::default_precision(100); + auto G = boost::math::constants::gauss(); + representations r(G, "G"); + + //using Real = boost::multiprecision::number>; + //auto G = gauss(); + //representations r(G, "G"); + std::cout << r << "\n"; +} diff --git a/example/representations/analyze.hpp b/example/representations/analyze.hpp new file mode 100644 index 000000000..44c0f564e --- /dev/null +++ b/example/representations/analyze.hpp @@ -0,0 +1,116 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#ifndef EXAMPLE_REPRESENTATIONS_ANALYZE_HPP +#define EXAMPLE_REPRESENTATIONS_ANALYZE_HPP +#include +#include +#include +#include +#include +#include +#include +#include + + +template +void analyze(Real x, std::string symbol) +{ + using std::abs; + using std::exp; + using std::log; + using boost::multiprecision::identify; + using boost::multiprecision::is_algebraic; + using boost::math::tools::centered_continued_fraction; + using boost::math::tools::simple_continued_fraction; + using boost::math::tools::luroth_expansion; + using boost::math::tools::engel_expansion; + using boost::multiprecision::checked_int1024_t; + constexpr int p = std::numeric_limits::max_digits10; + std::cout << std::setprecision(p); + if constexpr (p == 2147483647) { + std::cout << std::setprecision(x.backend().precision()); + } + + std::cout << symbol << " ≈ " << x << "\n"; + auto scf = simple_continued_fraction(x); + std::cout << " ≈ " << scf << "\n"; + auto ccf = centered_continued_fraction(x); + std::cout << " ≈ " << ccf << "\n"; + auto lur = luroth_expansion(x); + std::cout << " ≈ " << lur << "\n"; + auto eng = engel_expansion(x); + std::cout << " ≈ " << eng << "\n"; + std::cout << std::setprecision(11); + std::cout << "Empirical Khinchin geometric mean of simple cfrac : " << scf.khinchin_geometric_mean() << "\n"; + std::cout << "Expected Khinchin geometric mean if the value is 'random': 2.6854520010.\n"; + std::cout << "Empirical Khinchin harmonic mean of simple cfrac : " << scf.khinchin_harmonic_mean() << "\n"; + std::cout << "Expected Khinchin harmonic mean if the value is 'random' : 1.74540566240.\n"; + std::cout << "Empirical Khinchin geometric mean from centered continued fraction: " << ccf.khinchin_geometric_mean() << "\n"; + std::cout << "Expected Khinchin geometric mean of ccf if the value is 'random' : 5.454517244545585756966057724\n"; + std::cout << "Empirical digit geometric mean of Luroth expansion : " << lur.digit_geometric_mean() << "\n"; + std::cout << "Expected digit geometric mean if the value is 'random' : 2.200161058\n"; + + + Real expx = exp(x); + auto scf_exp = simple_continued_fraction(expx); + auto ccf_exp = centered_continued_fraction(expx); + auto lur_exp = luroth_expansion(expx); + std::cout << "exp(" << symbol << ")"; + std::cout << " ≈ " << scf_exp << "\n"; + std::cout << " ≈ " << ccf_exp << "\n"; + std::cout << " ≈ " << lur_exp << "\n"; + std::cout << std::setprecision(11); + std::cout << "Empirical Khinchin geometric mean of simple cfrac : " << scf_exp.khinchin_geometric_mean() << "\n"; + std::cout << "Expected Khinchin geometric mean if the value is 'random': 2.6854520010.\n"; + std::cout << "Empirical Khinchin harmonic mean of simple cfrac : " << scf_exp.khinchin_harmonic_mean() << "\n"; + std::cout << "Expected Khinchin harmonic mean if the value is 'random' : 1.74540566240.\n"; + std::cout << "Empirical Khinchin geometric mean from centered continued fraction: " << ccf_exp.khinchin_geometric_mean() << "\n"; + std::cout << "Expected Khinchin geometric mean of ccf if the value is 'random' : 5.454517244545585756966057724\n"; + std::cout << "Empirical digit geometric mean of Luroth expansion : " << lur_exp.digit_geometric_mean() << "\n"; + std::cout << "Expected digit geometric mean if the value is 'random' : 2.200161058\n"; + + Real lnx = log(x); + auto scf_log = simple_continued_fraction(lnx); + auto ccf_log = centered_continued_fraction(lnx); + auto lur_log = luroth_expansion(lnx); + std::cout << "ln(" << symbol << ")"; + std::cout << " ≈ " << scf_log << "\n"; + std::cout << " ≈ " << ccf_log << "\n"; + std::cout << " ≈ " << lur_log << "\n"; + std::cout << std::setprecision(11); + std::cout << "Empirical Khinchin geometric mean of simple cfrac : " << scf_log.khinchin_geometric_mean() << "\n"; + std::cout << "Expected Khinchin geometric mean if the value is 'random': 2.6854520010.\n"; + std::cout << "Empirical Khinchin harmonic mean of simple cfrac : " << scf_log.khinchin_harmonic_mean() << "\n"; + std::cout << "Expected Khinchin harmonic mean if the value is 'random' : 1.74540566240.\n"; + std::cout << "Empirical Khinchin geometric mean from centered continued fraction: " << ccf_log.khinchin_geometric_mean() << "\n"; + std::cout << "Expected Khinchin geometric mean of ccf if the value is 'random' : 5.454517244545585756966057724\n"; + std::cout << "Empirical digit geometric mean of Luroth expansion : " << lur_log.digit_geometric_mean() << "\n"; + std::cout << "Expected digit geometric mean if the value is 'random' : 2.200161058\n"; + + std::cout << "Is it algebraic?\n"; + std::string s = is_algebraic(x, symbol, Real(1e10)); + if (!s.empty()) + { + std::cout << s << "\n"; + } + else + { + std::cout << "It is probably not algebraic.\n"; + } + + std::cout << "Searching for closed forms.\n"; + s = identify(x, symbol, Real(1e10)); + if (!s.empty()) + { + std::cout << s << "\n"; + } + else + { + std::cout << "No relations found.\n"; + } +} +#endif diff --git a/example/representations/cfrac_patterns.cpp b/example/representations/cfrac_patterns.cpp new file mode 100644 index 000000000..a06adb834 --- /dev/null +++ b/example/representations/cfrac_patterns.cpp @@ -0,0 +1,167 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include "analyze.hpp" +#include +#include +#include +using boost::multiprecision::mpfr_float; +using boost::math::tools::continued_fraction_a; + +// This example comes from Steven R. Finch's book "Mathematical Constants", Section 6.2. +// "No one knows the outcome if we instead repeat each denominator in c_1, +// although numerically we find c_2 = 0.5851972651...." +// Mathematica notation: FromContinuedFraction[{0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10}] +template +class repeated_denom_generator +{ +public: + repeated_denom_generator() { i = 0; even = true;} + typedef T result_type; + + result_type operator()() + { + if (even) + { + i += 1; + even = false; + } + else + { + even = true; + } + return i; + } +private: + T i; + bool even; +}; +template +Real repeated_denom() +{ + auto gen = repeated_denom_generator(); + Real res = continued_fraction_a(gen, std::numeric_limits::epsilon()); + return res; +} + +// FromContinuedFraction[{0,1,1,4,9,16,25,36,49,64,81,100,121,144}] = 0.554226 +template +class denom_squared_generator +{ +public: + denom_squared_generator() { i = 1; first = true; } + typedef T result_type; + + result_type operator()() + { + if (first) + { + first = false; + return 1; + } + T res = i*i; + ++i; + return res; + } +private: + T i; + bool first; +}; + +template +Real denom_squared() +{ + auto gen = denom_squared_generator(); + Real res = continued_fraction_a(gen, std::numeric_limits::epsilon()); + return res; +} + + +// FromContinuedFraction[{0,1,1,2^3,3^3, 4^3,...}] = 0.52984 +template +class denom_cubed_generator +{ +public: + denom_cubed_generator() { i = 1; first = true; } + typedef T result_type; + + result_type operator()() + { + if (first) + { + first = false; + return 1; + } + T res = i*i*i; + ++i; + return res; + } +private: + T i; + bool first; +}; + +template +Real denom_cubed() +{ + auto gen = denom_cubed_generator(); + Real res = continued_fraction_a(gen, std::numeric_limits::epsilon()); + return res; +} + +// FromContinuedFraction[{0, 1, 1, 2^4, 3^4, 4^4,...}] = 0.51514 +template +class denom_quarted_generator +{ +public: + denom_quarted_generator() { i = 1; first = true; } + typedef T result_type; + + result_type operator()() + { + if (first) + { + first = false; + return 1; + } + T res = i*i*i*i; + ++i; + return res; + } +private: + T i; + bool first; +}; + +template +Real denom_quarted() +{ + auto gen = denom_quarted_generator(); + Real res = continued_fraction_a(gen, std::numeric_limits::epsilon()); + return res; +} + +int main() +{ + using Real = mpfr_float; + mpfr_float::default_precision(1000); + Real x = denom_quarted(); + std::string symbol = "[0; 1, 1⁴, 2⁴, 3⁴, 4⁴, 5⁴, ...]"; + analyze(x, symbol); + + x = denom_cubed(); + symbol = "[0; 1, 1³, 2³, 3³, 4³, 5³, ...]"; + analyze(x, symbol); + + x = denom_squared(); + symbol = "[0; 1, 1², 2², 3², 4², 5², ...]"; + analyze(x, symbol); + + + x = repeated_denom(); + symbol = "1|/|1 + 1|/|1 + 1|/|2 + 1|/|2 + "; + analyze(x, symbol); +} \ No newline at end of file diff --git a/example/representations/constants.cpp b/example/representations/constants.cpp new file mode 100644 index 000000000..be03c5b93 --- /dev/null +++ b/example/representations/constants.cpp @@ -0,0 +1,27 @@ +#include +#include "analyze.hpp" +#include + +using namespace boost::math::constants; +int main() { + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(1000); + using Real = mpfr_float; + + analyze(pi(), "π"); + analyze(euler(), "γ"); + analyze(khinchin(), "K0"); + analyze(glaisher(), "glaisher"); + analyze(plastic(), "plastic"); + analyze(gauss(), "gauss"); + analyze(first_feigenbaum(), "feigenbaum"); + analyze(rayleigh_kurtosis_excess(), "rayleigh_excess_kurtosis"); + analyze(rayleigh_skewness(), "rayleigh_skewness"); + analyze(extreme_value_skewness(), "extreme_value_skewness"); + analyze(catalan(), "catalan"); + analyze(zeta_three(), "zeta_three"); + analyze(sin_one(), "sin(1)"); + analyze(cos_one(), "cos(1)"); + analyze(e_pow_pi(), "e^pi"); + +} diff --git a/example/representations/dottie.cpp b/example/representations/dottie.cpp new file mode 100644 index 000000000..358b22f13 --- /dev/null +++ b/example/representations/dottie.cpp @@ -0,0 +1,25 @@ +#include +#include "analyze.hpp" + +template +Real dottie() +{ + using std::cos; + using std::abs; + using std::sin; + Real x{".739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246"}; + Real residual = cos(x) - x; + do { + x += residual/(sin(x)+1); + residual = cos(x) - x; + } while(abs(residual) > std::numeric_limits::epsilon()); + return x; +} + +int main() { + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(1300); + using Real = mpfr_float; + Real dot = dottie(); + analyze(dot, "dottie"); +} diff --git a/example/representations/golomb_dickman.cpp b/example/representations/golomb_dickman.cpp new file mode 100644 index 000000000..cfbf0ee7a --- /dev/null +++ b/example/representations/golomb_dickman.cpp @@ -0,0 +1,44 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include "analyze.hpp" +#include +#include +#include +#include +using boost::multiprecision::mpfr_float; +using boost::math::quadrature::tanh_sinh; +using std::exp; +using std::log; +using boost::math::expint; +using boost::math::constants::euler; + +template +Real golomb_dickman() +{ + auto integrator = tanh_sinh(18); + auto f = [](Real t) { return exp(expint(log(t)));}; + Real Q = integrator.integrate(f, Real(0), Real(1)); + return Q; +} + +int main() +{ + using Real = mpfr_float; + int p = 500; + mpfr_float::default_precision(p); + std::cout << std::setprecision(p); + Real lambda = golomb_dickman(); + std::cout << "Expected from OEIS = 0.624329988543550870992936383100837244179642620180529286973551902495638088855113254462460276195539868869\n"; + std::cout << "Computed = " << lambda << "\n"; + std::string symbol = "λ"; + analyze(lambda, symbol); + + lambda *= exp(-euler()); + symbol = "λexp(-𝛾)"; + analyze(lambda, symbol); + +} \ No newline at end of file diff --git a/example/representations/laplace_limit.cpp b/example/representations/laplace_limit.cpp new file mode 100644 index 000000000..5e820eb7d --- /dev/null +++ b/example/representations/laplace_limit.cpp @@ -0,0 +1,34 @@ +#include +#include "analyze.hpp" + +template +Real laplace_limit() +{ + using std::cos; + using std::abs; + using std::sin; + using std::sqrt; + using std::exp; + Real x{"0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168"}; + Real tmp = sqrt(1+x*x); + Real etmp = exp(tmp); + Real residual = x*exp(tmp) - 1 - tmp; + Real df = etmp -x/tmp + etmp*x*x/tmp; + do { + x -= residual/df; + tmp = sqrt(1+x*x); + etmp = exp(tmp); + residual = x*exp(tmp) - 1 - tmp; + df = etmp -x/tmp + etmp*x*x/tmp; + } while(abs(residual) > std::numeric_limits::epsilon()); + std::cout << "Residual = " << residual << "\n"; + return x; +} + +int main() { + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(1300); + using Real = mpfr_float; + Real ll = laplace_limit(); + analyze(ll, "laplace_limit"); +} diff --git a/example/representations/madelung.cpp b/example/representations/madelung.cpp new file mode 100644 index 000000000..8e1d768c0 --- /dev/null +++ b/example/representations/madelung.cpp @@ -0,0 +1,26 @@ +#include +#include "analyze.hpp" + +template +Real madelung() +{ + // Computed in Mathematica using: + /* + Quiet[12 Pi (Sech[Pi/Sqrt[2]]^2 + + NSum[Sum[ + Sech[Pi Norm[2 v + 1]/2]^2, {v, + FrobeniusSolve[{1, 1}, Round[m]]}, Method -> "Procedural"], {m, + 1, Infinity}, Compiled -> False, Method -> "WynnEpsilon", + NSumTerms -> 200, WorkingPrecision -> 500])] + */ + Real x{"1.7475645946331821906362120355443974034851614366247417581528253507650406235327611798907583626946078899308325815387537105932820299441838280130369330021565993632823766071722975686592380371672038104106034214556064382777786832173132243697558773426250474787821285086056791668167573992447684129703678251857628109371313372076707193197424971581157230969923096692739496577811072226715205474090115068915716583082820050184892117803134673122964985829"}; + return x; +} + +int main() { + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(400); + using Real = mpfr_float; + Real mad = madelung(); + analyze(mad, "M"); +} diff --git a/example/representations/nested_radical.cpp b/example/representations/nested_radical.cpp new file mode 100644 index 000000000..e3791407c --- /dev/null +++ b/example/representations/nested_radical.cpp @@ -0,0 +1,43 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include "analyze.hpp" +#include +using boost::multiprecision::mpfr_float; + +// This example comes from Steven R. Finch's book "Mathematical Constants", Section 1.2.1. +// Ideally, we'd know how many terms to use just on the precision of the Real number. +// But we don't know that! +template +Real nested_123(int64_t max_n) +{ + assert(max_n > 1); + using std::sqrt; + Real r = sqrt(max_n); + for (int64_t n = max_n - 1; n >= 1; --n) { + r = sqrt(n + r); + } + return r; +} + +int main() +{ + using Real = mpfr_float; + int p = 1000; + mpfr_float::default_precision(p); + int n = 30000; + Real x1 = nested_123(n/2); + std::cout << std::setprecision(p); + Real x2 = nested_123(n); + std::string symbol = "√(1 + √(2 + √(3 +...)"; + Real delta = abs(x1 - x2); + if (delta > 0) { + std::cerr << std::scientific << "delta = " << abs(x2 - x1) << "\n"; + std::cerr << "n must be increased to get all the bits correct.\n"; + return 1; + } + analyze(x2, symbol); +} \ No newline at end of file diff --git a/example/representations/zeta_related_integral.cpp b/example/representations/zeta_related_integral.cpp new file mode 100644 index 000000000..302a91c9a --- /dev/null +++ b/example/representations/zeta_related_integral.cpp @@ -0,0 +1,26 @@ +// Copyright Nick Thompson, 2020 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) +#include "analyze.hpp" +#include +#include +#include + +using std::log; +using std::sin; +using std::abs; +using boost::math::quadrature::tanh_sinh; +using boost::multiprecision::mpfr_float; +using boost::math::constants::pi; +using boost::math::constants::zeta_three; + +int main() { + using Real = mpfr_float; + mpfr_float::default_precision(1000); + Real x = -(7*zeta_three() - pi()*pi()*log(static_cast(4)))/16;; + std::string symbol = "(-7zeta(3) + pi^2ln(4))/16"; + analyze(x, symbol); + return 0; +} \ No newline at end of file diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp new file mode 100644 index 000000000..516032261 --- /dev/null +++ b/include/boost/multiprecision/pslq.hpp @@ -0,0 +1,665 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +// Mathematica has an implementation of PSLQ which has the following interface: +// FindIntegerNullVector[{E, Pi}, 100000] +// FindIntegerNullVector::norel: There is no integer null vector for {E,\[Pi]} with norm less than or equal to 100000. +// Or: +// FindIntegerNullVector[{E, \[Pi]}] +// FindIntegerNullVector::rnfu: FindIntegerNullVector has not found an integer null vector for {E,\[Pi]}. +// I don't like this, because it should default to telling us the norm, as it's coproduced by the computation. + +// Maple's Interface: +// with(IntegerRelations) +// v := [1.57079..., 1.4142135] +// u := PSLQ(v) +// u:= [-25474, 56589] +// Maple's interface is in fact worse, because it gives the wrong answer, instead of recognizing the precision provided. + +// David Bailey's interface in tpslqm2.f90 in https://www.davidhbailey.com/dhbsoftware/ in mpfun-fort-v19.tar.gz +// subroutine pslqm2(idb, n nwds, rb, eps, x, iq, r) +// idb is debug level +// n is the length of input vector and output relation r. +// nwds if the full precision level in words. +// rb is the log10 os max size Euclidean norm of relation +// eps tolerance for full precision relation detection. +// x input vector +// iq output flag: 0 (unsuccessful), 1 successful. +// r output integer relation vector, if successful. + +#ifndef BOOST_MULTIPRECISION_PSLQ_HPP +#define BOOST_MULTIPRECISION_PSLQ_HPP +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#if defined __has_include +# if __has_include () +# include +# else + #error This file has a dependency on Eigen; you must have Eigen (http://eigen.tuxfamily.org/index.php?title=Main_Page) in your include paths. +# endif +#endif + +#if defined(__APPLE__) || defined(__linux__) +#include +#include +#elif defined(_WIN32) || defined(__CYGWIN__) +#include +#endif + +namespace boost::multiprecision { + +// For debugging and unit testing: +template +auto small_pslq_dictionary() { + using std::sqrt; + using namespace boost::math::constants; + std::map m; + m.emplace(pi(), "π"); + m.emplace(e(), "e"); + m.emplace(root_two(), "√2"); + m.emplace(ln_two(), "ln(2)"); + Real euler_ = euler(); + m.emplace(euler_, "γ"); + m.emplace(euler_*euler_, "γ²"); + m.emplace(euler_*euler_*euler_, "γ³"); + m.emplace(1/euler_, "1/γ"); + m.emplace(1/(euler_*euler_), "1/γ²"); + m.emplace(1/(euler_*euler_*euler_), "1/γ³"); + m.emplace(-log(euler_), "-ln(γ)"); + m.emplace(exp(euler_), "exp(γ)"); + Real zeta_three_ = zeta_three(); + m.emplace(sqrt(zeta_three_), "√ζ(3)"); + m.emplace(zeta_three_, "ζ(3)"); + m.emplace(1/zeta_three_, "1/ζ(3)"); + m.emplace(1/(zeta_three_*zeta_three_), "1/ζ(3)²"); + m.emplace(1/(zeta_three_*zeta_three_*zeta_three_), "1/ζ(3)³"); + m.emplace(log(zeta_three_), "ln(ζ(3))"); + m.emplace(exp(zeta_three_), "exp(ζ(3))"); + m.emplace(zeta_three_*zeta_three_, "ζ(3)²"); + m.emplace(zeta_three_*zeta_three_*zeta_three_, "ζ(3)³"); + m.emplace(pow(zeta_three_, 4), "ζ(3)⁴"); + + + return m; +} + +template +auto standard_pslq_dictionary() { + using std::sqrt; + using std::log; + using std::exp; + using std::pow; + using std::cbrt; + using namespace boost::math::constants; + std::map m; + Real euler_ = euler(); + m.emplace(euler_, "γ"); + m.emplace(euler_*euler_, "γ²"); + m.emplace(euler_*euler_*euler_, "γ³"); + m.emplace(1/euler_, "1/γ"); + m.emplace(1/(euler_*euler_), "1/γ²"); + m.emplace(1/(euler_*euler_*euler_), "1/γ³"); + m.emplace(-log(euler_), "-ln(γ)"); + m.emplace(exp(euler_), "exp(γ)"); + Real zeta_three_ = zeta_three(); + m.emplace(sqrt(zeta_three_), "√ζ(3)"); + m.emplace(zeta_three_, "ζ(3)"); + m.emplace(1/zeta_three_, "1/ζ(3)"); + m.emplace(1/(zeta_three_*zeta_three_), "1/ζ(3)²"); + m.emplace(1/(zeta_three_*zeta_three_*zeta_three_), "1/ζ(3)³"); + m.emplace(log(zeta_three_), "ln(ζ(3))"); + m.emplace(exp(zeta_three_), "exp(ζ(3))"); + m.emplace(zeta_three_*zeta_three_, "ζ(3)²"); + m.emplace(zeta_three_*zeta_three_*zeta_three_, "ζ(3)³"); + m.emplace(pow(zeta_three_, 4), "ζ(3)⁴"); + + auto pi_ = pi(); + m.emplace(pi_, "π"); + m.emplace(1/pi_, "1/π"); + m.emplace(1/(pi_*pi_), "1/π²"); + m.emplace(sqrt(pi_), "√π"); + m.emplace(cbrt(pi_), "∛π"); + m.emplace(log(pi_), "ln(π)"); + m.emplace(pi_*pi_, "π²"); + m.emplace(pi_*pi_*pi_, "π³"); + + Real e_ = e(); + m.emplace(e_, "e"); + m.emplace(sqrt(e_), "√e"); + m.emplace(root_two(), "√2"); + m.emplace(cbrt(static_cast(2)), "∛2"); + m.emplace(cbrt(static_cast(3)), "∛3"); + m.emplace(root_three(), "√3"); + m.emplace(sqrt(static_cast(5)), "√5"); + m.emplace(sqrt(static_cast(7)), "√7"); + m.emplace(sqrt(static_cast(11)), "√11"); + + // φ is linearly dependent on √5; its logarithm is not. + m.emplace(log(phi()), "ln(φ)"); + m.emplace(exp(phi()), "exp(φ)"); + Real catalan_ = catalan(); + m.emplace(catalan_, "G"); + m.emplace(catalan_*catalan_, "G²"); + m.emplace(1/catalan_, "1/G"); + m.emplace(-log(catalan_), "-ln(G)"); + m.emplace(exp(catalan_), "exp(G)"); + m.emplace(sqrt(catalan_), "√G"); + + Real glaisher_ = glaisher(); + m.emplace(glaisher_, "A"); + m.emplace(glaisher_*glaisher_, "A²"); + m.emplace(1/glaisher_, "1/A"); + m.emplace(log(glaisher_), "ln(A)"); + m.emplace(exp(glaisher_), "exp(A)"); + Real khinchin_ = khinchin(); + m.emplace(khinchin_, "K₀"); + m.emplace(log(khinchin_), "ln(K₀)"); + m.emplace(exp(khinchin_), "exp(K₀)"); + m.emplace(1/khinchin_, "1/K₀"); + m.emplace(khinchin_*khinchin_, "K₀²"); + m.emplace(log(static_cast(1) + sqrt(static_cast(2))), "ln(1+√2)"); + // To recover multiplicative relations we need the logarithms of small primes. + Real ln2 = log(static_cast(2)); + m.emplace(ln2, "ln(2)"); + m.emplace(-log(ln2), "-ln(ln(2))"); + m.emplace(log(static_cast(3)), "ln(3)"); + m.emplace(log(static_cast(5)), "ln(5)"); + m.emplace(log(static_cast(7)), "ln(7)"); + m.emplace(log(static_cast(11)), "ln(11)"); + m.emplace(log(static_cast(13)), "ln(13)"); + m.emplace(log(static_cast(17)), "ln(17)"); + m.emplace(log(static_cast(19)), "ln(19)"); + // Omega constant = Lambert-W function evaluated at 1: + Real Omega_ = boost::math::lambert_w0(static_cast(1)); + m.emplace(Omega_, "Ω"); + m.emplace(Omega_*Omega_, "Ω²"); + m.emplace(1/Omega_, "1/Ω"); + + // Should we add the Madelung constant? + + // Now we add a few that will allow recovery of relations from 'Mathematical Constants' by Steven Finch. + // Is this a sensible way to go about this? How else should a standard dictionary be defined? + + // Looking at Mathematical Constants, 1.5.2, we need these to recover the relations stated there: + m.emplace(ln2*sqrt(pi_), "ln(2)√π"); + m.emplace(euler_*sqrt(pi_), "γ√π"); + // Section 1.5.3: + m.emplace(root_three()*pi_, "π√3"); + m.emplace(catalan_*pi_, "πG"); + + return m; +} + +// anonymous namespace: +namespace { + +class progress { + +public: + +progress(std::ostream & os) : os_{os} +{ + start_ = std::chrono::steady_clock::now(); +#if defined(__APPLE__) || defined(__linux__) + struct winsize size_; + ioctl(STDOUT_FILENO, TIOCGWINSZ, &size_); + bar_width_ = size_.ws_col; +#elif defined(_WIN32) + // From: https://stackoverflow.com/a/23370070/904050 + CONSOLE_SCREEN_BUFFER_INFO csbi; + GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi); + bar_width_ = csbi.srWindow.Right - csbi.srWindow.Left + 1; +#else + bar_width_ = 90; +#endif + bar_width_ -= 70; + os_ << std::setprecision(2); +} + +void display_progress(int64_t iteration, int64_t max_iterations, double norm_bound) +{ + double progress = iteration/double(max_iterations); + auto now = std::chrono::steady_clock::now(); + std::chrono::duration elapsed_milliseconds = now - start_; + // ETA: + double eta = (1-progress)*elapsed_milliseconds.count()/(1000*progress); + int pos = bar_width_ * progress; + os_ << "\033[0;32m["; + for (int i = 0; i < bar_width_; ++i) { + if (i < pos) os_ << "="; + else if (i == pos) os_ << ">"; + else os_ << " "; + } + os_ << "] " + << iteration * 100.0/max_iterations + << "%, iteration " << iteration << "/" << max_iterations; + os_ << ", ‖m‖₂≥" << norm_bound << ", ETA:"; + if (eta < 3600) { + os_ << eta << " s\r"; + } + else if (eta < 3600*24) { + os_ << eta/3600 << " hr\r"; + } + else if (eta < 3600*24*7) { + os_ << eta/(3600*24) << " days\r"; + } + else if (eta < 3600*24*7*4) { + os_ << eta/(3600*24*7) << " wks\r"; + } + else if (eta < 3.154e+7) { + os_ << eta/(2.628e+6) << " months\r"; + } + else { + os_ << eta/3.154e+7 << " years\r"; + } + os_.flush(); +} + +~progress() { + os_ << "\033[39m\n"; +} + +private: + std::chrono::steady_clock::time_point start_; + int bar_width_; + std::ostream & os_; +}; + +} + +// The PSLQ algorithm; partial sum of squares, lower trapezoidal decomposition. +// See: https://www.davidhbailey.com/dhbpapers/cpslq.pdf, section 3. +template +std::vector> pslq(std::vector & x, Real max_acceptable_norm_bound, Real gamma, std::ostream& os = std::cout) { + std::vector> relation; + /*if (!std::is_sorted(x.begin(), x.end())) { + std::cerr << "Elements must be sorted in increasing order.\n"; + return relation; + }*/ + + std::sort(x.begin(), x.end()); + using std::sqrt; + if (gamma <= 2/sqrt(3)) { + std::cerr << "γ > 2/√3 is required\n"; + return relation; + } + Real tmp = 1/Real(4) + 1/(gamma*gamma); + Real tau = 1/sqrt(tmp); + if (tau <= 1 || tau >= 2) { + std::cerr << "τ ∈ (1, 2) is required.\n"; + return relation; + } + + if (x.size() < 2) { + std::cerr << "At least two values are required to find an integer relation.\n"; + return relation; + } + + for (auto & t : x) { + if (t == 0) { + std::cerr << "Zero in the dictionary gives trivial relations.\n"; + return relation; + } + if (t < 0) { + std::cerr << "The algorithm is reflection invariant, so negative values in the dictionary should be removed.\n"; + return relation; + } + } + + // Are we computing too many square roots??? Should we use s instead? + std::vector s_sq(x.size()); + s_sq.back() = x.back()*x.back(); + int64_t n = x.size(); + for (int64_t i = n - 2; i >= 0; --i) { + s_sq[i] = s_sq[i+1] + x[i]*x[i]; + } + + using std::pow; + if (max_acceptable_norm_bound*max_acceptable_norm_bound*s_sq[0] > 1/std::numeric_limits::epsilon()) { + std::cerr << "The maximum acceptable norm bound " << max_acceptable_norm_bound << " is too large; spurious relations will be recovered.\n"; + std::cerr << "Either reduce the norm bound, or increase the precision of the variables.\n"; + std::cerr << "At this precision, your norm bound cannot exceed " << 1/sqrt(s_sq[0]*std::numeric_limits::epsilon()) << "\n"; + return relation; + } + Eigen::Matrix H(n, n-1); + for (int64_t i = 0; i < n - 1; ++i) { + for (int64_t j = 0; j < n - 1; ++j) { + if (i < j) { + H(i,j) = 0; + } + else if (i == j) { + H(i,i) = sqrt(s_sq[i+1]/s_sq[i]); + } + else { + // i > j: + H(i,j) = -x[i]*x[j]/sqrt(s_sq[j]*s_sq[j+1]); + } + } + } + for (int64_t j = 0; j < n - 1; ++j) { + H(n-1, j) = -x[n-1]*x[j]/sqrt(s_sq[j]*s_sq[j+1]); + } + + using std::ceil; + int64_t expected_iterations = ceil(boost::math::binomial_coefficient(n, 2)*log(pow(gamma, n-1)*max_acceptable_norm_bound)/log(tau)); + //std::cout << "Expected number of iterations = " << expected_iterations << "\n"; + // This validates that H is indeed lower trapezoidal, + // but that's trival and verbose: + //std::cout << "H = \n"; + //std::cout << H << "\n"; + + // Validate the conditions of Lemma 1 in the referenced paper: + // These tests should eventually be removed once we're confident that the code is correct. + auto Hnorm_sq = H.squaredNorm(); + if (abs(Hnorm_sq/(n-1) - 1) > sqrt(std::numeric_limits::epsilon())) { + std::cerr << "‖Hₓ‖² ≠ n - 1. Hence Lemma 1.ii of the reference has numerically failed; this is a bug.\n"; + return relation; + } + + // Notations now follows https://www.davidhbailey.com/dhbpapers/pslq-cse.pdf + Eigen::Matrix y(x.size()); + for (int64_t i = 0; i < n; ++i) { + y[i] = x[i]/sqrt(s_sq[0]); + } + + // Now check y: + for (int64_t i = 0; i < n; ++i){ + if (abs(y[i]) < std::numeric_limits::epsilon()) { + std::cerr << "Element y[" << i << "] = " << y[i] << " is too small; more precision is required.\n"; + return relation; + } + } + for (int64_t i = 1; i < n; ++i){ + // using the sorted assumption here: + if (abs(boost::math::float_distance(y[i], y[i-1])) <= 2) { + std::cerr << "Elements y[" << i << "] = " << y[i] << " and " << y << "[" << i + 1 << "] = " << y[i+1] << " are too close together.\n"; + return relation; + } + } + + auto v = y.transpose()*H; + for (int64_t i = 0; i < n - 1; ++i) { + if (abs(v[i])/(n-1) > sqrt(std::numeric_limits::epsilon())) { + std::cerr << "xᵀHₓ ≠ 0; Lemma 1.iii of the reference cpslq has numerically failed; this is a bug.\n"; + return relation; + } + } + using std::round; + + // Initialize: + // "1. Set the nxn matrices A and B to the identity." + Eigen::Matrix A = Eigen::Matrix::Identity(n, n); + Eigen::Matrix B = Eigen::Matrix::Identity(n, n); + for (int64_t i = 1; i < n; ++i) { + for (int64_t j = i - 1; j >= 0; --j) { + Real t = round(H(i,j)/H(j,j)); + int64_t tint = static_cast(t); + // This happens a lot because x_0 < x_1 < ...! + // Sort them in decreasing order and it almost never happens. + if (tint == 0) + { + continue; + } + for (int64_t k = 0; k <= j; ++k) + { + H(i,k) = H(i,k) - t*H(j,k); + } + for (int64_t k = 0; k < n; ++k) { + A(i,k) = A(i,k) - tint*A(j,k); + B(k,j) = B(k,j) + tint*B(k,i); + } + y[j] += t*y[i]; + } + } + //std::cout << "A, post-reduction = \n" << A << "\n"; + //std::cout << "B, post-reduction = \n" << B << "\n"; + //std::cout << "A*B, post-reduction = \n" << A*B << "\n"; + //std::cout << "H, post-reduction:\n" << H << "\n"; + Real max_coeff = 0; + for (int64_t i = 0; i < n - 1; ++i) { + if (abs(H(i,i)) > max_coeff) { + max_coeff = abs(H(i,i)); + } + } + Real norm_bound = 1/max_coeff; + Real last_norm_bound = norm_bound; + int64_t iteration = 0; + auto prog = progress(os); + while (norm_bound < max_acceptable_norm_bound) + { + // "1. Select m such that γ^{i+1}|H_ii| is maximal when i = m": + // (note my C indexing translated from DHB's Fortran indexing) + Real gammai = gamma; + Real max_term = 0; + int64_t m = -1; + for (int i = 0; i < n - 1; ++i) { + Real term = gammai*abs(H(i,i)); + if (term > max_term) { + max_term = term; + m = i; + } + gammai *= gamma; + } + // "2. Exchange the entries of y indexed m and m + 1" + if (m == n - 1) { + std::cerr << "OMG: m = n- 1, swap gonna segfault.\n"; + return relation; + } + if (m < 0) { + std::cerr << "OMG: m = - 1, swap gonna segfault.\n"; + return relation; + } + std::swap(y[m], y[m+1]); + // Swap the corresponding rows of A and H: + A.row(m).swap(A.row(m+1)); + H.row(m).swap(H.row(m+1)); + // Swap the corresponding columns of B: + B.col(m).swap(B.col(m+1)); + + //std::cout << "H, post-swap = \n" << H << "\n"; + // "3. Remove the corner on H diagonal:" + if (m < n - 2) { + Real t0 = H(m,m)*H(m,m) + H(m, m+1)*H(m, m+1); + using boost::math::rsqrt; + t0 = rsqrt(t0); + Real t1 = H(m,m)*t0; + Real t2 = H(m,m+1)*t0; + for (int64_t i = m; i < n; ++i) { + Real t3 = H(i,m); + Real t4 = H(i, m+1); + H(i,m) = t1*t3 + t2*t4; + H(i,m+1) = -t2*t3 + t1*t4; + } + //std::cout << "H, post corner reduce = \n" << H << "\n"; + } + + // "4. Reduce H:" + for (int64_t i = m+1; i < n; ++i) { + for (int64_t j = std::min(i-1, m+1); j >= 0; --j) { + Real t = round(H(i,j)/H(j,j)); + int64_t tint = static_cast(t); + if (tint == 0) + { + continue; + } + y[j] += t*y[i]; + for (int64_t k = 0; k <= j; ++k) { + H(i,k) = H(i,k) - t*H(j,k); + } + for (int64_t k = 0; k < n; ++k) { + A(i,k) = A(i,k) - tint*A(j,k); + B(k,j) = B(k,j) + tint*B(k,i); + } + } + } + + // Look for a solution: + for (int64_t i = 0; i < n; ++i) { + if (abs(y[i]) < pow(std::numeric_limits::epsilon(), Real(15)/Real(16))) { + auto bcol = B.col(i); + Real residual = 0; + Real absum = 0; + for (int64_t j = 0; j < n; ++j) { + residual += bcol[j]*x[j]; + absum += abs(bcol[j]*x[j]); + } + Real tolerable_residual = 16*std::numeric_limits::epsilon()*absum; + if (abs(residual) > tolerable_residual) + { + std::cout << std::endl; + std::cerr << "\033[31m"; + std::cerr << __FILE__ << ":" << __LINE__ << "\n\tBug in PSLQ: Found a relation with a large residual.\n"; + std::cerr << "\tThe residual " << abs(residual) << " is larger than the tolerable residual " << tolerable_residual << "\n"; + std::cerr << "\tThis is either a bug, or your constants are not specified to the full accuracy of the floating point type " << boost::core::demangle(typeid(Real).name()) << ".\n"; + } + + for (int64_t j = 0; j < n; ++j) + { + if (bcol[j] != 0) + { + relation.push_back({bcol[j], x[j]}); + } + } + return relation; + } + } + + max_coeff = 0; + for (int64_t i = 0; i < n - 1; ++i) { + if (abs(H(i,i)) > max_coeff) { + max_coeff = abs(H(i,i)); + } + } + norm_bound = 1/max_coeff; + //std::cout << "H = \n" << H << "\n"; + //std::cout << "A = \n" << A << "\n"; + //std::cout << "B = \n" << B << "\n"; + + //std::cout << "A*B = \n" << A*B << "\n"; + //std::cout << "y = \n" << y << "\n"; + // I've never observed this to happen; + // but I also haven't been able to find a proof that the norm is nondecreasing. + if (norm_bound < last_norm_bound) { + std::cerr << "Norm bound has decreased!\n"; + } + last_norm_bound = norm_bound; + ++iteration; + if (iteration % 250 == 0) { + prog.display_progress(iteration, expected_iterations, static_cast(norm_bound)); + } + } + return relation; +} + + +template +std::string pslq(std::map const & dictionary, Real max_acceptable_norm_bound, Real gamma, std::ostream& os = std::cout) { + using std::abs; + std::vector values(dictionary.size()); + size_t i = 0; + for (auto it = dictionary.begin(); it != dictionary.end(); ++it) { + values[i++] = it->first; + } + + auto m = pslq(values, max_acceptable_norm_bound, gamma, os); + if (m.size() > 0) { + std::ostringstream oss; + auto const & symbol = dictionary.find(m[0].second)->second; + oss << "As\n\t"; + Real sum = m[0].first*m[0].second; + oss << m[0].first << "⋅" << m[0].second; + for (size_t i = 1; i < m.size(); ++i) + { + if (m[i].first < 0) { + oss << " - "; + } else { + oss << " + "; + } + oss << abs(m[i].first) << "⋅" << m[i].second; + sum += m[i].first*m[i].second; + } + oss << " = " << sum << ",\nit is likely that\n\t"; + + oss << m[0].first << "⋅" << symbol; + for (size_t i = 1; i < m.size(); ++i) + { + if (m[i].first < 0) { + oss << " - "; + } else { + oss << " + "; + } + auto const & symbol = dictionary.find(m[i].second)->second; + oss << abs(m[i].first) << "⋅" << symbol; + } + oss << " = 0."; + + return oss.str(); + } + return ""; +} + +template +std::string pslq(std::map const & dictionary, Real max_acceptable_norm, std::ostream& os = std::cout) { + using std::sqrt; + Real gamma = 2/sqrt(3) + 0.01; + return pslq(dictionary, max_acceptable_norm, gamma, std::cout); +} + +template +std::string identify(Real x, std::string symbol, Real max_acceptable_norm, std::ostream& os = std::cout) { + auto dictionary = standard_pslq_dictionary(); + using std::sqrt; + Real gamma = 2/sqrt(3) + 0.01; + + using std::log; + using std::exp; + dictionary.emplace(x, symbol); + Real log_ = log(x); + if (log_ < 0) { + dictionary.emplace(-log(x), "-ln(" + symbol + ")"); + } else { + dictionary.emplace(log(x), "ln(" + symbol + ")"); + } + dictionary.emplace(exp(x), "exp(" + symbol + ")"); + dictionary.emplace(1/x, "1/" + symbol); + dictionary.emplace(x*x, symbol + "²"); + return pslq(dictionary, max_acceptable_norm, gamma, os); +} + +template +std::string is_algebraic(Real alpha, std::string symbol, Real max_acceptable_norm_bound) { + // TODO: Figure out this interface. + std::vector x(80, Real(1)); + for (size_t i = 1; i < x.size(); ++i) { + x[i] = x[i-1]*alpha; + } + using std::sqrt; + Real gamma = 2/sqrt(Real(3)) + 0.01; + std::vector> sol = pslq(x, max_acceptable_norm_bound, gamma); + if (sol.size() > 0) { + std::cout << "Solution has elements "; + for (auto [mi, xi] : sol) { + std::cout << mi << ", "; + } + return "Found a solution"; + } + return ""; +} + +} +#endif diff --git a/performance/pslq_performance.cpp b/performance/pslq_performance.cpp new file mode 100644 index 000000000..0c12b7ec0 --- /dev/null +++ b/performance/pslq_performance.cpp @@ -0,0 +1,69 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include +#include +#include +#include +#if __has_include() +#include +#endif +template +inline Real rsqrt(Real const & x) { + if constexpr (std::is_same_v) { + __m128 b, t; + b = _mm_set_ss (x); + t = _mm_rsqrt_ss (b); + float res; + _mm_store_ss (&res, t); + //res = res*(3-x*res*res)/2; + res = res + res*(1-x*res*res)/2; + return res; + } + if constexpr (std::is_same_v || std::is_same_v) { + double a = rsqrt(float(x)); + a = a + a*(1-x*a*a)/2; + return a; + //return 1/sqrt(x); + } + if constexpr (std::is_same_v) { + boost::multiprecision::float128 a = rsqrt(static_cast(x)); + // one newton iteration: + a = a*(3-x*a*a)/2; + // a second newton iterate: + return a*(3-x*a*a)/2; + } + else if constexpr (sizeof(Real) > 16) { + Real res; + mpfr_rec_sqrt (res.backend().data(), x.backend().data(), GMP_RNDN); + return res; + } +} + +template +void RSqrtBM(benchmark::State& state) +{ + Real x = 0.01; + for (auto _ : state) { + benchmark::DoNotOptimize(rsqrt(x)); + x += std::numeric_limits::epsilon(); + } + + state.SetComplexityN(8*sizeof(Real)); +} + +BENCHMARK_TEMPLATE(RSqrtBM, float)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, double)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, long double)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::float128)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::number>)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::number>)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::number>)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::number>)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::number>)->Complexity(); +BENCHMARK_TEMPLATE(RSqrtBM, boost::multiprecision::number>)->Complexity(); + +BENCHMARK_MAIN(); \ No newline at end of file diff --git a/test/test_pslq.cpp b/test/test_pslq.cpp new file mode 100644 index 000000000..0ba7398d8 --- /dev/null +++ b/test/test_pslq.cpp @@ -0,0 +1,70 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include "../../math/test/math_unit_test.hpp" +#include +#include +#include + +using std::pow; +using boost::multiprecision::is_algebraic; + +// A test: A := 2^1/6 + 3^1/5 +// From: Applications of integer relation algorithms, Jonathan M. Borwein, Petr Lisonek, Discrete Mathematics 217 (2000) 65{82 + + +// Getting algebraic numbers: https://www.davidhbailey.com/dhbpapers/pslq-cse.pdf +// B3=3.54409035955···[1], then B3 is the root of 4913 + 2108t^2 − 604t^3 − 977t^4 + 8t^5 + 44t^6 + 392t^7 − 193t^8 − 40t^9 + 48t^10 − 12t^11 + t^12 + + +// This is test 1 from: +// https://www.davidhbailey.com/dhbpapers/mpfun2015.pdf +template +void test_degree_30_algebraic() +{ + Real x = pow(3, Real(1)/Real(5)) - pow(2, Real(1)/Real(6)); + // DHB's code requires 240 decimal digits. + std::pair p{x, "{3^1/5 - 2^1/6}"}; + auto b = boost::multiprecision::is_algebraic(p, Real(1e10)); + if (!b.empty()) { + std::cout << "It is algebraic.\n"; + } +} + +// This is test 3 from: +// https://www.davidhbailey.com/dhbpapers/mpfun2015.pdf +template +void test_degree_56_algebraic() +{ + Real x = pow(3, Real(1)/Real(7)) - pow(2, Real(1)/Real(8)); + // DHB's code requires 640 decimal digits. + std::pair p{x, "{3^1/7 - 2^1/8}"}; + auto b = boost::multiprecision::is_algebraic(p, Real(1e10)); + if (!b.empty()) { + std::cout << "It is algebraic.\n"; + } +} + +template +void test_degree_72_algebraic() +{ + Real x = pow(3, Real(1)/Real(8)) - pow(2, Real(1)/Real(9)); + // DHB's code requires 1100 decimal digits. + std::pair p{x, "{3^1/8 - 2^1/9}"}; + auto b = boost::multiprecision::is_algebraic(p, Real(1e10)); + if (!b.empty()) { + std::cout << "It is algebraic.\n"; + } +} + + +int main() { + using Real = boost::multiprecision::number >; + test_degree_30_algebraic(); + test_degree_56_algebraic(); + test_degree_72_algebraic(); + return boost::math::test::report_errors(); +} \ No newline at end of file From a3508bdb5514d850eb4b0b929776d01660c6842a Mon Sep 17 00:00:00 2001 From: Nick Date: Fri, 24 Jul 2020 11:14:18 -0400 Subject: [PATCH 02/12] Add Fransen-Robinson constant to PSLQ examples. [CI SKIP] --- example/representations/fransen_robinson.cpp | 47 ++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 example/representations/fransen_robinson.cpp diff --git a/example/representations/fransen_robinson.cpp b/example/representations/fransen_robinson.cpp new file mode 100644 index 000000000..ccf4591ca --- /dev/null +++ b/example/representations/fransen_robinson.cpp @@ -0,0 +1,47 @@ +/* + * Copyright Nick Thompson 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ +#include "analyze.hpp" +#include +#include +#include + +using boost::multiprecision::mpfr_float; +using boost::math::quadrature::exp_sinh; +using std::exp; +using std::log; + +template +Real fransen_robinson() +{ + auto integrator = exp_sinh(18); + auto f = [](Real t)->Real { + Real y; + try { + y = boost::math::tgamma(t); + } catch (std::exception const & e) { + std::cout << "Problem at t= " << t << ":" << e.what() << "\n"; + return Real(0); + } + + return 1/y; + }; + Real Q = integrator.integrate(f); + return Q; +} + +int main() +{ + using Real = mpfr_float; + int p = 500; + mpfr_float::default_precision(p); + std::cout << std::setprecision(p); + Real F = fransen_robinson(); + std::cout << "Expected = 2.8077702420285\n"; + std::cout << "Computed = " < Date: Fri, 24 Jul 2020 11:47:26 -0400 Subject: [PATCH 03/12] Reciprocal Fibonacci constant. [CI SKIP] --- .../representations/reciprocal_fibonacci.cpp | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 example/representations/reciprocal_fibonacci.cpp diff --git a/example/representations/reciprocal_fibonacci.cpp b/example/representations/reciprocal_fibonacci.cpp new file mode 100644 index 000000000..9005528f3 --- /dev/null +++ b/example/representations/reciprocal_fibonacci.cpp @@ -0,0 +1,28 @@ +#include +#include "analyze.hpp" + +template +Real reciprocal_fibonacci() +{ + Real x0 = 1; + Real x1 = 1; + Real sum = 2; + Real diff = 1; + while (diff > std::numeric_limits::epsilon()) { + Real tmp = x1 + x0; + diff = 1/tmp; + sum += diff; + x0 = x1; + x1 = tmp; + } + return sum; +} + +int main() { + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(400); + using Real = mpfr_float; + Real rf = reciprocal_fibonacci(); + std::cout << "r = 3.35988566624317755317201130291892717968890513373196848649555381532513031899668338361541621645679008729704\n"; + analyze(rf, "φ"); +} From 5e646d0d7f8f9afea213974360114d014721def1 Mon Sep 17 00:00:00 2001 From: Nick Date: Thu, 30 Jul 2020 11:50:58 -0400 Subject: [PATCH 04/12] Remove Dottie and Laplace limit from examples, as they are now in boost::math::constants. [CI SKIP] --- doc/reference_pslq.qbk | 12 ++++-- example/representations/constants.cpp | 4 +- example/representations/dottie.cpp | 25 ------------ example/representations/laplace_limit.cpp | 34 ----------------- .../representations/reciprocal_fibonacci.cpp | 28 -------------- include/boost/multiprecision/pslq.hpp | 38 +++++++++++++++++++ 6 files changed, 49 insertions(+), 92 deletions(-) delete mode 100644 example/representations/dottie.cpp delete mode 100644 example/representations/laplace_limit.cpp delete mode 100644 example/representations/reciprocal_fibonacci.cpp diff --git a/doc/reference_pslq.qbk b/doc/reference_pslq.qbk index 8711e14eb..1286f58c7 100644 --- a/doc/reference_pslq.qbk +++ b/doc/reference_pslq.qbk @@ -59,8 +59,9 @@ Even though PSLQ does not produce proofs, pessimism is unwarranted: PSLQ is a to In addition, /verifying/ a suspected relation is much less computationally difficult than discovering it, so once a relation is found, we can compute the relevant terms of the dictionary to much higher precision to give additional evidence for the relation, and then search for a formal proof. Astute readers will notice that our example requires a very auspicious guess about the form the relation takes. -If we didn't know that β was somehow related to π⁶ and ζ(3)², we would've never found this relation. -This is a general property of the PSLQ algorithm: As its runtime scales quartically in the length of the dictionary, having a hunch about the form the relation will take is necessary. +If we didn't know that β was somehow related to π⁶ and ζ(3)², we would've never found it. +This is a general property of the PSLQ algorithm: +As its runtime scales quartically in the length of the dictionary, having a hunch about the form the relation will take is necessary. That said, simple modifications of the inputs allow us to recover different classes of relations. Let's use the Basel problem as an example: @@ -84,7 +85,7 @@ and if it is indeed algebraic with reasonable sized polynomial coefficients, the Communicating this idea has a somewhat different format than vanilla `pslq`, and hence we have provided a wrapper `is_algebraic` for precisely this purpose. Quartically scaling algorithms are rare, and users generally have little experience with them. -It is easy to construct a PSLQ dictionary that would only complete after a few thousand years. +It is easy to construct a PSLQ dictionary that induces a computation lasting a few thousand years. Hence we have defaulted to printing a progress bar to the terminal. @@ -105,6 +106,9 @@ We can silence the progress bar via by passing a `std::ostream` with the badbit ofs.setstate(std::ios_base::badbit); std::string s = boost::multiprecision::pslq(m, max_acceptable_norm_bound, ofs); +Communication of the result is done via strings both decimal strings and unicode characters. +We have attempted to maintain standard notation for the constants provided in the standard PSLQ dictionary, but if you can't determine the notation, all the constants can be looked up in the Online Encyclopedia of Integer Sequences. -[endsect] \ No newline at end of file + +[endsect] diff --git a/example/representations/constants.cpp b/example/representations/constants.cpp index be03c5b93..c7c6998f1 100644 --- a/example/representations/constants.cpp +++ b/example/representations/constants.cpp @@ -23,5 +23,7 @@ int main() { analyze(sin_one(), "sin(1)"); analyze(cos_one(), "cos(1)"); analyze(e_pow_pi(), "e^pi"); - + analyze(dottie(), "d"); + analyze(laplace_limit(), "λ"); + analyze(reciprocal_fibonacci(), "ψ"); } diff --git a/example/representations/dottie.cpp b/example/representations/dottie.cpp deleted file mode 100644 index 358b22f13..000000000 --- a/example/representations/dottie.cpp +++ /dev/null @@ -1,25 +0,0 @@ -#include -#include "analyze.hpp" - -template -Real dottie() -{ - using std::cos; - using std::abs; - using std::sin; - Real x{".739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246"}; - Real residual = cos(x) - x; - do { - x += residual/(sin(x)+1); - residual = cos(x) - x; - } while(abs(residual) > std::numeric_limits::epsilon()); - return x; -} - -int main() { - using boost::multiprecision::mpfr_float; - mpfr_float::default_precision(1300); - using Real = mpfr_float; - Real dot = dottie(); - analyze(dot, "dottie"); -} diff --git a/example/representations/laplace_limit.cpp b/example/representations/laplace_limit.cpp deleted file mode 100644 index 5e820eb7d..000000000 --- a/example/representations/laplace_limit.cpp +++ /dev/null @@ -1,34 +0,0 @@ -#include -#include "analyze.hpp" - -template -Real laplace_limit() -{ - using std::cos; - using std::abs; - using std::sin; - using std::sqrt; - using std::exp; - Real x{"0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168"}; - Real tmp = sqrt(1+x*x); - Real etmp = exp(tmp); - Real residual = x*exp(tmp) - 1 - tmp; - Real df = etmp -x/tmp + etmp*x*x/tmp; - do { - x -= residual/df; - tmp = sqrt(1+x*x); - etmp = exp(tmp); - residual = x*exp(tmp) - 1 - tmp; - df = etmp -x/tmp + etmp*x*x/tmp; - } while(abs(residual) > std::numeric_limits::epsilon()); - std::cout << "Residual = " << residual << "\n"; - return x; -} - -int main() { - using boost::multiprecision::mpfr_float; - mpfr_float::default_precision(1300); - using Real = mpfr_float; - Real ll = laplace_limit(); - analyze(ll, "laplace_limit"); -} diff --git a/example/representations/reciprocal_fibonacci.cpp b/example/representations/reciprocal_fibonacci.cpp deleted file mode 100644 index 9005528f3..000000000 --- a/example/representations/reciprocal_fibonacci.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include "analyze.hpp" - -template -Real reciprocal_fibonacci() -{ - Real x0 = 1; - Real x1 = 1; - Real sum = 2; - Real diff = 1; - while (diff > std::numeric_limits::epsilon()) { - Real tmp = x1 + x0; - diff = 1/tmp; - sum += diff; - x0 = x1; - x1 = tmp; - } - return sum; -} - -int main() { - using boost::multiprecision::mpfr_float; - mpfr_float::default_precision(400); - using Real = mpfr_float; - Real rf = reciprocal_fibonacci(); - std::cout << "r = 3.35988566624317755317201130291892717968890513373196848649555381532513031899668338361541621645679008729704\n"; - analyze(rf, "φ"); -} diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp index 516032261..2013137e6 100644 --- a/include/boost/multiprecision/pslq.hpp +++ b/include/boost/multiprecision/pslq.hpp @@ -166,6 +166,7 @@ auto standard_pslq_dictionary() { m.emplace(1/glaisher_, "1/A"); m.emplace(log(glaisher_), "ln(A)"); m.emplace(exp(glaisher_), "exp(A)"); + Real khinchin_ = khinchin(); m.emplace(khinchin_, "K₀"); m.emplace(log(khinchin_), "ln(K₀)"); @@ -202,6 +203,43 @@ auto standard_pslq_dictionary() { m.emplace(root_three()*pi_, "π√3"); m.emplace(catalan_*pi_, "πG"); + // Lemniscate constant: + auto lemniscate = pi_*gauss(); + m.emplace(lemniscate, "L"); + m.emplace(log(lemniscate), "ln(L)"); + m.emplace(exp(lemniscate), "exp(L)"); + m.emplace(lemniscate*lemniscate, "L²"); + m.emplace(1/lemniscate, "1/L"); + + // Plastic constant: + auto P = plastic(); + m.emplace(P, "P"); + m.emplace(log(P), "ln(P)"); + m.emplace(exp(P), "exp(P)"); + m.emplace(P*P, "P²"); + m.emplace(1/P, "1/P"); + + // Dottie number: x = cos(x). + auto dot_ = dottie(); + m.emplace(dot_, "d"); + m.emplace(-log(dot_), "-ln(d)"); + m.emplace(exp(dot_), "exp(d)"); + m.emplace(dot_*dot_, "d²"); + m.emplace(1/dot_, "1/d"); + + auto rfib = reciprocal_fibonacci(); + m.emplace(rfib, "ψ"); + m.emplace(log(rfib), "ln(ψ)"); + m.emplace(rfib*rfib, "ψ²"); + m.emplace(exp(rfib), "exp(ψ)"); + m.emplace(1/rfib, "1/ψ"); + + auto ll = laplace_limit(); + m.emplace(ll, "λ"); + m.emplace(ll*ll, "λ²"); + m.emplace(1/ll, "1/λ"); + m.emplace(-log(ll), "-ln(λ)"); + m.emplace(exp(ll), "exp(λ)"); return m; } From 6d2117ae4b5a58f683adcf3000dd9e025df505ad Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 2 Aug 2020 20:25:42 -0400 Subject: [PATCH 05/12] Analyze Soldner's constant [CI SKIP] --- example/representations/soldner.cpp | 38 +++++++++++++++++++++++++++ include/boost/multiprecision/pslq.hpp | 24 ++++++++++++++--- 2 files changed, 58 insertions(+), 4 deletions(-) create mode 100644 example/representations/soldner.cpp diff --git a/example/representations/soldner.cpp b/example/representations/soldner.cpp new file mode 100644 index 000000000..af76f79da --- /dev/null +++ b/example/representations/soldner.cpp @@ -0,0 +1,38 @@ +#include +#include +#include "analyze.hpp" + +// See: https://en.wikipedia.org/wiki/Ramanujan%E2%80%93Soldner_constant +template +T soldner() +{ + using boost::math::expint; + using std::log; + using std::abs; + // Initial guess from https://oeis.org/A070769/constant + T x{"1.45136923488338105028396848589202744949303228364801586309300455766242559575451783565953135771108682884"}; + T y = expint(log(x)); + // If x is the root, then + // li(x(1+μ)) = xμli'(x) + Ο(μ²) ≈ xμ/log(x), + // where μ = ε/2 is the unit roundoff, and we know x ≈ 1.4513. + // Plugging this in gets: + const T expected_residual = 1.9481078669535834*std::numeric_limits::epsilon(); + // This should that we cannot expect better than ~2ULPs of accuracy using Newton's method. + // Just a couple extra bits of accuracy should get this recover full precision, + // but variable precision always causes us pain . . . + while (abs(y) > expected_residual) { + // Use li(x) = Ei(ln(x)). + T ln = log(x); + y = expint(ln); + x -= ln*y; + } + return x; +} + +int main() { + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(400); + using Real = mpfr_float; + Real s = soldner(); + analyze(s, "μ"); +} diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp index 2013137e6..b9f0587a4 100644 --- a/include/boost/multiprecision/pslq.hpp +++ b/include/boost/multiprecision/pslq.hpp @@ -43,8 +43,9 @@ #include #include #include -#include #include +#include // for Ei(1), Ei(-1). +#include #include #if defined __has_include # if __has_include () @@ -127,7 +128,7 @@ auto standard_pslq_dictionary() { m.emplace(zeta_three_*zeta_three_, "ζ(3)²"); m.emplace(zeta_three_*zeta_three_*zeta_three_, "ζ(3)³"); m.emplace(pow(zeta_three_, 4), "ζ(3)⁴"); - + auto pi_ = pi(); m.emplace(pi_, "π"); m.emplace(1/pi_, "1/π"); @@ -175,7 +176,7 @@ auto standard_pslq_dictionary() { m.emplace(khinchin_*khinchin_, "K₀²"); m.emplace(log(static_cast(1) + sqrt(static_cast(2))), "ln(1+√2)"); // To recover multiplicative relations we need the logarithms of small primes. - Real ln2 = log(static_cast(2)); + Real ln2 = log(static_cast(2)); m.emplace(ln2, "ln(2)"); m.emplace(-log(ln2), "-ln(ln(2))"); m.emplace(log(static_cast(3)), "ln(3)"); @@ -240,6 +241,21 @@ auto standard_pslq_dictionary() { m.emplace(1/ll, "1/λ"); m.emplace(-log(ll), "-ln(λ)"); m.emplace(exp(ll), "exp(λ)"); + + auto ei1 = boost::math::expint(Real(1)); + m.emplace(ei1, "Ei(1)"); + m.emplace(ei1*ei1, "Ei(1)²"); + m.emplace(1/ei1, "1/Ei(1)"); + m.emplace(log(ei1), "ln(Ei(1))"); + m.emplace(exp(ei1), "exp(Ei(1))"); + + auto eim1 = boost::math::expint(-Real(1)); + m.emplace(-eim1, "-Ei(-1)"); + m.emplace(eim1*eim1, "Ei(-1)²"); + m.emplace(-1/eim1, "-1/Ei(-1)"); + m.emplace(-log(-eim1), "-ln(-Ei(-1))"); + m.emplace(exp(eim1), "exp(Ei(-1))"); + return m; } @@ -366,7 +382,7 @@ std::vector> pslq(std::vector & x, Real max_accep for (int64_t i = n - 2; i >= 0; --i) { s_sq[i] = s_sq[i+1] + x[i]*x[i]; } - + using std::pow; if (max_acceptable_norm_bound*max_acceptable_norm_bound*s_sq[0] > 1/std::numeric_limits::epsilon()) { std::cerr << "The maximum acceptable norm bound " << max_acceptable_norm_bound << " is too large; spurious relations will be recovered.\n"; From a76e8637b11e44ba1e5a3c7c6bc49c11f2cbc6fa Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Tue, 4 Aug 2020 14:16:51 -0400 Subject: [PATCH 06/12] Update documentation [CI SKIP] --- doc/reference_pslq.qbk | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/doc/reference_pslq.qbk b/doc/reference_pslq.qbk index 1286f58c7..5c551a26f 100644 --- a/doc/reference_pslq.qbk +++ b/doc/reference_pslq.qbk @@ -109,6 +109,26 @@ We can silence the progress bar via by passing a `std::ostream` with the badbit Communication of the result is done via strings both decimal strings and unicode characters. We have attempted to maintain standard notation for the constants provided in the standard PSLQ dictionary, but if you can't determine the notation, all the constants can be looked up in the Online Encyclopedia of Integer Sequences. +As an example, let's analyze the extreme value skewness: +The C++ code is + s = identify(extreme_value_skewness(), "extreme_value_skewness", Real(1e10)); + if (!s.empty()) + { + std::cout << s << "\n"; + } + else + { + std::cout << "No relations found.\n"; + } + +and the output is + + As + -2⋅0.130631 + 2⋅0.184034 + 5⋅0.693147 + 3⋅1.09861 - 6⋅1.14473 = 1.90276e-1000, + it is likely that + -2⋅ln(extreme_value_skewness) + 2⋅ln(ζ(3)) + 5⋅ln(2) + 3⋅ln(3) - 6⋅ln(π) = 0. + +We can collect the logarithms to find the correct relation. [endsect] From 52af1e89ccf5eb31be5814792d339e7a9aad964e Mon Sep 17 00:00:00 2001 From: Nick Date: Tue, 4 Aug 2020 15:02:01 -0400 Subject: [PATCH 07/12] Add Dirichlet L-function values. [CI SKIP] --- example/representations/analyze.hpp | 16 +-- example/representations/dirichlet_l.cpp | 156 ++++++++++++++++++++++++ 2 files changed, 164 insertions(+), 8 deletions(-) create mode 100644 example/representations/dirichlet_l.cpp diff --git a/example/representations/analyze.hpp b/example/representations/analyze.hpp index 44c0f564e..e205da97a 100644 --- a/example/representations/analyze.hpp +++ b/example/representations/analyze.hpp @@ -37,13 +37,13 @@ void analyze(Real x, std::string symbol) std::cout << symbol << " ≈ " << x << "\n"; auto scf = simple_continued_fraction(x); - std::cout << " ≈ " << scf << "\n"; + std::cout << symbol << " ≈ " << scf << "\n"; auto ccf = centered_continued_fraction(x); - std::cout << " ≈ " << ccf << "\n"; + std::cout << symbol << " ≈ " << ccf << "\n"; auto lur = luroth_expansion(x); - std::cout << " ≈ " << lur << "\n"; + std::cout << symbol << " ≈ " << lur << "\n"; auto eng = engel_expansion(x); - std::cout << " ≈ " << eng << "\n"; + std::cout << symbol << " ≈ " << eng << "\n"; std::cout << std::setprecision(11); std::cout << "Empirical Khinchin geometric mean of simple cfrac : " << scf.khinchin_geometric_mean() << "\n"; std::cout << "Expected Khinchin geometric mean if the value is 'random': 2.6854520010.\n"; @@ -61,8 +61,8 @@ void analyze(Real x, std::string symbol) auto lur_exp = luroth_expansion(expx); std::cout << "exp(" << symbol << ")"; std::cout << " ≈ " << scf_exp << "\n"; - std::cout << " ≈ " << ccf_exp << "\n"; - std::cout << " ≈ " << lur_exp << "\n"; + std::cout << "exp(" << symbol << ") ≈ " << ccf_exp << "\n"; + std::cout << "exp(" << symbol << ") ≈ " << lur_exp << "\n"; std::cout << std::setprecision(11); std::cout << "Empirical Khinchin geometric mean of simple cfrac : " << scf_exp.khinchin_geometric_mean() << "\n"; std::cout << "Expected Khinchin geometric mean if the value is 'random': 2.6854520010.\n"; @@ -79,8 +79,8 @@ void analyze(Real x, std::string symbol) auto lur_log = luroth_expansion(lnx); std::cout << "ln(" << symbol << ")"; std::cout << " ≈ " << scf_log << "\n"; - std::cout << " ≈ " << ccf_log << "\n"; - std::cout << " ≈ " << lur_log << "\n"; + std::cout << "ln(" << symbol << ") ≈ " << ccf_log << "\n"; + std::cout << "ln(" << symbol << ") ≈ " << lur_log << "\n"; std::cout << std::setprecision(11); std::cout << "Empirical Khinchin geometric mean of simple cfrac : " << scf_log.khinchin_geometric_mean() << "\n"; std::cout << "Expected Khinchin geometric mean if the value is 'random': 2.6854520010.\n"; diff --git a/example/representations/dirichlet_l.cpp b/example/representations/dirichlet_l.cpp new file mode 100644 index 000000000..8e466b0c9 --- /dev/null +++ b/example/representations/dirichlet_l.cpp @@ -0,0 +1,156 @@ +#include + +#include "analyze.hpp" + +int main() +{ + using boost::multiprecision::mpfr_float; + mpfr_float::default_precision(700); + using Real = mpfr_float; + + // See: Mathematical Constants II, Section 1.13 by Steven Finch. + // Mathematica code: N[DirichletL[5, 3, 3], 1000] + Real L53 { + "0.85482476664854301023569008353813769713839646493700528273070249938112\ +3833412689428128420956718709787428514052771595718349505229519170836211\ +0265215997608202788987421006122517302151023228314451262775716209021353\ +0330103464538993007686744724563146241599644743502492906957776462953182\ +9667841632135988402566499417190069150908562376936330929748257040245694\ +0947636595746622500581557129589490454239050805441470697914692697461925\ +5354856511779539292609535748785502504575456157533565495931814259664978\ +5624266800940770031423101453452565610482192953087267421460174314030751\ +3137341326100712002598205983618871129855306286794197785385570173753286\ +7373383080418131055820547604575044484604288354333927233062704861657761\ +1775573273825812512656536391973058944408433714228182173403526502740377\ +8019790317047796978242086364457235552140447881646673112125614034723521\ +9461836783049496989354309630355462809269796478766558651443420889056879\ +0653313860886540271978377580376324147334672995755120702038915803888794\ +683762686965986515690660620" + }; + analyze(L53, "L₅(3)"); + + // L_{12}(3) in Finch. + // Mathematica: N[DirichletL[12, 4, 3], 1005] + Real L12_3 { + "0.99004001943815994979181677686330405085688506765723614555366070034235\ +2053367181167785602231853515063548035652066831323529149037644412048770\ +1205677113678637008494970351453363227396656628133639036437820155148361\ +1523912758832677415308095540032679791548733465552729879160705970944068\ +3047247231869026571681015297783981794857372144450793139427295178127033\ +0454830070083925459978033457642066581717231626579344079820901719429178\ +5471870477025279791537416710415876608877730795883237914004459636823272\ +7702962909997958054571449426461001201274802599416886077744847672774823\ +1616819800264336530310270114499097578897361625552822358991324229683565\ +2440419849829503418120486666108009658007449371269272319830311342989368\ +9256988033289311190486011752634209285757168374562341354205902592319993\ +6753158050276247003042601732659629733370449380106911047741381993467677\ +9120957028579431667039129344224731359897428088131915928424606912198955\ +5791359848206367934669687971937601745357881392613210087688721247408254\ +223667773287141511792169284" + }; + + analyze(L12_3, "L₁₂(3)"); + + // L_8(3) in Finch. + // Mathematica: N[DirichletL[8, 2, 3], 1005] + Real L8_3 { + "0.95838045456309456205166940286157781882489531793977534075750450704707\ +5697484297936478252699777513729824780207090773293549634062455175352154\ +1050421643741242162948370712212459530531743124610365262262393362356340\ +7242366444356379467853814079517704296723654750177624408938046768924642\ +1346482676922365570029978820163032037998590615390431796035219020488577\ +9449673184860010370260657180847741482324488852267107290266951818939892\ +1621831906571565400608971913016646314300278934890667662255118229921345\ +1913684279204871333343241172881130708166426611155795048684486613113374\ +6476116522838928371408593020450368620451306094378769166720567324494314\ +2674178390819968326759193574012172232989305911322079647001333138576073\ +5927156685235137869515986545912527991281085063341437711573616218472682\ +4511230570005373326901436363158608189699371163775422771423459441606967\ +0456346938840266465193980885080436207621407981512000693167558302132361\ +2680561330362776119652142434648067880966579954470979005234252562531647\ +121115184737492615995093421" + }; + analyze(L8_3, "L₈(3)"); + + // L_1(1/2) in Finch. + // Mathematica: N[DirichletL[1, 1, 1/2], 1000] + Real mL1_12{ + "1.4603545088095868128894991525152980124672293310125814905428860878255\ +3052947450062527641937546335681951449637467986952958389234371035889426\ +1819232839753762925182633358649164127891229394154101197917310448108241\ +9409278816984288571768239557991845178836146554866593799168915231635216\ +0424275374940796571353042261006512411854164516150879658464017534002986\ +1665046125952490371908917840714940471565917574660272483532061693045430\ +4928626602599081880472531306067400305323926604512813065106539466001040\ +6651804549168506332768620827425551694390772400674037153616973833187775\ +8235704458653739288227657621030550745336237111852621087115124502607502\ +1266249685083684103157725190613240558606006778558775909756029975080223\ +1586212263122265454248392384908383946490628221225677591832636143132813\ +2287132579942988483462556863862335628109251185853308139377517130573507\ +6952743265652282989458373709112836747583785163266961089292683262398366\ +1271664820364213483392014697899141766404504399312667163511583342287816\ +228024986419668541471350141" + }; + analyze(mL1_12, "-L₁(1/2)"); + + // L_5(1/2) in Finch: + // Mathematica: N[DirichletL[5, 3, 1/2], 1005] + Real L5_12{ +"0.23175094750401575588338366176087722642788669640900589796611788733984\ +4380645423102142258987289424254174755735228412712566384062820509308705\ +4730991946118103108581508479940151667972848978549079821793104180029987\ +5808507633045155604499296665559640047056418881504273930695195039987374\ +9655387134942636424847345300222199526500963040449084966940305170478514\ +4300079352976507407869411645777801482430353783545520772060016388308645\ +6498069592935270987227131901916252715681404294521557860420838536190356\ +9730344478230686837854828680907742467760509871791222403042521404654094\ +6446559685196147827251082022649629787298601994056593163378178215867996\ +4027875252373167373305899703159404849664172429970895821195499062461733\ +9817714751567388988924686442324629989480323417740562232558232682450499\ +2053602641637595626222928012574680834901035016208006933977266581926435\ +0764872720074152097782617729596469368980369023619589084185543173668370\ +3533159284204843689808952892140190000233076861509907039016592954272464\ +521698841993288727173249743" + }; + + analyze(L5_12, "L₅(1/2"); + + // L_8(1/2) in Finch. + // Mathematica: N[DirichletL[8, 2, 1/2], 1005] + Real L8_12{"0.37369171291254730738158695002302311358644030221942900028304279605268\ +8183524339117594442803625671304702131077889794637378068813783175365256\ +4228969235084455400241311061776440622457975356457040743657088834425181\ +7357991137259848085703678164473027253401522847518587289315109434094171\ +7444363261741602577602717616861910421573365912961709337399657911462957\ +3639450771314130299904415250314469370672907171745570365727794190109146\ +2280149311202718143240016030773469890979198759091698822844939533674606\ +5679802580360150255051367798332198854668044180785360951193826398082141\ +8112635335305350497631910663227288800613882791338839931442299198287344\ +5159976742277467250474852234990657825867106235246197561277554211855248\ +3591978733700958283435834248424248779317875050506914790119476720910613\ +6511229200085886977434958519320673781798168806336012820273998971315408\ +4643183748615646439178661278652300905647925231942445919666472949195541\ +1380255030327201982860866296098884131318085230680147320324577294354703\ +029741887138677182829463823"}; + + analyze(L8_12, "L₈(1/2)"); + + // L_12(1/2) in Finch. + // Mathematica: N[DirichletL[12, 4, 1/2], 1005] + Real L12_12{"0.49855700245781543615675655398717998948832524973455721249708152216115\ +7514707901812918935494759798810333046573988572719765227147980916611124\ +7107730884702644877413152667518384205725484683504403135875688049833927\ +8082929333181342711772854005849616100527423608525207702010076677005495\ +2890704221336285435994536795081158255306386457761921529268304145568147\ +3567084558172992220716198750906923718587467434855320798360613623299817\ +0151673117280670418436275977492286419761757520829226040561543752207756\ +1413596377205660526731601068621787084655962694840091994230126051019486\ +2117946997729561254946124027945896584186401108861403515339247925948516\ +4794596728588038318241537842864085928506482284962410734766005977222421\ +3935905197490341443351734744472884358075773268497774668595703863116915\ +6326589879274823186502460462325430321601195107062271124593229787394103\ +8922259725038217623474768058873466381354423737295237504433036123093405\ +6016204032005908651267255460196545055772702201044857358408719970403031\ +562539200865086312925653630"}; + analyze(L12_12, "L₁₂(1/2)"); +} From f21aa5385d77ac047aa09b322c07a0f61a1227bd Mon Sep 17 00:00:00 2001 From: Nick Date: Tue, 4 Aug 2020 15:39:30 -0400 Subject: [PATCH 08/12] Add gamma(1/2), gamma(1/3). [CI SKIP] --- include/boost/multiprecision/pslq.hpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp index b9f0587a4..011de385b 100644 --- a/include/boost/multiprecision/pslq.hpp +++ b/include/boost/multiprecision/pslq.hpp @@ -45,6 +45,7 @@ #include #include #include // for Ei(1), Ei(-1). +#include // For Γ(1/3), Γ(1/4) #include #include #if defined __has_include @@ -256,6 +257,17 @@ auto standard_pslq_dictionary() { m.emplace(-log(-eim1), "-ln(-Ei(-1))"); m.emplace(exp(eim1), "exp(Ei(-1))"); + // These show up in a lot of identities in Finch: + auto gamma_14 = boost::math::tgamma(Real(1)/Real(4)); + m.emplace(gamma_14, "Γ(1/4)"); + m.emplace(log(gamma_14), "ln(Γ(1/4))"); + m.emplace(1/gamma_14, "1/Γ(1/4)"); + + auto gamma_13 = boost::math::tgamma(Real(1)/Real(3)); + m.emplace(gamma_13, "Γ(1/3)"); + m.emplace(log(gamma_13), "ln(Γ(1/3))"); + m.emplace(1/gamma_13, "1/Γ(1/3)"); + return m; } From 0e1bccbe4de0f83463b7bfffe2bfdf56edf3c5a9 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Tue, 4 Aug 2020 15:43:57 -0400 Subject: [PATCH 09/12] Slightly better verbiage [CI SKIP] --- example/representations/analyze.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/example/representations/analyze.hpp b/example/representations/analyze.hpp index e205da97a..b1671d1e0 100644 --- a/example/representations/analyze.hpp +++ b/example/representations/analyze.hpp @@ -91,7 +91,7 @@ void analyze(Real x, std::string symbol) std::cout << "Empirical digit geometric mean of Luroth expansion : " << lur_log.digit_geometric_mean() << "\n"; std::cout << "Expected digit geometric mean if the value is 'random' : 2.200161058\n"; - std::cout << "Is it algebraic?\n"; + std::cout << "Is " << symbol << " algebraic?\n"; std::string s = is_algebraic(x, symbol, Real(1e10)); if (!s.empty()) { @@ -102,7 +102,7 @@ void analyze(Real x, std::string symbol) std::cout << "It is probably not algebraic.\n"; } - std::cout << "Searching for closed forms.\n"; + std::cout << "Searching for closed forms of " << symbol << "\n"; s = identify(x, symbol, Real(1e10)); if (!s.empty()) { From c99b984446c560095a485ab4a3fc99f5da45941e Mon Sep 17 00:00:00 2001 From: NAThompson Date: Tue, 4 Aug 2020 15:57:06 -0400 Subject: [PATCH 10/12] gamma(1/4) cannot be included in the standard dictionary since its related to the lemniscate constant. [CI SKIP] --- include/boost/multiprecision/pslq.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp index 011de385b..484c47157 100644 --- a/include/boost/multiprecision/pslq.hpp +++ b/include/boost/multiprecision/pslq.hpp @@ -258,11 +258,6 @@ auto standard_pslq_dictionary() { m.emplace(exp(eim1), "exp(Ei(-1))"); // These show up in a lot of identities in Finch: - auto gamma_14 = boost::math::tgamma(Real(1)/Real(4)); - m.emplace(gamma_14, "Γ(1/4)"); - m.emplace(log(gamma_14), "ln(Γ(1/4))"); - m.emplace(1/gamma_14, "1/Γ(1/4)"); - auto gamma_13 = boost::math::tgamma(Real(1)/Real(3)); m.emplace(gamma_13, "Γ(1/3)"); m.emplace(log(gamma_13), "ln(Γ(1/3))"); From 7610a2fdd45a99c55c0ea06b75d951a30100b8e4 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Wed, 5 Aug 2020 12:25:55 -0400 Subject: [PATCH 11/12] Reintroduce gamma(1/4); but not ln(gamma(1/4)). [CI SKIP] --- example/representations/dirichlet_l.cpp | 4 ++-- example/representations/soldner.cpp | 2 +- include/boost/multiprecision/pslq.hpp | 6 ++++++ 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/example/representations/dirichlet_l.cpp b/example/representations/dirichlet_l.cpp index 8e466b0c9..14910fc10 100644 --- a/example/representations/dirichlet_l.cpp +++ b/example/representations/dirichlet_l.cpp @@ -5,7 +5,7 @@ int main() { using boost::multiprecision::mpfr_float; - mpfr_float::default_precision(700); + mpfr_float::default_precision(900); using Real = mpfr_float; // See: Mathematical Constants II, Section 1.13 by Steven Finch. @@ -113,7 +113,7 @@ int main() 521698841993288727173249743" }; - analyze(L5_12, "L₅(1/2"); + analyze(L5_12, "L₅(1/2)"); // L_8(1/2) in Finch. // Mathematica: N[DirichletL[8, 2, 1/2], 1005] diff --git a/example/representations/soldner.cpp b/example/representations/soldner.cpp index af76f79da..d212f8510 100644 --- a/example/representations/soldner.cpp +++ b/example/representations/soldner.cpp @@ -31,7 +31,7 @@ T soldner() int main() { using boost::multiprecision::mpfr_float; - mpfr_float::default_precision(400); + mpfr_float::default_precision(1000); using Real = mpfr_float; Real s = soldner(); analyze(s, "μ"); diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp index 484c47157..d040d70d9 100644 --- a/include/boost/multiprecision/pslq.hpp +++ b/include/boost/multiprecision/pslq.hpp @@ -258,6 +258,12 @@ auto standard_pslq_dictionary() { m.emplace(exp(eim1), "exp(Ei(-1))"); // These show up in a lot of identities in Finch: + auto gamma_14 = boost::math::tgamma(Real(1)/Real(3)); + m.emplace(gamma_14, "Γ(1/4)"); + // Don't put this in or the basis is linearly dependent with Gauss's constant: + //m.emplace(log(gamma_14), "ln(Γ(1/4))"); + m.emplace(1/gamma_14, "1/Γ(1/4)"); + auto gamma_13 = boost::math::tgamma(Real(1)/Real(3)); m.emplace(gamma_13, "Γ(1/3)"); m.emplace(log(gamma_13), "ln(Γ(1/3))"); From b89e233bbd5eefb0152718e0d44f86fad5d418f3 Mon Sep 17 00:00:00 2001 From: Nick Date: Fri, 7 Aug 2020 17:01:00 -0400 Subject: [PATCH 12/12] using declarations. [CI SKIP] --- include/boost/multiprecision/pslq.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/boost/multiprecision/pslq.hpp b/include/boost/multiprecision/pslq.hpp index d040d70d9..cbc838ac0 100644 --- a/include/boost/multiprecision/pslq.hpp +++ b/include/boost/multiprecision/pslq.hpp @@ -353,6 +353,11 @@ void display_progress(int64_t iteration, int64_t max_iterations, double norm_bou // See: https://www.davidhbailey.com/dhbpapers/cpslq.pdf, section 3. template std::vector> pslq(std::vector & x, Real max_acceptable_norm_bound, Real gamma, std::ostream& os = std::cout) { + using std::sqrt; + using std::log; + using std::ceil; + using std::pow; + using std::abs; std::vector> relation; /*if (!std::is_sorted(x.begin(), x.end())) { std::cerr << "Elements must be sorted in increasing order.\n"; @@ -360,7 +365,6 @@ std::vector> pslq(std::vector & x, Real max_accep }*/ std::sort(x.begin(), x.end()); - using std::sqrt; if (gamma <= 2/sqrt(3)) { std::cerr << "γ > 2/√3 is required\n"; return relation;