巴比倫演算法

# 以 N = 4567, x = 70 為例:
options(digits=22)
N <- 4567
x <- 70
y <- N/x
for (i in 1:30){x <- (x+y)/2; y <- N/x; cat("i=",i,"x=",x,"y=",y,"\n")}
## i= 1 x= 67.62142857142856655628 y= 67.53776275483258473287 
## i= 2 x= 67.57959566313057564457 y= 67.57956976785553138143 
## i= 3 x= 67.57958271549304640757 y= 67.57958271549058792971 
## i= 4 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 5 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 6 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 7 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 8 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 9 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 10 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 11 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 12 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 13 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 14 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 15 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 16 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 17 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 18 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 19 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 20 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 21 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 22 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 23 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 24 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 25 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 26 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 27 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 28 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 29 x= 67.57958271549182427407 y= 67.57958271549181006321 
## i= 30 x= 67.57958271549182427407 y= 67.57958271549181006321

修改有效位數：

# 以 N = 4567, x = 70 為例:
options(digits=12)
N <- 4567
x <- 70
y <- N/x
for (i in 1:30){x <- (x+y)/2; y <- N/x; cat("i=",i,"x=",x,"y=",y,"\n")}
## i= 1 x= 67.6214285714 y= 67.5377627548 
## i= 2 x= 67.5795956631 y= 67.5795697679 
## i= 3 x= 67.5795827155 y= 67.5795827155 
## i= 4 x= 67.5795827155 y= 67.5795827155 
## i= 5 x= 67.5795827155 y= 67.5795827155 
## i= 6 x= 67.5795827155 y= 67.5795827155 
## i= 7 x= 67.5795827155 y= 67.5795827155 
## i= 8 x= 67.5795827155 y= 67.5795827155 
## i= 9 x= 67.5795827155 y= 67.5795827155 
## i= 10 x= 67.5795827155 y= 67.5795827155 
## i= 11 x= 67.5795827155 y= 67.5795827155 
## i= 12 x= 67.5795827155 y= 67.5795827155 
## i= 13 x= 67.5795827155 y= 67.5795827155 
## i= 14 x= 67.5795827155 y= 67.5795827155 
## i= 15 x= 67.5795827155 y= 67.5795827155 
## i= 16 x= 67.5795827155 y= 67.5795827155 
## i= 17 x= 67.5795827155 y= 67.5795827155 
## i= 18 x= 67.5795827155 y= 67.5795827155 
## i= 19 x= 67.5795827155 y= 67.5795827155 
## i= 20 x= 67.5795827155 y= 67.5795827155 
## i= 21 x= 67.5795827155 y= 67.5795827155 
## i= 22 x= 67.5795827155 y= 67.5795827155 
## i= 23 x= 67.5795827155 y= 67.5795827155 
## i= 24 x= 67.5795827155 y= 67.5795827155 
## i= 25 x= 67.5795827155 y= 67.5795827155 
## i= 26 x= 67.5795827155 y= 67.5795827155 
## i= 27 x= 67.5795827155 y= 67.5795827155 
## i= 28 x= 67.5795827155 y= 67.5795827155 
## i= 29 x= 67.5795827155 y= 67.5795827155 
## i= 30 x= 67.5795827155 y= 67.5795827155

亦可寫成函數的形式：

Babylon <- function(N,x) {
options(digits=16)
y <- N/x
for (i in 1:30){x <- (x+y)/2; y <- N/x; cat("i=",i,"x=",x,"y=",y,"\n")}
}

Babylon(6789,80)
## i= 1 x= 82.43125000000001 y= 82.35954204261127 
## i= 2 x= 82.39539602130563 y= 82.39538041961124 
## i= 3 x= 82.39538822045844 y= 82.39538822045769 
## i= 4 x= 82.39538822045807 y= 82.39538822045806 
## i= 5 x= 82.39538822045807 y= 82.39538822045806 
## i= 6 x= 82.39538822045807 y= 82.39538822045806 
## i= 7 x= 82.39538822045807 y= 82.39538822045806 
## i= 8 x= 82.39538822045807 y= 82.39538822045806 
## i= 9 x= 82.39538822045807 y= 82.39538822045806 
## i= 10 x= 82.39538822045807 y= 82.39538822045806 
## i= 11 x= 82.39538822045807 y= 82.39538822045806 
## i= 12 x= 82.39538822045807 y= 82.39538822045806 
## i= 13 x= 82.39538822045807 y= 82.39538822045806 
## i= 14 x= 82.39538822045807 y= 82.39538822045806 
## i= 15 x= 82.39538822045807 y= 82.39538822045806 
## i= 16 x= 82.39538822045807 y= 82.39538822045806 
## i= 17 x= 82.39538822045807 y= 82.39538822045806 
## i= 18 x= 82.39538822045807 y= 82.39538822045806 
## i= 19 x= 82.39538822045807 y= 82.39538822045806 
## i= 20 x= 82.39538822045807 y= 82.39538822045806 
## i= 21 x= 82.39538822045807 y= 82.39538822045806 
## i= 22 x= 82.39538822045807 y= 82.39538822045806 
## i= 23 x= 82.39538822045807 y= 82.39538822045806 
## i= 24 x= 82.39538822045807 y= 82.39538822045806 
## i= 25 x= 82.39538822045807 y= 82.39538822045806 
## i= 26 x= 82.39538822045807 y= 82.39538822045806 
## i= 27 x= 82.39538822045807 y= 82.39538822045806 
## i= 28 x= 82.39538822045807 y= 82.39538822045806 
## i= 29 x= 82.39538822045807 y= 82.39538822045806 
## i= 30 x= 82.39538822045807 y= 82.39538822045806

計算結果做比較：

options(digits=16)
sqrt(6789)
## [1] 82.39538822045807