7
\$\begingroup\$

For a while now I have been working to construct a program to calculate a lot of mathematical constants. Before I explain, here's the code:

Code

#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/math/special_functions.hpp>
#include <complex>
#include <iostream>

/*
    Only change PRECISION
*/

const unsigned long long PRECISION = 100, INACC = 4;

typedef boost::multiprecision::number<
    boost::multiprecision::cpp_dec_float<PRECISION + INACC> > arb;

arb alladi_grinstead()
{
    arb c;

    for(unsigned long long n = 1;; ++n) {
        const arb last = c;
        c += (boost::math::zeta<arb>(n + 1) - 1) / n;
        if(c == last) break;
    }

    return exp(c - 1);
}

arb aperys()
{
    return boost::math::zeta<arb>(3);
}

arb buffons()
{
    return 2 / boost::math::constants::pi<arb>();
}

arb catalans()
{
    const arb PI = boost::math::constants::pi<arb>();
    return 0.125 * (boost::math::trigamma<arb>(0.25) - PI * PI);
}

arb champernowne(const unsigned long long b = 10)
{
    arb c;
    for(unsigned long long n = 1;; ++n) {
        arb sub;

        for(unsigned long long k = 1; k <= n; ++k) {
            sub += floor(log(k) / log(b));
        }

        const arb last = c;
        c += n / pow(b, n + sub);
        if(c == last) break;
    }

    return c;
}

arb delian()
{
    return boost::math::cbrt<arb>(2);
}

arb dottie()
{
    arb x;
    std::string precomp, postcomp;

    for(x = 1; x == 1 || precomp != postcomp;) {
        precomp = static_cast<std::string>(x);
        precomp.resize(PRECISION);

        x -= (cos(x) - x) / (-sin(x) - 1);

        postcomp = static_cast<std::string>(x);
        postcomp.resize(PRECISION);
    }

    return x;
}

arb e()
{
    return boost::math::constants::e<arb>();
}

arb erdos_borwein()
{
    arb e;

    for(unsigned long long n = 1;; ++n) {
        const arb last = e;
        e += 1 / (pow(static_cast<arb>(2), n) - 1);
        if(e == last) break;
    }

    return e;
}

arb euler_mascheroni()
{
    return boost::math::constants::euler<arb>();
}

arb favard(const unsigned long long r = 2)
{
    if(r % 2 == 0) {
        return (-4 / boost::math::constants::pi<arb>()) *
            (pow(-2, static_cast<arb>(-2) * (r + 1)) / boost::math::tgamma<arb>(r + 1)) *
            (boost::math::polygamma<arb>(r, static_cast<arb>(0.25)) - boost::math::polygamma<arb>(r, static_cast<arb>(0.75)));
    } else {
        return (4 / boost::math::constants::pi<arb>()) *
            ((1 - pow(2, -(static_cast<arb>(r) + 1))) * boost::math::zeta<arb>(r + 1));
    }
}

arb gauss()
{
    const arb ROOT_TWO = boost::math::constants::root_two<arb>(),
        PI = boost::math::constants::pi<arb>();
    return pow(boost::math::tgamma<arb>(0.25), 2)
        / (2 * ROOT_TWO * pow(PI, 3 / static_cast<arb>(2)));
}

arb gelfond_schneider()
{
    return pow(2, boost::math::constants::root_two<arb>());
}

arb gelfonds()
{
    return pow(boost::math::constants::e<arb>(), boost::math::constants::pi<arb>());
}

arb giesekings()
{
    return (9 - boost::math::trigamma<arb>(2 / static_cast<arb>(3)) +
        boost::math::trigamma<arb>(4 / static_cast<arb>(3))) /
        (4 * boost::math::constants::root_three<arb>());
}

arb glaisher_kinkelin()
{
    return boost::math::constants::glaisher<arb>();
}

arb golden_ratio()
{
    return boost::math::constants::phi<arb>();
}

std::complex<arb> i()
{
    return std::complex<arb>(0,1);
}

arb inverse_golden_ratio()
{
    return boost::math::constants::phi<arb>() - 1;
}

arb khinchin()
{
    return boost::math::constants::khinchin<arb>();
}

arb khinchin_levy()
{
    return pow(boost::math::constants::pi<arb>(), 2) / (12 * log(static_cast<arb>(2)));
}

arb kinkelin()
{
    return 1 / static_cast<arb>(12) - log(boost::math::constants::glaisher<arb>());
}

arb knuth()
{
    return (1 - (1 / boost::math::constants::root_three<arb>())) / 2;
}

arb levys()
{
    return exp(pow(boost::math::constants::pi<arb>(), 2) / (12 * log(static_cast<arb>(2))));
}

arb liebs()
{
    return (8 * boost::math::constants::root_three<arb>()) / 9;
}

arb lochs()
{
    return (6 * log(static_cast<arb>(2)) * log(static_cast<arb>(10))) /
        pow(boost::math::constants::pi<arb>(), 2);
}

arb niven()
{
    arb c;
    for(unsigned long long j = 2;; ++j) {
        const arb last = c;
        c+= 1 - 1/boost::math::zeta<arb>(j);
        if(c == last) break;
    }
    return c + 1;
}

arb nortons()
{
    const arb PI = boost::math::constants::pi<arb>(),
        EULER = boost::math::constants::euler<arb>(),
        GLAISHER = boost::math::constants::glaisher<arb>(),
        PI_SQR = pow(PI, 2),
        LOG_TWO = log(static_cast<arb>(2));
    return -((PI_SQR - 6 * LOG_TWO * (-3 + 2 * EULER + LOG_TWO + 24 * log(GLAISHER) - 2 * log(PI))) / PI_SQR);
}

arb omega()
{
    arb omega;
    std::string precomp, postcomp;

    for(omega = 0; omega == 0 || precomp != postcomp;) {
        precomp = static_cast<std::string>(omega);
        precomp.resize(PRECISION);

        omega -= ((omega * exp(omega)) - 1) /
            (exp(omega) * (omega + 1) - ((omega + 2) * (omega * exp(omega) - 1) / ((2 * omega) + 2)));

        postcomp = static_cast<std::string>(omega);
        postcomp.resize(PRECISION);
    }

    return omega;
}

arb one()
{
    return 1;
}

arb pi()
{
    return boost::math::constants::pi<arb>();
}

arb plastic_number()
{
    return (boost::math::cbrt<arb>(108 + 12 * sqrt(static_cast<arb>(69))) + 
        boost::math::cbrt<arb>(108 - 12 * sqrt(static_cast<arb>(69)))) / static_cast<arb>(6);
}

arb pogsons()
{
    return pow(100, 1 / static_cast<arb>(5));
}

arb polyas_random_walk()
{
    const arb PI = boost::math::constants::pi<arb>(),
        PI_CBD = pow(PI, 3),
        ROOT_SIX = sqrt(static_cast<arb>(6));
    return 1 - 1/((ROOT_SIX / (32 * PI_CBD)) *
        boost::math::tgamma<arb>(1 / static_cast<arb>(24)) *
        boost::math::tgamma<arb>(5 / static_cast<arb>(24)) *
        boost::math::tgamma<arb>(7 / static_cast<arb>(24)) *
        boost::math::tgamma<arb>(11 / static_cast<arb>(24)));
}

arb porters()
{
    const arb PI = boost::math::constants::pi<arb>(),
        GLAISHER = boost::math::constants::glaisher<arb>();
    return ((6 * log(static_cast<arb>(2)) * (48 * log(GLAISHER) - log(static_cast<arb>(2)) - 4 * log(PI) - 2))
        / pow(PI, 2)) - (1 / static_cast<arb>(2));
}

arb prince_ruperts_cube()
{
    return (3 * boost::math::constants::root_two<arb>()) / 4;
}

arb pythagoras()
{
    return boost::math::constants::root_two<arb>();
}

arb robbins()
{
    const arb PI = boost::math::constants::pi<arb>(),
        ROOT_TWO = boost::math::constants::root_two<arb>(),
        ROOT_THREE = boost::math::constants::root_three<arb>();
    return ((4 + 17 * ROOT_TWO - 6 * ROOT_THREE - 7 * PI) / 105)
        + (log(1 + ROOT_TWO) / 5) + ((2 * log(2 + ROOT_THREE)) / 5);
}

arb sierpinski_k()
{
    const arb PI = boost::math::constants::pi<arb>(),
        E = boost::math::constants::e<arb>(),
        EULER = boost::math::constants::euler<arb>();
    return PI * log((4 * pow(PI, 3) * pow(E, 2 * EULER)) / pow(boost::math::tgamma<arb>(0.25), 4));
}

arb sierpinski_s()
{
    const arb PI = boost::math::constants::pi<arb>(),
        E = boost::math::constants::e<arb>(),
        EULER = boost::math::constants::euler<arb>();
    return log((4 * pow(PI, 3) * pow(E, 2 * EULER)) / pow(boost::math::tgamma<arb>(0.25), 4));
}

arb silver_ratio()
{
    return boost::math::constants::root_two<arb>() + 1;
}

arb theodorus()
{
    return boost::math::constants::root_three<arb>();
}

arb twenty_vertex_entropy()
{
    return (3 * boost::math::constants::root_three<arb>()) / 2;
}

arb weierstrass()
{
    const arb PI = boost::math::constants::pi<arb>(),
        E = boost::math::constants::e<arb>();
    return (pow(2, static_cast<arb>(1.25)) * sqrt(PI) * pow(E, PI / 8)) /
        pow(boost::math::tgamma<arb>(0.25), 2);
}

arb wylers() {
    const arb PI = boost::math::constants::pi<arb>();
    return (9 / (8 * pow(PI, 4))) * pow(pow(PI, 5) / 1920, 0.25);
}

arb zero()
{
    return 0;
}

int main()
{
    std::cout << std::fixed << std::setprecision(PRECISION)
        << wylers() << '\n';
    return 0;
}

Code breakdown

From top to bottom (-ish):

  1. I use Boost for the following purposes:

    1. To allow for arbitrary precision, via cpp_dec_float.
    2. To take advantage of Boost's pre-built constants - such as Pi, the Khinchin constant, etc.
    3. To take advantage of Boost's "special functions."
  2. To use Boost's cpp_dec_float, I define a type arb. arb's precision is determined at compile-time, so I use const PRECISION. PRECISION is meant to be manually changed depending on what's needed. I also have the variable INACC, which adds precision to arb. By having INACC, I prevent larger inaccuracy (last ≥4 digits), but the constants may still be inaccurate in the last two digits because of rounding. So it goes.

  3. Each constant's function is independent of one another - and this is how I want it to be. I want it to be so that someone can copy the function code, the declaration of the arbitrary type (arb), and the required headers - that's it.

  4. Often, some constants require summation or successive approximation (such as Newton's method or Halley's method). To achieve this, I use a for-loop. The problem: it will go forever (& exceeding the accuracy of PRECISION) if unchecked. This results in two situations:

    1. I use const arb last to compare the variable before and after calculations. You can see this in action in alladi_grinstead().
    2. I use std::strings, by statically casting the variable, and then resizing it to PRECISION. You can see this in action in dottie(). The reason I have to use this method is because without resizing the variable, it will work with extended precision (beyond PRECISION, because of INACC) and continue forever (or, at least way longer). If I attempt to use the first method described, it doesn't work.
  5. Constants can be printed in conjunction with std::setprecision(PRECISION) - std::fixed isn't required, but it helps.

Conclusion and Misc.

For further detail on the mathematics I use, my web page expands on such things.

I tested this on Linux, using g++ with the -O3 flag. I have checked that all functions work properly to calculate the constants to any precision - well, except the Niven constant, which I couldn't find a way to verify its accuracy past 100 digits.


How can I optimize my code, by improving the code's logic, structure, readability, or by improving performance. My code goes beyond the 80-char width, which is annoying because it makes me scroll horizontally a lot, however I am not sure where to break it so that it still is readable.

Honestly, all improvements are useful.

\$\endgroup\$
  • \$\begingroup\$ Is that pow() supposed to be std::pow()? If it is, you'll need to include <cmath> and add using std::pow(); into the functions that use it (or some equivalent change). I would avoid std::pow() for small integer exponents - hand-code the multiplication instead. \$\endgroup\$ – Toby Speight Aug 23 '17 at 17:13
4
\$\begingroup\$

I see some things that may help you improve your program.

Use static const for performance

The way it works right now, each function recalculates the value every time. However, these are constants, right? Each value needs only to be calculated once. Doing so will use some memory, but as I'll show later, there are also ways to mitigate that. Here's a simple modification that shows the concept:

arb wylers() {
    static const arb PI = boost::math::constants::pi<arb>();
    static const arb result = (9 / (8 * pow(PI, 4))) * pow(pow(PI, 5) / 1920, 0.25);
    return result;
}
int main()
{
    auto wy = wylers();
    for (std::size_t trials{1'000'000}; trials; --trials) {
        wy = wylers();
    }
    std::cout << std::fixed << std::setprecision(PRECISION) << wy << '\n';
}

Note that I'm using C++14 in the test code (only) hence, the quote character as a digit separator. Other than that, it should compile and run in any C++11 compliant environment. This test code calculating Wyler's constant one million times took over 3 minutes on my machine with the original version but just 140 ms with this version because, of course, the value was only actually calculated once and the only overhead was involved in copying result one million times.

Use a template

There is a potential problem with the code above in that if the user wants, say, both a 50 digit and 104 digit version, there's not really a way to do that because it's a single function. I'd fix that by using a template. The rewritten version would look like this:

template<typename T>
T wylers() {
        static const T PI = boost::math::constants::pi<arb>();
        static const T result = (9 / (8 * pow(PI, 4))) * pow(pow(PI, 5) / 1920, 0.25);
        return result;
}

In addition to allowing the creation of multiple version, each with its own precision, the fact that its a template means that it generates no code at all unless it's actually invoked somewhere in the code. This makes the resulting binary smaller and faster.

"Eat your own dogfood"

One way to help assure that your software is usable is to use it yourself. The phrase "Eat your own dogfood" refers to a company using its own products. For example, why store an additional private copy of PI? So starting with the proposal above:

static const T PI = boost::math::constants::pi<arb>();
static const T result = (9 / (8 * pow(PI, 4))) * pow(pow(PI, 5) / 1920, 0.25);

we can instead write it like this:

template<typename T>
const T &wylers() {
        static const T result = (9 / (8 * pow(pi<T>(), 4))) * pow(pow(pi<T>(), 5) / 1920, 0.25);
        return result;
}

Using this version, our million calculations take 9 ms on my machine. Naturally, we also need the templated version of pi():

template<typename T>
T pi()
{
    return boost::math::constants::pi<T>();
}

In this case, I decided not to use the static const trick because I would expect that the boost version of a templated pi is probably already very efficient.

Note that this also requires that pi must be defined before any use which means your neat alphabetic order would be disturbed.

Return a const reference where practical

The current routines always make a copy even when none is needed. It would probably make more sense to return a const reference instead which would then only force a copy when one is really needed. For instance the boost code for number includes this operator:

cpp_dec_float& operator*=(const cpp_dec_float& v);

Because the other argument is const &, we might guess that the multiplication can be done without making a copy. Indeed, examining the source code for that function verifies this hunch. The effect is then that we can write things like this:

arb pi_5 = 1;
for (int i=5; i; --i)
    pi_5 *= pi<arb>();   // no need to make a copy of pi here

and they won't be forced to make multiple copies of pi.

Break templates into pieces where needed

For functions like omega(), the static const trick won't work because the value is calculated in a loop and can't be const. One simple way of dealing with this is by splitting the function into two pieces:

template<typename T>
T omega_helper()
{
        T omega;
        std::string precomp, postcomp;

        for(omega = 0; omega == 0 || precomp != postcomp;) {
                precomp = static_cast<std::string>(omega);
                precomp.resize(PRECISION);

                omega -= ((omega * exp(omega)) - 1) /
                        (exp(omega) * (omega + 1) - ((omega + 2) * (omega * exp(omega) - 1) / ((2 * omega) + 2)));

                postcomp = static_cast<std::string>(omega);
                postcomp.resize(PRECISION);
        }
        return omega;
}

template<typename T>
const T& omega()
{
        static const T omega = omega_helper<T>();
        return omega;
}

Add your constants to boost

Since you've invested some time and effort into this already, you might want to actually add the constants to boost (either your local copy or submit it for inclusion). That mechanism is detailed in the boost docs.

\$\endgroup\$
  • \$\begingroup\$ Thanks. I especially liked the thought of using templates, such that the constants can be computed independent of Boost's arbitrary types. \$\endgroup\$ – esote Jul 23 '17 at 19:40
  • \$\begingroup\$ A trick I sometimes use for one-off initializations like omega_helper() is to put the code into an immediately-executed lambda expression: static const T omega = []{ /* calculation */ return omega; }();. On a separate note, would it be worth declaring these functions as constexpr? \$\endgroup\$ – Toby Speight Aug 23 '17 at 17:08
1
\$\begingroup\$

This is a really interesting project. I don't have a whole lot to add, but here are some things that I see.

Put the condition in your loop

All of your for loops terminate when an iteration of the loop produces the same result as previously. So why not put the condition into the loop construct? Either write it as:

arb c = 0;
arb last = -1;
for(unsigned long long n = 1; c != last; ++n) {
    last = c;
    //... whatever calculation
}

Either that or write it as a do loop:

arb c = 0;
arb last = -1;
unsigned long long n = 1;
do {
    last = c;
    // Do calculation
    n++;
} while (last != c);

This makes the flow of control much more obvious to someone reading your code.

Do You Need These Functions?

There are a few functions that don't calculate anything and appear to be returning a constant, such as one(), pi(), and zero(). The only reason I can see to have these is to avoid the syntax needed to write them out. Why not replace it with a constant, like:

const auto kArbPi = boost::math::constants::pi<arb>();
\$\endgroup\$
  • \$\begingroup\$ By implementing your first suggestion, I was able to improve run-time by a noticeable margin. Regarding your 2nd suggestion, I know that it is rather redundant to have stuff like 1, 0, and pi as their own functions - and I only really did it for completeness' sake. \$\endgroup\$ – esote Jul 23 '17 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.