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..5c551a26f
--- /dev/null
+++ b/doc/reference_pslq.qbk
@@ -0,0 +1,134 @@
+[/
+ 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 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:
+
+[$../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 induces a computation lasting 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);
+
+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]
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..b1671d1e0
--- /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 << symbol << " ≈ " << scf << "\n";
+ auto ccf = centered_continued_fraction(x);
+ std::cout << symbol << " ≈ " << ccf << "\n";
+ auto lur = luroth_expansion(x);
+ std::cout << symbol << " ≈ " << lur << "\n";
+ auto eng = engel_expansion(x);
+ 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";
+ 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 << "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";
+ 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 << "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";
+ 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 " << symbol << " 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 of " << symbol << "\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..c7c6998f1
--- /dev/null
+++ b/example/representations/constants.cpp
@@ -0,0 +1,29 @@
+#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");
+ analyze(dottie(), "d");
+ analyze(laplace_limit(), "λ");
+ analyze(reciprocal_fibonacci(), "ψ");
+}
diff --git a/example/representations/dirichlet_l.cpp b/example/representations/dirichlet_l.cpp
new file mode 100644
index 000000000..14910fc10
--- /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(900);
+ 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)");
+}
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 = " <
+#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/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/soldner.cpp b/example/representations/soldner.cpp
new file mode 100644
index 000000000..d212f8510
--- /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(1000);
+ using Real = mpfr_float;
+ Real s = soldner();
+ analyze(s, "μ");
+}
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..cbc838ac0
--- /dev/null
+++ b/include/boost/multiprecision/pslq.hpp
@@ -0,0 +1,736 @@
+/*
+ * 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