Re: GMP precision problem - only getting 15 places of PI

[mensanator]

nmenuj@gmail.com wrote:
> The output is correct after incorporating the change David pointed out.
> Mensanator, please see below for the changed program using precomputed
> mpf's for the float literals in the loop. There is some speedup I
> think. I thought the type promotion for a const would occur just once
> in a C++ program, had no idea it runs inside the loop code.
>
> I appreciate your comments.

I couldn't get GMP to link with my g++ compiler under Cygwin,
so I couldn't actually try your program. But I noticed the
precomputed problem in both the C and Python using GMP.
The time saved will depend on how many iterations are in your
loop.

For instance, in C, I was running a Collatz sequence which
looped 2.5 million times. The precomputed constants changed
the run time from 300 sec to 17 sec. Another thing I noticed
was although Java's BigNumber library outperformed Python/GMP
on small problems, it rapidly fell behind when serious problems
had to be run. Serious being >2.5 million loops using numbers
that _start_ with over 50000 decimal digits.

>
> Before (changing constants to mpfs):
> real    0m1.236s
> user    0m1.191s
> sys     0m0.001s
>
> After:
> real    0m1.225s
> user    0m1.177s
> sys     0m0.002s
>
> I know the loop bounds are way more than required for just 100 digits,
> but I suppose they are good for performance checking.

I wrote a Python program to compute Pi digits. To eliminate
floating point imprecision, I used a Pi formula based on
rationals (since GMP can do unlimited precision rationals):

# better pi/2 = 1 + 1/3 + (1*2)/(3*5) + (1*2*3)/(3*5*7) + ...

By keeping the running value as a rational and discarding
the floating point version (after checking if the number of
digits has converged to the requested number), I don't have
to preset my FP precision and hope it's enough. Converting
the rational to a FP will automatically set the FP to the
necessary precision.

import gmpy
# unlimited precision math library
def pialso(n,b):
# input number of digits (n) in requested base (b)
        p = gmpy.mpq(1,1)
        # gmpy rationals are unlimited precision
        sn = 1
        sd = 3
        c = gmpy.mpq(sn,sd)
        # create the next term to be summed
        num = gmpy.mpf(p.numer())
        # seperately convert the numerator
        den = gmpy.mpf(p.denom())
        # and denominator to gmpy floats
        f = (num/den) * 2
        # to get unlimited precision float
        olds = gmpy.fdigits(f,b,n,0,1,2)
        # extract the requested digits and do base conversion
        done = 0
        while done==0:
                p = p + c
                # sum next term of pi equation
                sn += 1
                # set numerator
                sd += 2
                # and demoniator
                c = c * gmpy.mpq(sn,sd)
                # to create term for next iteration
                num = gmpy.mpf(p.numer())
                # meanwhile, convert this iteration
                den = gmpy.mpf(p.denom())
                f = (num/den) * 2
                # to an unlimited precision float
                s = gmpy.fdigits(f,b,n,0,1,2)
                if s[0]==olds[0]:
                # we're done when the number of digits
                        done = 1
                # we want stops changing
                else:
                # otherwise, keep iterating until we reach the
                        olds = s
                # the desired convergence
        print s[0],sn
        # prints the digits and how many iterations it took

pialso(100,10)

31415926535897932384626433832795028841971693993751
05820974944592307816406286208998628034825342117067

327

With this algorithm, it took 327 iterations to converge to
100 digits.

Running again to compute the first 100 base 12 digits, it
takes more iterations since base 12 digits are 'bigger'
than base 10 digits.

pialso(100,12)

3184809493b918664573a6211bb151551a05729290a7809a49
2742140a60a55256a0661a03753a3aa54805646880181a3682

353

So not only did I not have to guess how much precision to
assign, I didn't have to guess how big to make the loop either.

>
> #include <iostream>
> #include <gmpxx.h>
>
> using namespace std;
> const int PREC=20000;
>
> mpf_class pival(0.0, PREC);
> mpf_class d2(1.0, PREC), d3(1.0,PREC);
>
> const mpf_class const_1(1.0, PREC);
> const mpf_class const_2(2.0, PREC);
> const mpf_class const_3(3.0, PREC);
> const mpf_class const_4(4.0, PREC);
> const mpf_class const_9(9.0, PREC);
> const mpf_class const_16(16.0, PREC);
> const mpf_class const_81(81.0, PREC);
>
> int main()
> {
>         int i = 0;
>         mpf_set_default_prec(PREC);
>         d2 = const_1/const_2;
>         pival = d2 - (d2/const_4)/const_3;
>         for (i=5;i<20000;i+=4)
>         {
>                 d2 /= const_16;
>                 pival = pival + d2/i;
>                 pival = pival - (d2/const_4)/(i + 2);
>         }
>         cout << "def prec " << mpf_get_default_prec() << '\n';
>         d3 = mpf_class(1.0, PREC)/const_3;
>         pival = pival +  d3 - (d3/const_9)/const_3;
>         for (i=5;i<20000;i+=4)
>         {
>                 d3 /= const_81;
>                 pival = pival + d3/i;
>                 pival = pival - (d3/const_9)/(i + 2);
>         }
>         cout.precision(100);
>         cout << pival * 4.0 << '\n';
>         return 0;
> }