CODEWITHSUNDEEP
ralgorithm
amadeusamadeus
amadeusamadeus

Reputation: 121

Implementing an algorithm to compute pi in R

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?

Upvotes: 1

Views: 793

Answers (1)

AdamO
AdamO

Reputation: 4930

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

Upvotes: 12

Related Questions