Deranged !Combinatorics: Compute the Subfactorial

The subfactorial or rencontres numbers (A000166) are a sequence of numbers similar to the factorial numbers which show up in the combinatorics of permutations. In particular the nth subfactorial !n gives the number of derangements of a set of n elements. A derangement is a permutation in which no element remains in the same position. The subfactorial can be defined via the following recurrence relation:

!n = (n-1) (!(n-1) + !(n-2))


In fact, the same recurrence relation holds for the factorial, but for the subfactorial we start from:

!0 = 1
!1 = 0


(For the factorial we'd have, of course, 1! = 1.)

Rules

Like the factorial, the subfactorial grows very quickly. It is fine if your program can only handle inputs n such that !n can be represented by your language's native number type. However, your algorithm must in theory work for arbitrary n. That means, you may assume that integral results and intermediate value can be represented exactly by your language. Note that this excludes the constant e if it is stored or computed with finite precision.

The result needs to be an exact integer (in particular, you cannot approximate the result with scientific notation).

You may write a program or a function and use any of the standard methods of receiving input and providing output.

You may use any programming language, but note that these loopholes are forbidden by default.

This is , so the shortest valid answer – measured in bytes – wins.

Test Cases

n     !n
0     1
1     0
2     1
3     2
4     9
5     44
6     265
10    1334961
12    176214841
13    2290792932
14    32071101049
20    895014631192902121
21    18795307255050944540
100   34332795984163804765195977526776142032365783805375784983543400282685180793327632432791396429850988990237345920155783984828001486412574060553756854137069878601


Funciton, 336 bytes

Byte count assumes UTF-16 encoding with BOM.

┌─╖┌─╖  ┌─╖
│f╟┤♭╟┐┌┤♭╟┐
╘╤╝╘═╝├┘╘═╝├────┐
│┌─╖ │ ┌┐┌┘╔═╗╓┴╖
││f╟─┴┐└┴┼─╢0║║f║
│╘╤╝  │  │ ╚═╝╙─╜
│┌┴╖ ┌┴╖┌┴╖ ╔═╗
││+╟┐│×╟┤?╟┐║1║
│╘╤╝│╘╤╝╘╤╝┘╚╤╝
└─┘ └─┘  └───┘


This defines a function f which takes one integer and outputs another integer at a 90 degree turn to the left. It works for arbitrarily large inputs.

Try it online!

Considering this is Funciton it's even reasonably fast (n = 20 takes about 14 seconds on TIO). The main slowdown comes from the double-recursion, as I don't think the Funciton interpreter automatically memoises functions.

Unfortunately, some monospaced fonts don't correctly monospace the ♭ and/or insert small gaps between the lines. Here's a screenshot of the code from TIO in all its beauty:

I think it might be possible to golf this some more, e.g. by changing the condition from >0 to <1 and swapping the branches of the conditional, so that I could reuse the number literal, or maybe by using a completely different formula, but I'm quite happy with how compact it is already.

Explanation

This basically implements the recursion given in the challenge, although it uses the base case !(-1) = !0 = 1. n-1 and n-2 are computed with the predecessor function ♭, and the intermediate result n-1 is reused in three places. There isn't much more to it, so I'll just quickly go through the control flow:

               ─┐
╓┴╖
║f║
╙─╜


This is the function header which emits the function's input n long the attached line. This immediately reaches T-junction, which simply duplicates the value.

        ┌┐┌┘╔═╗
└┴┼─╢0║
│ ╚═╝


The 0 box is just a numeric literal. A 4-way junction computes two functions: the path that leads off the bottom computes 0 < n, which we'll use to determine the base case. The path that goes left separately computes 0 << n (a left-shift), but we discard this value with the Starkov construct.

         ┌┴╖ ╔═╗
┤?╟┐║1║
╘╤╝┘╚╤╝
└───┘


We lead this into the three-way conditional ?. If the value is false, we return the constant result 1. The loose end right of the ? is the function output. I'm twisting it around by 180 degrees here, so that the relative orientation of input and output of f is more convenient in the rest of the program.

If the condition was true, then the other value will be used. Let's look at the path that leads to this branch. (Note that Funciton's evaluation is actually lazy so that this branch will never be evaluated if it's not needed, which makes the recursion possible in the first place.)

        ┌─╖
┐┌┤♭╟┐
├┘╘═╝
│
─┴┐


In the other branch we first compute n-1 and then split the path twice so we get three copies of the value (one for the coefficient of the recurrence, one for the first subfactorial, the last for n-2).

┌─╖┌─╖
│f╟┤♭╟
╘╤╝╘═╝
│┌─╖
││f╟
│╘╤╝
│┌┴╖
││+╟
│╘╤╝
└─┘


Like I said, we decrement one copy again with another ♭, then we feed both n-1 and n-2 recursively to f and finally add the two results together in the +.

       ┐
│
┌┴╖
┐│×╟
│╘╤╝
└─┘


All that's left is to multiply n-1 by !(n-1) + !(n-2).

Oasis, 5 bytes

Uses the formula given by Martin. Code:

+n<*X


Dissected version:

+n<*


with a(0) = 1 and a(1) = 0.

Explanation, a(n) =:

+       # Add the previous two terms, a(n - 1) + a(n - 2).
n<     # Compute n - 1.
*    # Multiply the top two elements.


Try it online!

• Nice trick using X :-) BTW, did you implement this yet? One of these days we won't be able to get away with just changing the initial values – Luis Mendo Mar 26 '17 at 17:18
• @LuisMendo Yes I did! It's used as a command flag (here is a link to the info page). Thank you for the suggestion :). – Adnan Mar 26 '17 at 17:24

f 0=1
f n=n*f(n-1)+(-1)^n


Try it online! Uses the other recurrence from the OEIS page.

Jelly, 7 bytes

R=Œ!Ḅċ0


This approach constructs the derangements, so it's rather slow.

Try it online!

How it works

R=Œ!Ḅċ0  Main link. Argument: n

R        Range; yield [1, ..., n].
Œ!     Yield all permutations of [1, ..., n].
=       Perform elementwise comparison of [1, ..., n] and each permutation.
Ḅ    Unbinary; convert each result from base 2 to integer. This yields 0 for
derangements, a positive value otherwise.
ċ0  Count the number of zeroes.


Brachylog (2), 11 bytes

⟦₁{p:?\≠ᵐ}ᶜ


Try it online!

Explanation

This is basically just a direct translation of the spec from English to Brachylog (and thus has the advantage that it could easily be modified to handle small changes to the specification, such as finding the number of derangements of a specific list).

⟦₁{p:?\≠ᵐ}ᶜ
{      }ᶜ  Then count the number of ways to
p         permute that list
\      such that taking corresponding elements
:?       in {the permutation} and the list of distinct elements
≠     gives different elements
ᵐ    at every position


Languages with built-in solutions

Following xnor's suggestion this is a CW answer into which trivial solutions based on a single built-in to compute the subfactorial or generate all derangements should be edited.

Mathematica, 12 bytes

Subfactorial

• sigh Mathematica... – epicbob57 Mar 28 '17 at 2:41

Python 3, 35 32 bytes

f=lambda n:n<1or(-1)**n+n*f(n-1)


This uses the recurrence relation !n = n !(n-1) + (-1)n from @Laikoni's Haskell answer, with base case !0 = 1.

Try it online!

• I think you can also use the other equation given here, which would save two bytes: f=lambda n:n<1or n*f(n-1)+(-1)**n. – Adnan Mar 26 '17 at 17:21
• Three bytes with a little reordering. ;) – Dennis Mar 26 '17 at 17:25
• The fun part of this recurrence is that if you push the base case back to n=-1, it doesn't matter at all which value you use. That could be useful for some languages (e.g. in Mathematica you could actually leave it undefined if that saved any bytes). – Martin Ender Mar 26 '17 at 22:07

M, 9 bytes

o2!÷Øe+.Ḟ


As you can see by removing the Ḟ, M uses symbolic math, so there will be no precision issues.

Try it online! Not the shortest solution that has been posted, but fast.

How it works

o2!÷Øe+.Ḟ  Main link. Argument: n

o2         Replace input 0 with 2, as the following formula fails for 0.
!        Compute the factorial of n or 2.
÷Øe     Divide the result by e, Euler's natural number.
+.   Add 1/2 to the result.
Ḟ  Floor; round down to the nearest integer.


MATL, 9 8 bytes

:tY@-!As


Similarly to @Dennis' Jelly answer, this actually builds the permutations and counts how many of them are derangements; so it is slow.

Try it online!

:     % Input n implicitly: Push [1 2 ... n]
t     % Duplicate
Y@    % Matrix of all permutations, each on a row
-     % Element-wise subtract. A zero in a row means that row is not a derangement
!     % Transpose
A     % True for columns that don't contain zeros
s     % Sum. Implicitly display


Mathics, 21 bytes

Round@If[#>0,#!/E,1]&


I am very new to this and have no idea what I'm doing...

Try it online!

• Two alternatives at the same byte count: Round[(#/. 0->2)!/E]& and ±0=1;±n_:=Round[n!/E] (although I don't know whether Mathics supports single-byte encodings for source files like Mathematica does). – Martin Ender Mar 26 '17 at 18:09
• The first one works well (I think I know what it does), but Mathics doesn't appear to support ± in the second one. It would work with f, but at the cost of two bytes. – Dennis Mar 26 '17 at 18:23
• Another at the same byte count: Round[#!/E]+1-Sign@#&. Annoying initial values...! – Greg Martin Mar 26 '17 at 19:48

Ruby, 27 bytes

f=->n{n<1?1:n*f[n-1]+~0**n}


~0**n is shorter than (-1)**n!

CJam (10 bytes)

1qi{~*)}/z


This uses the recurrence !n = n !(n-1) + (-1)^n, which I derived from n! / e and then discovered was already in OEIS.

Dissection

The loop calculates (-1)^n !n, so we need to take the absolute value at the end:

1     e# Push !0 to the stack
qi{   e# Read an integer n and loop from 0 to n-1
~   e#   Bitwise not takes i to -(i+1), so we can effectively loop from 1 to n
*   e#   Multiply
)   e#   Increment
}/
z     e# Take the absolute value


05AB1E, 8 bytes

ÎƒN*®Nm+


Try it online!

Explanation

Î         # initialize stack with 0 and input
ƒ        # for N in range [0 ... input]:
N*      # multiply top of stack with N
®Nm   # push (-1)^N


MATLAB, 33 bytes

@(n)(-1)^n*hypergeom([1 -n],[],1)


Anonympus function that uses the formula in Section 3 of Derangements and applications by Mehdi Hassani.

Example use:

>> @(n)(-1)^n*hypergeom([1 -n],[],1)
ans =
@(n)(-1)^n*hypergeom([1,-n],[],1)
>> ans(6)
ans =
265


JavaScript (ES6), 26 bytes

f=n=>!n||n*f(n-1)-(~n%2|1)


Uses the recurrence relation from @Laikoni's answer. In ES7 you can save a byte by using +(-1)**n instead of -(~n%2|1).

PostScript, 8176 69 bytes

Here are implementations of both formulae.

n*f(n-1)+(-1)^n

/f{dup 0 eq{pop 1}{dup dup 1 sub f mul exch 2 mod 2 mul 1 sub sub}ifelse}def

/f{dup 0 eq{pop 1}{dup dup 1 sub f mul -1 3 2 roll exp add}ifelse}def


That version outputs a float. If it's necessary to output an integer:

/f{dup 0 eq{pop 1}{dup dup 1 sub f mul -1 3 2 roll exp cvi add}ifelse}def


which weighs in at 73 bytes.

The other formula is a little longer: 81 bytes.

(n-1)*(f(n-1)+f(n-2))

/f{dup 1 le{1 exch sub}{1 sub dup f exch dup 1 sub f 3 -1 roll add mul}ifelse}def


These functions get their argument from the stack, and leave the result on the stack.

You can test the functions, either in a file or at an interactive PostScript prompt (eg GhostScript) with

0 1 12{/i exch def [i i f] ==}for


output

[0 1]
[1 0.0]
[2 1.0]
[3 2.0]
[4 9.0]
[5 44.0]
[6 265.0]
[7 1854.0]
[8 14833.0]
[9 133496.0]
[10 1334961.0]
[11 14684570.0]
[12 176214848.0]


Here's a complete PostScript file which renders the output to the screen or a printer page. (Comments in PostScript start with %).

%!PS-Adobe-3.0

% (n-1)*(f(n-1)+f(n-2))
% /f{dup 1 le{1 exch sub}{1 sub dup f exch dup 1 sub f 3 -1 roll add mul}ifelse}def

% n*f(n-1)+(-1)^n
/f{dup 0 eq{pop 1}{dup dup 1 sub f mul -1 3 2 roll exp add}ifelse}def

% 0 1 12{/i exch def [i i f] ==}for

/FS 16 def              %font size
/LM 5 def               %left margin
/numst 12 string def    %numeric string buffer

/Newline{currentpoint exch pop FS sub LM exch moveto}def
/Courier findfont FS scalefont setfont
LM 700 moveto

(Subfactorials) Newline
0 1 12{
dup numst cvs show (: ) show f numst cvs show Newline
}for
showpage
quit

• -1 3 2 roll exp is a fair bit shorter than exch 2 mod 2 mul 1 sub. – Peter Taylor Mar 27 '17 at 9:52
• @PeterTaylor So it is! :) I forgot about exp :oops: However, it returns a float, and I think I need to output an integer to conform to the question. – PM 2Ring Mar 27 '17 at 10:30
• I think you can still chuck in a cvi and make a saving. (Note: untested, but from reading the doc I think it should work). – Peter Taylor Mar 27 '17 at 10:38
• @PeterTaylor Yes, cvi works, and it's still shorter than my previous version. – PM 2Ring Mar 27 '17 at 11:04

function f($i){return$i>1?$i*f($i-1)+(-1)**$i:1-$i;}echo f($argv[1]);  use this way a(n) = n*a(n-1) + (-1)^n • You only need to give the function, not the full program, so you can drop the last 17 chars. There's a further saving by not special-casing input 1. I think the two savings take it down to 47 bytes. – Peter Taylor Mar 28 '17 at 10:49 PHP, 50 44 for(;$i++<$argn;)$n=++$n*$i-$i%2*2;echo$n+1;


Run with echo <n> | php -nR '<code>

The beauty of a(n) = n*a(n-1) + (-1)^n is that it depends only on the previous value. This allows it to be implemented iteratively instead of recursively. This saves a long function declarition.

-6 bytes by @Titus. Thanks !

• -1 byte: $i++<$argv[1]. -2 bytes: for(;$i++<$argv[1];)$n=++$n*$i-$i%2*2;echo$n+1;. (-3 bytes with -R and $argn.) – Titus Mar 30 '17 at 8:36
• @Titus someone got bored ? :D would you mind giving me an example for -R and $argn? – Christoph Mar 30 '17 at 10:42 • Not bored, but eager to golf. See php.net/manual/de/features.commandline.options.php: echo <input> | php -nR '<code>'. example: codegolf.stackexchange.com/a/113046 – Titus Mar 30 '17 at 10:56 • @Titus did I get it right ? ;-) – Christoph Mar 30 '17 at 11:03 Mathematica, 40 bytes ±0=1;±1=0;±n_:=(n-1)(±(n-1)+±(n-2))  Which comes in at 31 bytes under the default ISO 8859-1 encoding. C, 34 bytes a(n){return n?n*a(n-1)-n%2*2+1:1;}  Explanation: a(n){ } define a function called a of n return ; make the function evaluate to... n? :1 set the base case of 1 when n is 0 n*a(n-1) first half of the formula on the page -n%2*2+1 (-1)**n  R, 47 bytes n=scan();if(!n,1,floor(gamma(n+1)/exp(1)+.5))  Uses the same formula as Mego's answer. Alternate method, 52 bytes using the PerMallows library n=scan();if(!n,1,PerMallows::count.perms(n,n,'h'))  Actually, 18 bytes ;!@ur⌠;!@0Dⁿ/⌡MΣ*≈  Try it online! Explanation: ;!@ur⌠;!@0Dⁿ/⌡MΣ*≈ ; duplicate input ! n! @ur range(0, n+1) (yields [0, n]) ⌠;!@0Dⁿ/⌡M for each i in range: ; duplicate i ! i! @0Dⁿ (-1)**i / (-1)**i/i! Σ sum * multiply sum by n! ≈ floor into int  A 12-byte version that would work if Actually had more precision: ;!╠@/1½+L@Y+  Try it online! Unlike all of the other answers (as of posting), this solution uses neither the recursive formula nor the summation formula. Instead, it uses the following formula: This formula is relatively easy to implement in Actually: !╠@/1½+L ! n! ╠ e @/ divide n! by e 1½+ add 0.5 L floor  Now, the only problem is that the formula only holds for positive n. If you try to use n = 0, the formula incorrectly yields 0. This is easily fixed, however: by applying boolean negation to the input and adding that to the output of the formula, the correct output is given for all non-negative n. Thus, the program with that correction is: ;!╠@/1½+L@Y+ ; duplicate input ! n! ╠ e @/ divide n! by e 1½+ add 0.5 L floor @Y boolean negate the other copy of the input (1 if n == 0 else 0) + add  • Keeps giving negative answers for me... – Leaky Nun Apr 24 '17 at 7:30 • @LeakyNun That's because of precision limits. For large inputs (around n = 100), (-1)**n/n! can't be represented with double-precision IEEE 754 floats. That's acceptable according to the challenge. – Mego Apr 24 '17 at 7:36 • Even for n=4... – Leaky Nun Apr 24 '17 at 7:40 • @LeakyNun Oh. I don't know why I was using floored division. Fixing it now. – Mego Apr 24 '17 at 7:42 • 16 bytes – Leaky Nun Apr 24 '17 at 8:56 Stacked, 30 bytes [:[#-::#-f\f+*]$not,\2<# !]@:f


Try it online!

Simple recursive function. Uses the fact that !n = not n for n < 2. Call as n f.

Alice, 20 18 bytes

1/o
k\i@/&wq*eqE+]


Try it online!

Explanation

This uses the recursion from Laikoni's answer, !n = n !(n-1) + (-1)n. Similar to the Funciton answer, I'm using the base case !(-1) = 1. We put that 1 on the stack with 1.. Then this...

.../o
...\i@/...


... is just the usual decimal I/O framework. The main code is actually

&wq*eqE+]k


Broken down:

&w    Push the current IP address N times to the return address stack, which
effectively begins a loop which will be executed N+1 times.
q     Push the position of the tape head, which we're abusing as the
iterator variable n.
*     Multiply !(n-1) by n.
e     Push -1.
q     Retrieve n again.
E     Raise -1 to the nth power.