# Lifting the Recurse of Factorial

## Thanks for the memoize

Let’s start by revisiting a recursive Python implementation:

`def factorial(n):    return 1 if n == 0 else n * factorial(n - 1)`
`@cachedef factorial(n):    return 1 if n == 0 else n * factorial(n - 1)`
`factorials = [1]def factorial(n):    while len(factorials) <= n:        factorials.append(factorials[-1] * len(factorials))    return factorials[n]`

## Factorial unrolled

We may find that, in practice, our optimisation won’t demand much space — factorial has the kind of growth that makes native integer representations whimper. We can easily get a sense of how rapidly factorial grows:

`>>> for n in range(22): print(f'{n:2}! = {factorial(n)}') 0! = 1 1! = 1 2! = 2 3! = 6 4! = 24 5! = 120 6! = 720 7! = 5040 8! = 40320 9! = 36288010! = 362880011! = 3991680012! = 47900160013! = 622702080014! = 8717829120015! = 130767436800016! = 2092278988800017! = 35568742809600018! = 640237370572800019! = 12164510040883200020! = 243290200817664000021! = 51090942171709440000`
`static const long long factorials[] ={    1,    1,    2,    6,    24,    120,    720,    5040,    40320,    362880,    3628800,    39916800,    479001600,    6227020800,    87178291200,    1307674368000,    20922789888000,    355687428096000,    6402373705728000,    121645100408832000,    2432902008176640000,};long long factorial(int n){    return factorials[n];}`

## Generations

We can precalculate the values by hand or by code and copy them into the table… or just copy them from the code above (go ahead, those numbers are not subject to copyright)… or use code generation:

`print('static const long long factorials[] =\n{')for i in range(21):    print(f'    {factorial(i)},')print('};')`
`#include "factorials.h"long long factorial(int n){    return factorials[n];}`
`print('''\static const long long factorials[] ={''')for i in range(21):    print(f'    {factorial(i)},')print('''\};long long factorial(int n){    return factorials[n];}''')`

## Getting meta

If you’re unhappy with magical constants appearing in your code, no matter what their provenance, another approach is to unroll the calculation for each table entry.

`static const long long factorials[] ={    1L,    1L,    2L*1,    3L*2*1,    4L*3*2*1,    5L*4*3*2*1,    6L*5*4*3*2*1,    7L*6*5*4*3*2*1,    8L*7*6*5*4*3*2*1,    9L*8*7*6*5*4*3*2*1,    10L*9*8*7*6*5*4*3*2*1,    11L*10*9*8*7*6*5*4*3*2*1,    12L*11*10*9*8*7*6*5*4*3*2*1,    13L*12*11*10*9*8*7*6*5*4*3*2*1,    14L*13*12*11*10*9*8*7*6*5*4*3*2*1,    15L*14*13*12*11*10*9*8*7*6*5*4*3*2*1,    16L*15*14*13*12*11*10*9*8*7*6*5*4*3*2*1,    17L*16*15*14*13*12*11*10*9*8*7*6*5*4*3*2*1,    18L*17*16*15*14*13*12*11*10*9*8*7*6*5*4*3*2*1,    19L*18*17*16*15*14*13*12*11*10*9*8*7*6*5*4*3*2*1,    20L*19*18*17*16*15*14*13*12*11*10*9*8*7*6*5*4*3*2*1,};long long factorial(int n){    return factorials[n];}`
`namespace factorials{    constexpr long long calculate(int n)    {        return n == 0 ? 1 : n * calculate(n - 1);    }    constexpr std::array calculated    {        calculate(0),        calculate(1),        calculate(2),        calculate(3),        calculate(4),        calculate(5),        calculate(6),        calculate(7),        calculate(8),        calculate(9),        calculate(10),        calculate(11),        calculate(12),        calculate(13),        calculate(14),        calculate(15),        calculate(16),        calculate(17),        calculate(18),        calculate(19),        calculate(20),    };}constexpr long long factorial(int n){    return factorials::calculated[n];}`

## Transparent maps

Taking a slightly different path, let’s implement factorial in Ada to highlight another language-related consideration. Here’s the iterative version:

`function Factorial(N: Integer) return Long_Integer is    Result: Long_Integer := 1;begin    for I in 1..N loop        Result := Result * Long_Integer(I);    end loop;    return Result;end;`
`function Factorial(N: Integer) return Long_Integer is    (if N = 0 then 1 else Long_Integer(N) * Factorial(N - 1));`
`for N in 0..21 loop    Put_Line(N'Image & "! = " & Factorial(N)'Image);end loop;`
`function Factorial(N: Integer) return Long_Integer is    Factorials: constant array(0..20) of Long_Integer :=    (        1,        1,        2,        6,        24,        120,        720,        5040,        40320,        362880,        3628800,        39916800,        479001600,        6227020800,        87178291200,        1307674368000,        20922789888000,        355687428096000,        6402373705728000,        121645100408832000,        2432902008176640000    );begin    return Factorials(N);end;`
`Factorials: constant array(0..20) of Long_Integer :=(    1,    1,    2,    6,    24,    120,    720,    5040,    40320,    362880,    3628800,    39916800,    479001600,    6227020800,    87178291200,    1307674368000,    20922789888000,    355687428096000,    6402373705728000,    121645100408832000,    2432902008176640000);function Factorial(N: Integer) return Long_Integer is    (Factorials(N));`
`Factorial: constant array(0..20) of Long_Integer :=(    1,    1,    2,    6,    24,    120,    720,    5040,    40320,    362880,    3628800,    39916800,    479001600,    6227020800,    87178291200,    1307674368000,    20922789888000,    355687428096000,    6402373705728000,    121645100408832000,    2432902008176640000);`

## A sterling approximation

So maybe you don’t have a convenient big integer type to hand and are limited to the complimentary two’s complement native integers and floating-point numbers.

`def factorialish(n):    return (n / e)**n * sqrt(2 * pi * n)`
`>>> for n in range(22): print(f'{n:2}! = {factorialish(n):.2f}') 0! = 0.00 1! = 0.92 2! = 1.92 3! = 5.84 4! = 23.51 5! = 118.02 6! = 710.08 7! = 4980.40 8! = 39902.40 9! = 359536.8710! = 3598695.6211! = 39615625.0512! = 475687486.4713! = 6187239475.1914! = 86661001740.6015! = 1300430722199.4716! = 20814114415223.1417! = 353948328666101.1218! = 6372804626194313.0019! = 121112786592294192.0020! = 2422786846761135104.0021! = 50888617325509746688.00`
`def factorialish(n):    return ceil((n / e)**n * sqrt(2 * pi * n))`
`>>> for n in range(22): print(f'{n:2}! = {factorialish(n)}') 0! = 0 1! = 1 2! = 2 3! = 6 4! = 24 5! = 119 6! = 711 7! = 4981 8! = 39903 9! = 35953710! = 359869611! = 3961562612! = 47568748713! = 618723947614! = 8666100174115! = 130043072220016! = 2081411441522417! = 35394832866610218! = 637280462619431319! = 12111278659229419220! = 242278684676113510421! = 50888617325509746688`

## Holding all the cards

You may recall we opened with the observation that shuffling a pack of cards creates something unique in the universe. The probability that two randomly shuffled decks — or two games of 52-card pickup — have the same sequence is 1 in 52!. You now have many ways to calculate how unlikely that is.

--

--

## More from Kevlin Henney

consultant · father · husband · itinerant · programmer · speaker · trainer · writer

## Get the Medium app

consultant · father · husband · itinerant · programmer · speaker · trainer · writer