Implementing an algorithm to compute pi in R

Question

I am trying to implement a variation of the Brent-Salamin algorithm in R. It works well for the first 25 iterations, but then, it behaves unexpectedly, returning negative results.

The algorithm I want to implement is:

initial values:
x_0 = 1; y_0 = 1/sqrt(2); z_0 = 1/2

x_n = (x_n-1 + y_n-1)/2 
y_n = sqrt(x_n-1 * y_n-1)
z_n = z_n-1 - 2^n * (x_n^2-y_n^2)

p_n = (2 * x_n^2) / z_n

where n is the current iteration.

A more beautifully formatted formula is here.

The code I figured out is:

mypi <- function(n){

  x = 1
  y = 1/sqrt(2)
  z = 1/2
  iteration = 0

  while(iteration < n){
    iteration = iteration + 1

    newx = (x + y) / 2
    y = sqrt(x * y)
    x = newx
    z = z-(2^iteration * (x^2 - y^2))
    p = (2 * x^2) / z
  }

  return(p)
}

Output:

> mypi(10)
[1] 3.141593
> mypi(20)
[1] 3.141593
> mypi(50)
[1] -33.34323

So as I am new to R, is there a bug in my code or is it the algorithm?

After playing around with your code for about 20 minutes, I don't have an exact answer, and I'm not sure there is one, other than the explanation is probably due to the limitations of floating point arithmetic. I got it to work for almost 50 iterations before I started seeing negative numbers. I think after so many iterations, the 2^iteration term in z has become so large, and the x^2 - y^2 term has become so small, that rounding etc. starts to kick in. The negative number you see is just an artifact. — Tim Biegeleisen
– Tim Biegeleisen, Commented Dec 27, 2017 at 15:48
I edited the question to fix the i issue, and also to provide some example output. — Claus Wilke
– Claus Wilke, Commented Dec 27, 2017 at 15:52
@TimBiegeleisen it's pretty surprising though. In my implementation there is a factor p which multiplies itself by two each iteration. Is it perhaps something wrong with the ^ operator? — AdamO
– AdamO, Commented Dec 27, 2017 at 23:39

AdamO · Accepted Answer · 2017-12-27 15:51:12Z

13

Your code simply messes up because it does not agree with the algorithm as written in the wiki page. A correct version looks like this:

mypi <- function(n){

  x = 1
  y = 1/sqrt(2)
  z = 1/4
  p <- 1

  iteration = 0

  while(iteration < n){
    iteration = iteration + 1

    newx = (x + y) / 2
    y = sqrt(x * y)
    # x = newx
    # z = z-(2^iteration * (x^2 - y^2))
    z = z- p* (x-newx)^2
    p = 2*p
    x = newx
  }

  (newx + y)^2/(4*z)
}

Gives

> mypi(10)
[1] 3.141593
> mypi(20)
[1] 3.141593
> mypi(50)
[1] 3.141593

answered Dec 27, 2017 at 15:51

AdamO

5,0832 gold badges33 silver badges47 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

Tim Biegeleisen Over a year ago

Nice catch...I wouldn't have thought to check for the accuracy of the code itself +1

amadeusamadeus Over a year ago

Thank you for your answer. As I wrote in the first sentence of my question, I was trying to implement a (given) variation of the algorithm, but because of the results I wasn't sure about its accuracy. This is why I asked at the end of my question if there are any bugs in the R code or if the problem is the algorithm itself! So if the variation is wrong for more than ~25 iterations, but the implementation of the variation is ok, this would answer my question and be very helpful to me. So did you find any problems within the implementation of the 'wrong' algorithm, too?

AdamO Over a year ago

@amadeusamadeus If you want to debug the algorithm, I think the folks at math.stackexchange.com should be able to help you. I am not 100% certain how you derived this... it kind of looks like Householder's algorithm for some arctangent value that gives a fraction of $\pi$. Anyway, your algorithm matches the beautifully formatted formula, so no issue there. One last thing to consider: most algorithms terminate not based on running fixed iterations like you've done, but when the value achieves a target precision. 50 iterations could be reaching the limits of floating point precision...

AdamO Over a year ago

@amadeusamadeus to check this, I would return the various arithmetic combinations (x,y values) in an array. Geometrically, good approximations will come from functions that are very well behaved around the target value.

Collectives™ on Stack Overflow

Implementing an algorithm to compute pi in R

1 Answer 1

4 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Your Answer

Sign up or log in

Post as a guest

Related