1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

use rug::Integer;
use std::time::Instant;

fn factorial(input: u64) -> Integer {
    let mut acc = 1.into();

    for i in 1..=input {
        acc *= i;
    }

    acc
}

fn multinomial_coef(k: &[u64]) -> Integer {
    let mut sum = 0;

    for &n in k {
        sum += n;
    }

    let nominator = factorial(sum);

    let mut denominator = Integer::from(1);
    for &n in k {
        denominator *= factorial(n);
    }

    nominator / denominator
}

fn binomial_coef(n: u64, k: u64) -> Integer {
    multinomial_coef(&[n - k, k])
}

pub fn sim(a: u64, b: u64) -> (Integer, Integer) {
    fn enumerate_ras<F: FnMut(Vec<u64>)>(a: u64, b: u64, result: &mut F) {
        fn enumerate_ras_inner<F: FnMut(Vec<u64>)>(
            rem_items: u64,
            rem_sum: u64,
            result: &mut F,
            current: Vec<u64>,
        ) {
            if rem_items == 0 && rem_sum == 0 {
                result(current);
                return;
            }
            if rem_items == 0 {
                return;
            }
            for next in 0..=rem_sum {
                let mut new_current = current.clone();
                new_current.push(next);
                enumerate_ras_inner(rem_items - 1, rem_sum - next, result, new_current);
            }
        }
        enumerate_ras_inner(a, b, result, Vec::new());
    }

    let lhs_single_item = |ras: Vec<u64>| {
        let mut acc = Integer::from(1);
        for r in ras {
            acc *= binomial_coef(b, r);
        }
        acc
    };

    let mut lhs = Integer::new();
    enumerate_ras(a, b, &mut |ras| {
        lhs += lhs_single_item(ras);
    });

    let rhs = binomial_coef(a * b, b);

    (lhs, rhs)
}

pub fn verify() -> (bool, usize) {
    let begin = Instant::now();
    for a in 1..10 {
        for b in 1..10 {
            let (lhs, rhs) = sim(a, b);
            if !(lhs == rhs) {
                return (false, begin.elapsed().as_millis() as usize);
            }
        }
    }
    (true, begin.elapsed().as_millis() as usize)
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

use rug::Integer;
use std::time::Instant;

fn factorial(input: u64) -> Integer {
    (1..=input).product()
}

fn multinomial_coef(k: &[u64]) -> Integer {
    factorial(k.iter().sum::<u64>())
        / k.iter()
            .map(|&n| factorial(n))
            .fold(Integer::from(1), |acc, n| acc * n)
}

fn binomial_coef(n: u64, k: u64) -> Integer {
    multinomial_coef(&[n - k, k])
}

pub fn sim(a: u64, b: u64) -> (Integer, Integer) {
    fn enumerate_ras<F: FnMut(Vec<u64>)>(a: u64, b: u64, result: &mut F) {
        fn enumerate_ras_inner<F: FnMut(Vec<u64>)>(
            rem_items: u64,
            rem_sum: u64,
            result: &mut F,
            current: Vec<u64>,
        ) {
            if rem_items == 0 && rem_sum == 0 {
                result(current);
                return;
            }
            if rem_items == 0 {
                return;
            }
            for next in 0..=rem_sum {
                let mut new_current = current.clone();
                new_current.push(next);
                enumerate_ras_inner(rem_items - 1, rem_sum - next, result, new_current);
            }
        }
        enumerate_ras_inner(a, b, result, Vec::new());
    }

    let lhs_single_item = |ras: Vec<u64>| {
        ras.iter()
            .map(|r| binomial_coef(b, *r))
            .fold(Integer::from(1), |acc, n| acc * n)
    };

    let mut lhs = Integer::new();
    enumerate_ras(a, b, &mut |ras| {
        lhs += lhs_single_item(ras);
    });

    let rhs = binomial_coef(a * b, b);

    (lhs, rhs)
}

pub fn verify() -> (bool, usize) {
    let begin = Instant::now();
    for a in 1..10 {
        for b in 1..10 {
            let (lhs, rhs) = sim(a, b);
            if !(lhs == rhs) {
                return (false, begin.elapsed().as_millis() as usize);
            }
        }
    }
    (true, begin.elapsed().as_millis() as usize)
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

import System.CPUTime (getCPUTime)

-- there is already a product function built in
-- that computes the product of a list of numbers
-- we only need to feed it the sequence [1..n]
-- no need to write any explicit loops or recursion
factorial :: (Integral n) => n -> n
factorial n = product [1 .. n]

mCoef :: (Integral n) => [n] -> n
-- to apply a function to each element of a list
-- we simply `map` the function over the list
-- no need to write any explicit loops or recursion
mCoef xs = factorial (sum xs) `div` product (map factorial xs)

-- this is a 'curried' function
-- which means it takes its arguments one at a time and returns a function
-- that takes the next argument, and so on
-- this is useful for 'partial application' where we supply some of the arguments
-- ahead of time and then supply the rest later, as we will use it below
bCoef :: (Integral n) => n -> n -> n
bCoef n k = mCoef [k, n - k]

sim :: Int -> Int -> (Integer, Integer)
sim a b =
  -- we can define some helper functions inside a function
  let
    -- here we define a function that takes a and b
    -- and generates all possible combinations of
    -- {r_1 ... r_a} such that r_1 + ... + r_a = b
    -- we notice that r_1 ... r_a can be chosen
    -- recursively where we choose r_1 first and then
    -- choose r_2 from the remaining sum, and so on
    -- all until we have chosen a times.
    makeRas :: Int -> Int -> [[Int]]
    -- this is called 'pattern matching'
    -- where the first pattern that matches the input is used
    -- for example the next line states that if a and b are both 0
    -- then return a list containing an empty list
    -- (because there is one way to choose 0 items from 0 items)
    makeRas 0 0 = [[]]
    -- this line states that if a is 0 (and b is not 0 is implied)
    -- then return an empty list (because there is no way to choose
    -- any items from 0 items)
    makeRas 0 _ = []
    -- then we have a catch-all recursive case
    -- that choosing {r_i ... r_a} is equivalent to choosing r_i
    -- first and then choosing {r_{i+1} ... r_a} from the remaining sum
    makeRas remItems remSum =
      -- this is called a 'list comprehension'
      -- it is a way to generate a list by iterating over multiple lists
      -- here we iterate over all possible choices of r_i as 'x'
      -- and then combine it with all possible choices of {r_{i+1} ... r_a}
      [x : xs | x <- [0 .. remSum], xs <- makeRas (remItems - 1) (remSum - x)]

    -- this computes a single summation term on the left hand side
    -- given a list of r values
    lhsSingleItem :: [Int] -> Integer
    -- this is a bit complex and may be confusing if you are not familiar with functional programming
    -- firstly 'fromIntegral' is simply used to convert between machine and arbitrary precision integers
    -- notice here we used the 'partial application' of 'bCoef'
    -- we first supplied the first argument to be 'b', and then
    -- we pass the partially applied function to 'map' to compute
    -- [bCoef b r_1, bCoef b r_2, ...]
    -- then we use 'product' to compute the product of all these values
    lhsSingleItem ras = product $ map (bCoef $ fromIntegral b) (map fromIntegral ras)

    -- compute the summation by summing the mapped values of 'lhsSingleItem' over all possible 'ras'
    lhs = sum $ map lhsSingleItem $ makeRas a b
    -- right hand side is trivial to compute, we just give it a name for consistency
    rhs = bCoef (fromIntegral (a * b)) (fromIntegral b)
   in
    -- this function return the pair (lhs, rhs)
    (lhs, rhs)

main :: IO ()
-- the 'do' block would be the closest resemblance to imperative programming
-- it actually strings operations together sequentially, like statements in an imperative language
-- (side note: even the 'do' notation is technically still just a composition of functions, it is just syntactic sugar)
main = do
  begin <- getCPUTime
  putStr
    $ show
    -- all is a function that takes a list of boolean values and returns whether
    -- all of them are True
    $ all
      -- uncurry as its name suggests convert a curried function to a function
      -- that takes a single pair of arguments
      -- notice that operators in Haskell are functions too
      -- so we can use 'uncurry' on the equality operator to
      -- get back a function that checks whether a pair of values are equal
      (uncurry (==))
    -- we use list comprehension again to generate all pairs of numbers from 1 to 10
    $ map (uncurry sim) [(x, y) | x <- [1 .. 10], y <- [1 .. 10]]
  end <- getCPUTime
  putStrLn $ "," ++ show (fromIntegral (end - begin) / 10 ^ 9)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61

import java.math.BigInteger

fun factorial(n: Long) : BigInteger {
    return (1..n).fold(BigInteger.ONE) 
        { acc, i -> acc * BigInteger.valueOf(i.toLong()) }
}

fun mCoef(ks: Array<Long>): BigInteger {
    val n = ks.sum()
    val num = factorial(n)
    val den = ks.fold(BigInteger.ONE) { acc, ki -> acc * factorial(ki) }
    return num / den
}

fun bCoef(n: Long, k: Long): BigInteger {
    return mCoef(arrayOf(n - k, k))
}

fun makeRas(a: Long, b:Long): Sequence<Array<Long>> {
    fun makeRasInner(rem_items: Long, rem_sum: Long, current: Array<Long>) : Sequence<Array<Long>> {
        return sequence {
            if (rem_items == 0L) {
                if (rem_sum == 0L) {
                    yield(current)
                }
            } else {
                for (i in 0..rem_sum) {
                    val new_current = current.toMutableList()
                    new_current.add(i)
                    yieldAll(makeRasInner(rem_items - 1, rem_sum - i, new_current.toTypedArray()))
                }
            }
        }
    }

    return makeRasInner(a, b, arrayOf())
}

fun sim(a: Long, b:Long): Pair<BigInteger, BigInteger> {
    fun lhsSingleItem(ras: Array<Long>) : BigInteger {
        return ras.fold(BigInteger.ONE) { acc, ki -> acc * bCoef(b, ki) }
    }
    
    val lhs = makeRas(a, b).map { lhsSingleItem(it) }.fold(BigInteger.ZERO) { acc, i -> acc + i }
    val rhs = bCoef(a * b, b)
    return Pair(lhs, rhs)
}

fun main() {
    val start = System.currentTimeMillis()
    for (a in 1..10) {
        for (b in 1..10) {
            val (lhs, rhs) = sim(a.toLong(), b.toLong())
            if (lhs != rhs) {
                throw Exception("Error: $a, $b, $lhs, $rhs")
            }
        }
    }
    val end = System.currentTimeMillis()
    println("true,${end - start}")
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

lines <- read_csv("lines.csv")
timings <- read_csv("timings.csv")

p1 <- lines %>%
  ggplot(aes(x = paradigm, y = lines, fill = paradigm)) +
  geom_col() +
  geom_text(aes(label = lines), vjust = -0.5) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1), add = c(0, 10))) +
  labs(title = "Lines of Code in Different Paradigms",
       x = "Paradigm",
       y = "Lines of Code")

p2 <- timings %>%
  ggplot(aes(x = paradigm, y = time_ms, fill = paradigm)) +
  geom_boxplot() +
  labs(title = "Timing of Different Paradigms",
       x = "Paradigm",
       y = "Time (ms)")
plot_grid(p1, p2, ncol = 1)