--- title: "The binary splitting with the R gmp package - Application to the Gauss hypergeometric function" author: "Stéphane Laurent" date: "2017-06-04" tags: R, special-functions, maths rbloggers: yes output: md_document: variant: markdown toc: true preserve_yaml: true html_document: default highlighter: pandoc-solarized editor_options: chunk_output_type: console --- {r setup, echo=FALSE} options(width=60)  In this article you will firstly see how to get rational numbers arbitrary close to $\pi$ by performing the *binary splitting algorithm* with the gmp package. The *binary splitting algorithm* fastly calculates the partial sums of a rational hypergeometric series by manipulating only integer numbers. But these integer numbers are generally gigantic hence they cannot be handled by ordinary arithmetic computing. After describing the binary splitting algorithm we will show how to implement it in R with the gmp package which allows *arithmetic without limitation*. Our main application is the evaluation of the Gauss hypergeometric function. Introductory example: Euler approximation of $\pi$ ---------------------------------------- The following formula is due to Euler $$\frac{\pi}{2} = 1 + \frac{1}{3} + \frac{1\times 2}{3\times 5} + \frac{1\times 2 \times 3}{3\times 5 \times 7} + \cdots + \frac{n!}{3\times 5 \times 7 \times \cdots \times (2n+1)} + \cdots,$$ that is, $\frac{\pi}{2} = \lim S_n$ where \begin{aligned} S_n & = 1 + \frac{u_1}{v_1} + \frac{u_1 u_2}{v_1v_2} + \frac{u_1u_2 u_3}{v_1v_2v_3} + \cdots + \frac{u_1u_2\ldots u_{n-1}u_n}{v_1v_2\ldots v_{n-1}v_n} \\ & = 1 + \sum_{k=1}^n \prod_{i=1}^k\frac{u_i}{v_i} \\ \end{aligned} with $u_i=i$ and $v_i=2i+1$. Using new notations $\alpha_i = \delta_i = u_i$ and $\beta_i=v_i$ and then writing $$S_n -1 = \frac{\alpha_1}{\beta_1} + \frac{\delta_1 \alpha_2}{\beta_1\beta_2} + \frac{\delta_1\delta_2 \alpha_3}{\beta_1\beta_2\beta_3} + \cdots + \frac{\delta_1\delta_2\ldots\delta_{n-1}\alpha_n}{\beta_1\beta_2\ldots\beta_{n-1}\beta_n}$$ could sound silly at first glance. But now assume $\boxed{n=2^m}$. Then, by summing each $(2i-1)$-st term with the $(2i)$-th term, we can write $S_n-1$ as a sum of $n/2$ terms with a similar expression: $$S_n - 1 = \frac{\alpha'_1}{\beta'_1} + \frac{\delta'_1 \alpha'_2}{\beta'_1\beta'_2} + \frac{\delta'_1\delta'_2 \alpha'_3}{\beta'_1\beta'_2\beta'_3} + \cdots + \frac{\delta'_1\delta'_2\ldots\delta'_{\frac{n}{2}-1}\alpha'_\frac{n}{2}}{\beta'_1\beta'_2\ldots\beta'_{\frac{n}{2}-1}\beta'_{\frac{n}{2}}}$$ where $\alpha'_i$, $\delta'_i$ and $\beta'_i$ are given by \begin{aligned} \alpha'_i = \alpha_{2i-1}\beta_{2_i} + \alpha_{2i}\delta_{2i-1}, \quad \delta'_i = \delta_{2i-1}\delta_{2i} \qquad \text{and } \quad \beta'_i = \beta_{2i-1}\beta_{2i} \end{aligned} for all $i \in \{1, \ldots, n/2\}$. Continuing so on, after $m$ steps we obtain $$S_n - 1 = \frac{\alpha^{(m)}}{\beta^{(m)}}$$ where $\alpha^{(m)}$ and $\beta^{(m)}$ are integer numbers obtained by applying above formulas The above method is the *binary splitting algorithm* for evaluating $S_n$ with $n=2^m$, summarized as follows: 1. Initialization: put $\alpha^{(0)}_i = \delta^{(0)}_i = u_i$ and $\beta^{(0)}_i=v_i$ for $i \in \{1,n\}$; 2. Compute recursively for $k$ going from $1$ to $m$ \begin{aligned} \alpha^{(k)}_i = \alpha^{(k-1)}_{2i-1}\beta^{(k-1)}_{2_i} + \alpha^{(k-1)}_{2i}\delta^{(k-1)}_{2i-1}, \quad \delta^{(k)}_i = \delta^{(k-1)}_{2i-1}\delta^{(k-1)}_{2i} \qquad \text{and } \quad \beta^{(k)}_i = \beta^{(k-1)}_{2i-1}\beta^{(k-1)}_{2i} \end{aligned} for $i \in \{1,n/2^k\}$; 3. Evaluate $S_n = 1 + \frac{\alpha^{(m)}}{\beta^{(m)}}$. The advantage of the binary splitting as compared to a direct evaluation of $S_n$ by summing its $2^m$ terms is twofold: * the binary splitting only performs operations on integer numbers; * it returns an exact expression of $S_n$ as a ratio of two integer numbers. {r} ## example: rational approximation of pi ## bs.pi <- function(m){ u <- function(i) as.numeric(i) v <- function(i) 2*i+1 n <- 2^m indexes <- c(1:n) delta <- alpha <- u(indexes) beta <- v(indexes) j <- 1; l <- n while(jtol){i <- i+1} return(i) }) } irred.frac <- function(x, rnd=n.decimals(x)){ b <- 10^rnd a <- as.bigz(b*round(x,rnd)) num <- a/gcd.bigz(a,b) den <- b/gcd.bigz(a,b) return(list(num=num, den=den)) }  For example: {r} irred.frac(pi) irred.frac(pi, rnd=7)  Finally, here is a user-friendly function for evaluating ${}_2\!F_1$ with the binary splitting: {r} Hypergeometric2F1 <- function(a, b, c, z, m=7, rnd.params=max(n.decimals(c(a,b,c))), rnd.z=n.decimals(z), check.cv=FALSE){ frac.a <- irred.frac(a,rnd.params) frac.b <- irred.frac(b,rnd.params) frac.c <- irred.frac(c,rnd.params) a1 <- frac.a$num; a2 <- frac.a$den b1 <- frac.b$num; b2 <- frac.b$den c1 <- frac.c$num; c2 <- frac.c$den frac.z <- irred.frac(z,rnd.z) p <- frac.z$num; q <- frac.z$den out <- hypergeo_bs(a1,a2, b1,b2, c1,c2, p,q, m) if(check.cv){ x <- hypergeo_bs(a1,a2, b1,b2, c1,c2, p,q, m+1) cv <- x==out out <- list(result=out, convergence=cv) if(!cv){ out$convergence <- paste(out$convergence, " - m=", m, " need to be increased", sep="") } } return(out) }  For example: {r} a <- 20.5; b <- 11.92; c <- 19 z <- 0.5 Hypergeometric2F1(a,b,c,z) Hypergeometric2F1(a,b,c,z, m=3, check.cv=TRUE) Hypergeometric2F1(a,b,c,z, m=7, check.cv=TRUE)  Note that Robin Hankin's gsl package does an excellent job: {r} library(gsl) hyperg_2F1(a,b,c,z)  # Update 2018-11-13 - Converting a bigq rational number to a decimal number with as.numeric is not a good idea. It is better to use the mpfr function of the Rmpfr package, or the q2d function of the package rcdd: {r, message=FALSE} halfpi_bigq <- bs.pi.gmp(8)\$Sn library(Rmpfr) mpfr(halfpi_bigq, 128) library(rcdd) q2d(as.character(halfpi_bigq))  - My function irred.frac is not good. To convert a decimal number to a bigq rational number, it is better to use the d2q function of the rcdd package: {r} as.bigq(d2q(pi))