2

I am trying to calculate pi for one of my university project using Ramanujan's formula for arbitrary number of digits after the floating point. For the job, I am using boost::multiprecision library that is simply wrapper around mpfr and mpir I already have installed on my machine.

So far so good, but I am either missing something or doing something wrong since my calculate function (the one calculating each iteration of the sum in the formula) throws StackOverflow exception every now and then. (when it starts it is consistent)

This is how it looks

boost::multiprecision::mpf_float calculatePi(int start, int end)
{
    boost::multiprecision::mpf_float partition = 0;
    for (; start < end; ++start)
    {
        boost::multiprecision::mpf_float n = 
            factorial(4 * start) * boost::multiprecision::mpf_float(1103 + 26390 * start);
        boost::multiprecision::mpf_float d = 
            boost::multiprecision::pow(factorial(start), 4) * pow((boost::multiprecision::mpf_float)396, 4 * start); <-- this is where stackoverflow exception is thrown

        partition += (n / d);
    }

    return (2 * boost::multiprecision::sqrt((boost::multiprecision::mpf_float)2) / 9801) * partition;
}

I admit I am far too green into whole arbitrary precision calculation/conversions and libraries so I might be missing something.

I can't show you full call stack because it is too long but after

ParallelPi.exe!boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>::**do_multiplies**<boost::multiprecision::detail::expression<boost::multiprecision::detail::function,boost::multiprecision::detail::number_kind_floating_pointpow_funct<boost::multiprecision::backends::gmp_float<0> >,boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>,int,void>,boost::multiprecision::detail::function>(const boost::multiprecision::detail::expression<boost::multiprecision::detail::function,boost::multiprecision::detail::number_kind_floating_pointpow_funct<boost::multiprecision::backends::gmp_float<0> >,boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>,int,void> & e, const boost::multiprecision::detail::function & __formal) Line 1754    C++

Then there are 3603 consequent calls to

ParallelPi.exe!boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>::**operator=**<boost::multiprecision::detail::function,boost::multiprecision::detail::number_kind_floating_pointpow_funct<boost::multiprecision::backends::gmp_float<0> >,boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>,int,void>(const boost::multiprecision::detail::expression<boost::multiprecision::detail::function,boost::multiprecision::detail::number_kind_floating_pointpow_funct<boost::multiprecision::backends::gmp_float<0> >,boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>,int,void> & e) Line 216  C++
ParallelPi.exe!boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>::number<boost::multiprecision::backends::gmp_float<0>,1><boost::multiprecision::detail::function,boost::multiprecision::detail::number_kind_floating_pointpow_funct<boost::multiprecision::backends::gmp_float<0> >,boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>,int,void>(const boost::multiprecision::detail::expression<boost::multiprecision::detail::function,boost::multiprecision::detail::number_kind_floating_pointpow_funct<boost::multiprecision::backends::gmp_float<0> >,boost::multiprecision::number<boost::multiprecision::backends::gmp_float<0>,1>,int,void> & e, void * __formal) Line 324   C++

only to lead into stackoverflow exceptions few calls later in

ParallelPi.exe!boost::multiprecision::backends::gmp_float<0>::precision() Line 630  C++

Funny as it sounds, I don't recall changing calculatePi function before the time it was "working".

Could you guys help me decipher what is going on cause I seem to be lost?

My default precision for boost::multiprecision::mpf_float is 150 digits after floating point which previously was not a problem (I recall I calculated some number with 10 000 floating point precision with not SO or whatsoever)


Requests

1.Code of my factorial

boost::multiprecision::mpf_float factorial(int n)
{
    boost::multiprecision::mpf_float fact = 1;
    for (int i = 1; i <= n; ++i)
        fact *= i;

    return fact;
}
kuskmen
  • 3,648
  • 4
  • 27
  • 54
  • 1
    I do not have an answer, but you can greatly improve the performance be reusing some terms. `[(n+1)!]^4 = (n+1)^4 * (n!)^4`, `396^[4(n+1)] = 394^4 * 396^(4n)`. `150!^4` has over 1000digits, `394^(4n)` over 1500. Maybe that will hide the issue. – Quimby May 21 '19 at 20:43
  • 2
    I am currently in process in making the dumbest solution work, then proceed with the optimizations, nevertheless, thanks for the suggestion I will keep it in mind when I reach that far. Regarding digits overflow, I did not have this problem before, plus this is happening on the very first iteration of the loop when start is 0. – kuskmen May 21 '19 at 20:47
  • @kuskmen: Where is your definition of factorial? – user14717 May 22 '19 at 02:45
  • 2
    Probably not related to your problem, but i suggest you use mpfr_float instead of mpf_float. It also has the advantage that it already provides pi, so you can compare with your result. – Marc Glisse May 22 '19 at 05:57
  • @MarcGlisse changing mpf_float to mpfr_float "fix" the problem and now I don't get SO exception, but I suppose the question remains open why there is SO with mpf_float – kuskmen May 22 '19 at 06:11

1 Answers1

0

This is a bug in recent versions of Boost.Multiprecision - I'm trying to figure this out now, but see https://github.com/boostorg/multiprecision/issues/164.

Note that using mpfr rather than mpf would fix the issue.

John Maddock
  • 111
  • 1
  • 2