Skip to content

problem with the precision of numbers

8 messages · kayj, (Ted Harding), William Dunlap +2 more

kayj
# 
Hi All,

I was wondering if R can deal with high precsion numbers, below is an
example that I tried on using R and Maple where I got different results. I
was not able to use the r-bc package in R, instead I used the Rmpfr
package, below are both R and Maple results
+ p<-as.numeric(i)
+ c<-choose(80,i)
+ s=s+((j^i)*c*(1-(i+1)*1/200)^200)
+  
+ }
1 'mpfr' number of precision  500   bits 
[1]
4.6484825769379918039202329231788093343835337105287241996630873753794810915791302829474613042869191092322033403649086087872119205043462728841279067042348e-6


Maple result

6.656852479*10^(-20)


why are the two results different?

I appreciate your help
(Ted Harding)
#
On 24-Jan-10 22:24:10, kayj wrote:
The result from bc (run "raw" and not from R, and assuming I
programmed it correctly) is (with some formatting for clarity):

s
0.
0000000000 0000000000 0000000000 0000000000 0000000000
0000000000 0000000000 0000000000 0000049078 2995761647
7217765991 5850974377 7588454525 1985604552 6041517025
6783796473 3395739785 9774327566 5079076761 9753976534

i.e.
4.90782995761647.......... * 10^(-86)

The 'bc' program was:

scale=200
define f(x) {
  if(x <= 1) return (1);
  return ( f(x-1)*x );
}
define choose(n,m) {
  return ( f(n)/(f(m)*f(n-m)) )
}
s=0
j=-1
for(i=0;i<=200;i++){
  c=choose(200,i);
  s=s+((j^i)*c*(1-(i+1)*1/200)^200);
}
s

(and, despite the recursive definition, the loop only took
about 3 seconds).

So my attempt came out different from your mpfr and your Maple.
It might be worth someone checking my 'bc' code! And, just in
case anyone suspects the accuracy of numbers like f(200) = 200!,
in R I get:

   print(lgamma(201),15)
  # [1] 863.231987192405

and in bc:
  l(f(200))
  #     863.2319871924054734957......................

so there is agreement there.

Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 24-Jan-10                                       Time: 23:54:54
------------------------------ XFMail ------------------------------
William Dunlap
# 
You are computing (1-(i+1)*1/200)^200 with ordinary
32-bit integers and 64-bit double precision numbers.
The same goes for choose(80,i).  Try something like

f <- function (precision) 
{
    require(Rmpfr)
    s <- mpfr(0, precision)
    j <- mpfr(-1, precision)
    c <- mpfr(1, precision)
    for (i in 0:80) {
        i <- mpfr(i, precision)
        s <- s + ((j^i) * c * (1 - (i + 1) * 1/200)^200)
        c <- (c * (80 - i))/(i + 1)
    }
    s
}
1 'mpfr' number of precision  500   bits 
[1] 6.656852478966203769967328e-20

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
288905p1288905.html
(Ted Harding)
#
On 25-Jan-10 00:12:29, William Dunlap wrote:
Bill, using your code strategy I agree with your mpfr result
(and with Maple) using 'bc:

s=0
j = -1
c = 1
for (i=0; i <= 80; i++ ){
        s = s + ((j^i) * c * (1 - (i + 1) * 1/200)^200)
        c = (c * (80 - i))/(i + 1)
    }
s
0.
0000000000 0000000006 6568524789 6620376996 7327577118
2135789699 5101929940 0816838328 9634271996 1257038416 
6039147205 2215051658 2822460868 2286772321 8214397124
2739979506 9906258378 1667597096 6899329803 4460676207

(time to look and see what I did wrong before ...)
Ted.
--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 25-Jan-10                                       Time: 00:47:53
------------------------------ XFMail ------------------------------
Gabor Grothendieck
# 
I have added a FAQ to the r-bc home page which should help address
your problems in getting bc to work.
On Sun, Jan 24, 2010 at 5:24 PM, kayj <kjaja27 at yahoo.com> wrote:
kayj
# 
Hi All,

thank you all for your help. I have tried Bill's script and it works! so I
am trying to think what was the problem and it looks like it i sthe
precision. so you defined a function of the precision and evaluates at
precision=500. Bill, I was wondering how did you choose 500? is it
arbitrary?

I have looked at the Rmpfr documentation and they defined the mfr object to
be
 Usage
mpfr(x, precBits, base = 10)
Const(name = c("pi", "gamma", "catalan"), prec = 120L)

so is the precision defined in Bill's script the same as precBits?

my bc package works , I have used bc ( from run "raw" and not from R) so is
it possible to use it from R?

thanks again for your help
William Dunlap
# 
The problem was that you were computing several quantities
in the machine's native double precision format, not using
the mpfr format.   I chose to parameterize the function
with precision so you could see how the precision affected
the result.  I showed precision=500 mainly because that is
what you used for the variables that you did use in mpfr
format.  I did not do any analysis to see what precision
would be needed to get a given final precision, but that
would be possible.  Here is how precision affects the results:
require(Rmpfr)
    s <- mpfr(0, precision)
    j <- mpfr(-1, precision)
    c <- mpfr(1, precision)
    for (i in 0:80) {
        i <- mpfr(i, precision)
        s <- s + ((j^i) * c * (1 - (i + 1) * 1/200)^200)
        c <- (c * (80 - i))/(i + 1)
    }
    s
}
work.
1 'mpfr' number of precision  2   bits 
[1] 7.2535549176877750482e24
1 'mpfr' number of precision  4   bits 
[1] 1099511627776
1 'mpfr' number of precision  8   bits 
[1] -99090432
1 'mpfr' number of precision  16   bits 
[1] 1815744
1 'mpfr' number of precision  32   bits 
[1] 40.655612170696258545
1 'mpfr' number of precision  64   bits 
[1] 3.3910635476710543806e-9
1 'mpfr' number of precision  128   bits 
[1] 6.6568524698115877284e-20
1 'mpfr' number of precision  256   bits 
[1] 6.6568524789662037700e-20
1 'mpfr' number of precision  512   bits 
[1] 6.6568524789662037700e-20
1 'mpfr' number of precision  1024   bits 
[1] 6.6568524789662037700e-20 

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
288905p1289483.html
Hans W Borchers
# 
kayj <kjaja27 <at> yahoo.com> writes:
Your question was "I was wondering if it is possible to increase the precision
using R", but you never told which accuracy you really need, As a rule of thumb,
the number of correct decimal digits is about 0.3 * precision, so for 50 digits
you could set the precision to 256 to Be on the safe side.

And if you are wondering how to compute the binomial coefficients with 'Rmpfr',
here's one possibility to do that:

    require(Rmpfr)

    mpfrChoose <- function(a, b, prec=128) {
        m1 <- mpfr(1, prec=prec)
        # r <- gamma(m1+a) / (gamma(m1+b) * gamma(m1+a-b))
        # if (is.whole(r)) return(r) else return(round(r))
        gamma(m1+a) / (gamma(m1+b) * gamma(m1+a-b))
    }

An advantage of using 'Rmpfr' is that the power of R can be applied, for
example vectorization, so avoid 'for' loops if possible:

    pr  <- 256
    m_1 <- mpfr(-1, pr)
    m1  <- mpfr(1, pr)

    i <- mpfr(0:80, pr)
    s <- sum((m_1^i) * mpfrChoose(80, i, prec=pr) * (m1-(i+1)*1/200)^200)

    print(s, digits=50)
    # 1 'mpfr' number of precision  256   bits
    # [1] 6.6568524789662037699673275771182135789699510194000e-20

Hans Werner