The binary splitting with Haskell

Posted on June 5, 2017 by Stéphane Laurent

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.

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.