agm: Arithmetic-geometric Mean

Description Usage Arguments Details Value Note References See Also Examples

View source: R/agm.R

Description

The arithmetic-geometric mean of real or complex numbers.

Usage

1
agm(a, b)

Arguments

a, b

real or complex numbers.

Details

The arithmetic-geometric mean is defined as the common limit of the two sequences a_{n+1} = (a_n + b_n)/2 and b_{n+1} = √(a_n b_n).

Value

Returnes one value, the algebraic-geometric mean.

Note

The calculation of the AGM is continued until the two values of a and b are identical (in machine accuracy).

References

Borwein, J. M., and P. B. Borwein (1998). Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity. Second, reprinted Edition, A Wiley-interscience publication.

See Also

Arithmetic, geometric, and harmonic mean.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
##  Gauss constant
1 / agm(1, sqrt(2))  # 0.834626841674073

##  Calculate the (elliptic) integral 2/pi \int_0^1 dt / sqrt(1 - t^4)
f <- function(t) 1 / sqrt(1-t^4)
2 / pi * integrate(f, 0, 1)$value
1 / agm(1, sqrt(2))

##  Calculate pi with quadratic convergence (modified AGM)
#   See algorithm 2.1 in Borwein and Borwein
y <- sqrt(sqrt(2))
x <- (y+1/y)/2
p <- 2+sqrt(2)
for (i in 1:6){
  cat(format(p, digits=16), "\n")
  p <- p * (1+x) / (1+y)
  s <- sqrt(x)
  y <- (y*s + 1/s) / (1+y)
  x <- (s+1/s)/2
  }

## Not run: 
##  Calculate pi with arbitrary precision using the Rmpfr package
require("Rmpfr")
vpa <- function(., d = 32) mpfr(., precBits = 4*d)
# Function to compute \pi to d decimal digits accuracy, based on the 
# algebraic-geometric mean, correct digits are doubled in each step.
agm_pi <- function(d) {
    a <- vpa(1, d)
    b <- 1/sqrt(vpa(2, d))
    s <- 1/vpa(4, d)
    p <- 1
    n <- ceiling(log2(d));
    for (k in 1:n) {
        c <- (a+b)/2
        b <- sqrt(a*b)
        s <- s - p * (c-a)^2
        p <- 2 * p
        a <- c
    }
    return(a^2/s)
}
d <- 64
pia <- agm_pi(d)
print(pia, digits = d)
# 3.141592653589793238462643383279502884197169399375105820974944592
# 3.1415926535897932384626433832795028841971693993751058209749445923 exact

## End(Not run)

Example output

[1] 0.8346268
[1] 0.8346268
[1] 0.8346268
3.414213562373095 
3.142606753941623 
3.141592660966044 
3.141592653589793 
3.141592653589793 
3.141592653589793 
Loading required package: Rmpfr
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following objects are masked from 'package:stats':

    dbinom, dnorm, dpois, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

1 'mpfr' number of precision  256   bits 
[1] 3.141592653589793238462643383279502884197169399375105820974944592

numbers documentation built on May 15, 2021, 1:08 a.m.