--- author: Stéphane Laurent date: '2017-06-05' highlighter: hscolour output: html_document: default md_document: variant: markdown prettify: True tags: 'haskell, special-functions, maths' title: The binary splitting with Haskell --- At the first line of each script in this article, we'll load the following small Haskell module:

-- BinarySplitting.hs
module BinarySplitting
  where
import Data.Ratio ((%))
 
split0 :: ([Rational], [Rational]) -> [Rational]
split0 u_v = map (\i -> (u !! (2*i)) * (v !! (2*i+1))) [0 .. m]
  where (u, v) = u_v
        m = div (length u) 2 - 1
 
split1 :: ([Rational], [Rational], [Rational]) ->
               ([Rational], [Rational], [Rational])
split1 adb = split adb (length alpha)
  where (alpha, _, _) = adb
        split :: ([Rational], [Rational], [Rational]) -> Int ->
                             ([Rational], [Rational], [Rational])
        split u_v_w n =
          if n == 1
            then u_v_w
            else split (x, split0 (v,v), split0 (w,w)) (div n 2)
          where (u, v, w) = u_v_w
                x  = zipWith (+) (split0 (u, w)) (split0 (v, u))
 
bsplitting :: Int -> [Rational] -> [Rational] -> Rational
bsplitting m u v = num / den + 1
  where ([num], _, [den]) = split1 (take (2^m) u, take (2^m) u, take (2^m) v)

The `bsplitting` function performs the [binary splitting algorithm](https://laustep.github.io/stlahblog/posts/hypergeometric.html). Given an integer $m \geq 0$ and two sequences $(u_i)$ and $(v_i)$ of rational numbers, it calculates the sum $$ A_m = 1 + \sum_{k=1}^{2^m} \prod_{i=1}^k\frac{u_i}{v_i}. $$ Approximation of $\pi$ ---------------------- For example, $A_m \to \frac{\pi}{2}$ when $u_i = i$ and $v_i = 2i+1$. So we get a rational approximate of $\pi$ as follows:

> :load BinarySplitting.hs
> 
> :{
> approxPi :: Int -> Rational
> approxPi m = 2 * bsplitting m u v
>   where u = [i | i <- [1 ..]]
>         v = [2*i+1 | i <- [1 ..]]
> :}
> 
> let x = approxPi 5
> x
12774464002301303455744 % 4066238182722121490175
> fromRational x
3.141592653519747

Kummer hypergeometric function ------------------------------ Consider the confluent hypergeometric series $$ {}_1\!F_1(a, b; x) = \sum_{n=0}^{\infty}\frac{{(a)}_{n}}{{(b)}_{n}}\frac{x^n}{n!} = 1 + \sum_{n=1}^{\infty}\frac{{(a)}_{n}}{{(b)}_{n}}\frac{x^n}{n!}. $$ Here ${(a)}_n:=a(a+1)\cdots(a+n-1)$ is the Pocchammer symbol denoting the ascending factorial. The sum from $n=0$ to $n=2^m$ is evaluated by the `bsplitting` function by taking the sequences\ $u_i = (a+i-1)x$ and $v_i = (b+i-1)i$. Below we evaluate this sum for $a=8.1$, $b=10.1$ and $x=100$. We compare the result to the value of ${}_1\!F_1(a, b; x)$ given by Wolfram.

> :load BinarySplitting.hs
> 
> :{
> hypergeo1F1 :: Int -> Rational -> Rational -> Rational -> Double
> hypergeo1F1 m a b x = fromRational $ bsplitting m u v
>   where u = [(a+i)*x | i <- [0 ..]]
>         v = [(b+i)*(i+1) | i <- [0 ..]]
> :}
> 
> let wolfram = 1.7241310759926883216143646e41
> 
> wolfram - hypergeo1F1 6 (81%10) (101%10) 100
1.7238238908740056e41
> wolfram - hypergeo1F1 7 (81%10) (101%10) 100
3.0481841873624932e38
> wolfram - hypergeo1F1 8 (81%10) (101%10) 100
0.0

We find a good approximate for $m=8$ (so $2^m=256$), and the evaluation is really fast.