﻿--- title: The binary splitting with the R gmp package - Application to the Gauss hypergeometric function author: Stéphane Laurent date : 2012-11-30 --- &lead {r setup, echo=FALSE} opts_chunk$set(fig.path="assets/fig/BS2F1-") 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's 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,$\pi = \lim S_nwhere \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} withu_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'_iare 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 alli \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} fori \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) 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); return(a) }  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)