Monday, February 29, 2016

Factorials and Rationals in R

At the recent NorDevCon, Richard Astbury used the example of factorial to compare some “modern” languages (the theme of his talk being that they were all invented decades ago). Among them was Scheme, which impressed me by having native support for very large numbers and rational numbers.

I felt like my go-to-language for maths-stuff, R, ought to be able to do that, too. The shortest R way to calculate factorials (using built-in functions) I can find is:

prod(1:5)

which gives 120. But it gives me Inf if I try prod(1:50).

BTW, I hoped this might work:

`*`(1:5)

but * operator only takes two arguments (which is ironic for a language where everything is a vector, i.e. an array).

It seems I need to use the gmp package to get big number and rational number support. Here is how factorial can be written:

library(gmp)
prod(as.bigz(1:50))

which outputs
“30414093201713378043612608166064768844377641568960512000000000000”

That is not bad: still fairly short and, as you can see, vectorized operations are still trivial (as.bigz(1:50) creates a vector of 50 gmp numbers).

as.bigq(2, 3) is how you create a rational (⅔). Here is a quick example of creating vectors of rational numbers:

as.bigq(1, 1:50)

which outputs:

Big Rational ('bigq') object of length 50:
[1] 1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
[16] 1/16 1/17 1/18 1/19 1/20 1/21 1/22 1/23 1/24 1/25 1/26 1/27 1/28 1/29 1/30
[31] 1/31 1/32 1/33 1/34 1/35 1/36 1/37 1/38 1/39 1/40 1/41 1/42 1/43 1/44 1/45
[46] 1/46 1/47 1/48 1/49 1/50

Rational operations also work nicely:

as.bigq(1,1:50) + as.bigq(1,3)

giving:
4/3 5/6 2/3 ... 49/138 50/141 17/48 52/147 53/150

Of course that is not as nice as having built-in rationals, but good enough for those relatively few times when you need them.