Calculate pi in R quikstart
R中也可以用Rmpfr
包实现多精度的计算。例如,我们可以用如下代码实现AGM算法计算Pi到小数点后256位。
library(Rmpfr)
<- function(prec=256, itermax = 100, verbose=TRUE) {
piMpfr <- mpfr(2, prec) # '2' as mpfr number
m2 ## -> all derived numbers are mpfr (with precision 'prec')
<- m2 + sqrt(m2) # 2 + sqrt(2) = 3.414..
p <- sqrt(sqrt(m2)) # 2^ {1/4}
y <- (y+1/y) / m2
x <- 0L
it repeat {
<- p
p.old <- it+1L
it <- p * (1+x) / (1+y)
p if(verbose) cat(sprintf("it=%2d, pi^ = %s, |.-.|/|.|=%e\n",
formatMpfr(p, min(50, prec/log2(10))), 1-p.old/p))
it, if (abs(p-p.old) <= m2^(-prec))
break
if(it > itermax) {
warning("not converged in", it, "iterations") ; break
}
## else
<- sqrt(x)
s <- (y*s + 1/s) / (1+y)
y <- (s+1/s)/2
x }
p}
piMpfr(prec = 256)
## it= 1, pi^ = 3.1426067539416226007907198236183018919713562462772, |.-.|/|.|=-8.642723e-02
## it= 2, pi^ = 3.1415926609660442304977522351203396906792842568645, |.-.|/|.|=-3.227958e-04
## it= 3, pi^ = 3.1415926535897932386457739917571417940347896238675, |.-.|/|.|=-2.347934e-09
## it= 4, pi^ = 3.1415926535897932384626433832795028841972241204666, |.-.|/|.|=-5.829228e-20
## it= 5, pi^ = 3.1415926535897932384626433832795028841971693993751, |.-.|/|.|=-1.741826e-41
## it= 6, pi^ = 3.1415926535897932384626433832795028841971693993751, |.-.|/|.|=0.000000e+00
## 1 'mpfr' number of precision 256 bits
## [1] 3.141592653589793238462643383279502884197169399375105820974944592307816406286163