---
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_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(j<n){
       l <- l/2
       odd <- 2*c(1:l); even <- odd-1
       alpha <- beta[odd]*alpha[even] + delta[even]*alpha[odd]
       j <- 2*j
       beta <- beta[odd]*beta[even]
       delta <- delta[even]*delta[odd]
       }
     Sn <- alpha/beta + 1
     out <- list(alpha=alpha, beta=beta, Sn=Sn)
     return(out)
}

```

The method very well performs while $m\leq 7$ :
```{r}
print(bs.pi(7),digits=22)
print(pi/2,digits=22)
```

But the numerator and the denominator become too gigantic when $m=8$:
```{r}
bs.pi(8)
```




Second example: exponential  of a rational number
----------------------------------------------------------

It is well known that $\exp(x)=\lim S_n(x)$ where $S_n(x)=\sum_{k=0}^n\frac{x^n}{n!}$. 
Thus, when $x=p/q$ for some integers $p$ and $q$, we can write as before
$$
S_n(x) = 1 + \sum_{k=1}^n \prod_{i=1}^k\frac{u_i}{v_i}
$$
where $u_i \equiv p$ and $v_i= i q$ are integer numbers. Thus, we can use the binary splitting algorithm to compute $S_{2^m}$:
```{r}
## example: rational approximation of exp(p/q) ##
bs.exp <- function(p,q,m){
    v <- function(i) i*q
     n <- 2^m
    indexes <- 1:n
     delta <- alpha <- rep(p,n)
    beta <- v(indexes)
     j <- 1; l <- n
     while(j<n){
       l <- l/2
       odd <- 2*c(1:l); even <- odd-1
       alpha <- beta[odd]*alpha[even] + delta[even]*alpha[odd]
       j <- 2*j
       beta <- beta[odd]*beta[even]
       delta <- delta[even]*delta[odd]
       }
     Sn <- alpha/beta + 1
     out <- list(alpha=alpha, beta=beta, Sn=Sn)
     return(out)
}

```

Let us try to evaluate $\exp(1)$. For $m=7$, the approximation is not entirely satisfactory:
```{r}
print(bs.exp(1,1,7), digits=22)
print(exp(1), digits=22)
```
And for $m=8$, it crashes:
```{r}
bs.exp(1,1,8)
```



The `gmp` package comes to our rescue
----------------------------------------

As we noted above, the binary splitting manipulates only *integer* numbers. 
The evaluation of $\exp(1)$ has crashed because the numerator and the denominator were too big integers. 
The crantastic [`gmp`](http://www.inside-r.org/packages/cran/gmp) package overcomes this problem because it allows ``arithmetic without limitations'' using the [C library GMP (GNU Multiple Precision Arithmetic)](http://gmplib.org/).

Let us show how the `gmp` works on the $\pi$ example. This is very easy: we only have to convert the two input sequences of integers $(u_i)$ and $(v_i)$ to sequences of `bigz` integers:
  
```{r, message=FALSE}
library(gmp)
## rational approximation of pi with gmp ##
bs.pi.gmp <- function(m){
    u <- function(i) as.numeric(i)
    v <- function(i) 2*i+1
     n <- 2^m
    indexes <- 1:n
     delta <- alpha <- as.bigz(u(indexes))
     beta <- as.bigz(v(indexes))
     j <- 1; l <- n
     while(j<n){
       l <- l/2
       odd <- 2*c(1:l); even <- odd-1
       alpha <- beta[odd]*alpha[even] + delta[even]*alpha[odd]
       j <- 2*j
       beta <- beta[odd]*beta[even]
       delta <- delta[even]*delta[odd]
       }
     Sn <- alpha/beta + 1
  out <- list(Sn=Sn, eval.Sn=format(as.numeric(Sn),digits=22))
return(out)
}

```

The evaluation of $S_n$ with $n=2^3$ illustrates the first advantage of the `gmp` package:

```{r}
bs.pi.gmp(3)
bs.pi(3)
```
As you can see, $S_n$ is written as an irreducible fraction with the `gmp` approach. 
But this is not the main strength of the  `gmp` package. Now we have (almost) no limitation on $m$ for evaluating $S_{2^m}$:
```{r, tidy=FALSE}
bs.pi.gmp(8)
```

Obviously the first limitation is the width of your screen. The more serious limitations of the `gmp` package are beyond the scope of this article. 

Let us come back to the exponential example:
```{r}
## rational approximation of exp(p/q) with gmp ##
bs.exp.gmp <- function(p,q,m){
      v <- function(i) i*q
     n <- 2^m
    indexes <- 1:n
     delta <- alpha <- as.bigz(rep(p,n))
    beta <- as.bigz(v(indexes))
     j <- 1; l <- n
     while(j<n){
       l <- l/2
       odd <- 2*c(1:l); even <- odd-1
       alpha <- beta[odd]*alpha[even] + delta[even]*alpha[odd]
       j <- 2*j
       beta <- beta[odd]*beta[even]
       delta <- delta[even]*delta[odd]
       }
     Sn <- alpha/beta + 1
  out <- list(Sn=Sn, eval.Sn=format(as.numeric(Sn),digits=22))
return(out)
}

```

```{r}
 bs.exp.gmp(1,1,8)
```

Very well.



A general function for the binary splitting algorithm 
----------------------------------------------------------------

Before turning to the Gauss hypergeometric function we write a general function for the binary splitting taking as arguments the two sequences $(u_i)$ and $(v_i)$:
```{r}
bs.gmp <- function(u,v,m=7,value="eval"){
     n <- 2^m
    indexes <- 1:n
     delta <- alpha <- as.bigz(u(indexes))
    beta <- as.bigz(v(indexes))
     j <- 1; l <- n
     while(j<n){
       l <- l/2
       odd <- 2*c(1:l); even <- odd-1
       alpha <- beta[odd]*alpha[even] + delta[even]*alpha[odd]
       j <- 2*j
       beta <- beta[odd]*beta[even]
       delta <- delta[even]*delta[odd]
       }
     Sn <- alpha/beta + 1
     eval.Sn <- format(as.numeric(Sn) ,digits=22)
     out <- switch(value, "eval"=eval.Sn, "exact"=Sn, "both"=list(Sn=Sn, eval.Sn=eval.Sn))
return(out)
}

```






The Gauss hypergeometric function 
-----------------------------------------------------------------


Now consider the *Gauss hypergeometric function*  ${}_2\!F_1$. 
This is the function ${}_2\!F_1(\alpha,\beta,\gamma; \cdot)$  
with complex parameters $\alpha$, $\beta$,  $\gamma \not\in \mathbb{Z}^-$ and complex variable $z$  defined for $|z|<1$ as the sum of an absolute convergent series
$${}_2\!F_1(\alpha,\beta,\gamma; z) = \sum_{n=0}^{\infty}\frac{ {(\alpha)}_{n}{(\beta)}_n}{ {(\gamma)}_{n}}\frac{z^n}{n!},$$
and extended by analytical continuation in the complex plane with the cut 
along $(1,+\infty)$. Here ${(a)}_n:=a(a+1)\cdots(a+n-1)$ denotes Pochhammer's symbol used 
to represent the $n$-th ascending factorial of $a$. . 

The binary splitting allows to evaluate ${}_2\!F_1(\alpha,\beta,\gamma; z)$ for rational values of 
$\alpha,\beta,\gamma, z$ by manipulating only integer numbers. 
This is performed by the R function below

```{r}
## rational approximation of  2F1(a1/a2, b1/b2, c1/c2; p/q) with gmp ##
hypergeo_bs <- function(a1,a2, b1,b2, c1,c2, p,q, m){
  u <- function(i) c2*(a1+(i-1)*a2)*(b1+(i-1)*b2)*p
  v <- function(i) a2*b2*i*(c1+(i-1)*c2)*q
	bs.gmp(u,v,m)
}

```


For more convenience I have firstly written the function below which returns the irreducible rational 
notation of a given number $x$. The user can also specify a rounding order for $x$. 
```{r}
n.decimals <- function(x, tol=.Machine$double.eps){
  sapply(x, function(x) {
		i <- 0
		while(abs(x-round(x,i))>tol){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)

```