Deranged !Combinatorics: Compute the Subfactorial

Question

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.)}

Your task is to compute !n, given n.

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 code-golf, 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

Answer 1

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).

Answer 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!

Answer 3

Haskell, 25 bytes

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

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

Answer 4

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.

Answer 5

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:?\≠ᵐ}ᶜ
⟦₁           Start with a list of {the input} distinct elements
  {      }ᶜ  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

Answer 6

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

Answer 7

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)ⁿ from @Laikoni's Haskell answer, with base case !0 = 1.

Try it online!

Answer 8

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

Answer 9

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.

Answer 10

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!

Answer 11

Ruby, 27 bytes

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

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

Answer 12

CJam (10 bytes)

1qi{~*)}/z

Online demo.

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

Answer 13

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
       +  # add

Answer 14

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

Answer 15

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).

Answer 16

PostScript, 81 76 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

Answer 17

PHP, 69 Bytes

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

Answer 18

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 !

Answer 19

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

Answer 20

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

Answer 21

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.

Answer 22

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'))

Answer 23

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.

Answer 24

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)ⁿ. 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.
  +     Add it to n*!(n-1).
  ]     Move the tape head to the right.
k     Jump back to the w, as long as there is still a copy of the return
      address on the return address stack. Otherwise, do nothing and exit
      the loop.