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.


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

31 Answers 31


Funciton, 336 bytes

Byte count assumes UTF-16 encoding with BOM.

┌─╖┌─╖  ┌─╖ 
 │┌─╖ │ ┌┐┌┘╔═╗╓┴╖
 │╘╤╝  │  │ ╚═╝╙─╜
 │┌┴╖ ┌┴╖┌┴╖ ╔═╗
 └─┘ └─┘  └───┘

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:

enter image description here

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.


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:


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.

          │ ╚═╝

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.

         ┌┴╖ ╔═╗

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


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:


Dissected version:


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!

  • \$\begingroup\$ 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 \$\endgroup\$
    – Luis Mendo
    Mar 26, 2017 at 17:18
  • \$\begingroup\$ @LuisMendo Yes I did! It's used as a command flag (here is a link to the info page). Thank you for the suggestion :). \$\endgroup\$
    – Adnan
    Mar 26, 2017 at 17:24

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.


Jelly, 7 bytes


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


Try it online!


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

⟦₁           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

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

  • 1
    \$\begingroup\$ sigh Mathematica... \$\endgroup\$
    – epicbob57
    Mar 28, 2017 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!

  • \$\begingroup\$ 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. \$\endgroup\$
    – Adnan
    Mar 26, 2017 at 17:21
  • 1
    \$\begingroup\$ Three bytes with a little reordering. ;) \$\endgroup\$
    – Dennis
    Mar 26, 2017 at 17:25
  • 1
    \$\begingroup\$ 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). \$\endgroup\$ Mar 26, 2017 at 22:07

M, 9 bytes


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


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


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

Try it online!

  • 2
    \$\begingroup\$ 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). \$\endgroup\$ Mar 26, 2017 at 18:09
  • \$\begingroup\$ 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. \$\endgroup\$
    – Dennis
    Mar 26, 2017 at 18:23
  • \$\begingroup\$ Another at the same byte count: Round[#!/E]+1-Sign@#&. Annoying initial values...! \$\endgroup\$ Mar 26, 2017 at 19:48

05AB1E, 8 bytes


Try it online!


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

Ruby, 27 bytes


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


CJam (10 bytes)


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.


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

Quipu, 73 bytes


Attempt This Online!

Based on Aaroneous Miller's factorial answer.

Equivalent pseudocode:

  a = [0, 0, 0, 0, 0]    // implicitly
  a[0] = input + 1
  goto 2
  a[1] = (a[1] + 1) * a[2] - a[3]
  a[2] = (a[2] + 1) % a[0]
  if a[2] == 0 goto 4
  a[3] = 2 - a[3]
  goto 1
  a[4] = a[1] + 1
  print a[4]

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 = 
>> ans(6)
ans =

JavaScript (ES6), 26 bytes


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, 81 76 69 bytes

Here are implementations of both formulae.


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


/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


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


% (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
  • 1
    \$\begingroup\$ -1 3 2 roll exp is a fair bit shorter than exch 2 mod 2 mul 1 sub. \$\endgroup\$ Mar 27, 2017 at 9:52
  • \$\begingroup\$ @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. \$\endgroup\$
    – PM 2Ring
    Mar 27, 2017 at 10:30
  • 1
    \$\begingroup\$ 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). \$\endgroup\$ Mar 27, 2017 at 10:38
  • \$\begingroup\$ @PeterTaylor Yes, cvi works, and it's still shorter than my previous version. \$\endgroup\$
    – PM 2Ring
    Mar 27, 2017 at 11:04

Raku, 22 bytes


Try it online!

This is an expression for the infinite sequence of subfactorials.


Python 3, 27 bytes

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

Try it online!

Based on @Laikoni's formula and the observation that n and !n always have opposite parity.

Broken down:

n odd -> n x !(n-1) also odd -> adding (-1)^n=-1 is the same as unsetting the lowest bit

n even -> n x !(n-1) also even -> adding (-1)^n=1 is the same as setting the lowest bit

therefore flipping the lowest bit works in either case


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

  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Mar 28, 2017 at 10:49

PHP, 50 44


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 !

  • \$\begingroup\$ -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.) \$\endgroup\$
    – Titus
    Mar 30, 2017 at 8:36
  • \$\begingroup\$ @Titus someone got bored ? :D would you mind giving me an example for -R and $argn? \$\endgroup\$
    – Christoph
    Mar 30, 2017 at 10:42
  • 1
    \$\begingroup\$ 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 \$\endgroup\$
    – Titus
    Mar 30, 2017 at 10:56
  • 1
    \$\begingroup\$ @Titus did I get it right ? ;-) \$\endgroup\$
    – Christoph
    Mar 30, 2017 at 11:03

Alice, 20 18 bytes


Try it online!


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


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


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.

Julia 1.0, 21 bytes


Try it online!

with the following formula from the OEIS page: $$ a(0)=1,\ a(n) = n*a(n-1) + (-1)^n $$


TI-Basic, 12 19 bytes


with the following formula from the OEIS page: $$ a(n) = n!*\sum_{k=0..n}(-1)^k/k! $$

screenshot of demo on calculator

previous answer (12 bytes)


this is invalid because 𝓮 (euler number) is stored with limited precision, which is explicitly forbiden by the challenge


05AB1E, 8 bytes


I know there is already a great 8 bytes 05AB1E answer, but this uses a completely different approach.

Try it online or verify all test cases.


  λ       # Start a recursive environment,
   è      # to output the (implicit) input'th value
0L        # Starting with a(0)=1 and a(1)=0
          # Where every following a(n) is calculated as:
    +     #  Add two values together
          #  (since + requires 2 arguments, it implicitly uses a(n-2) and a(n-1))
     N<*  #  Multiplied by n-1
          # (after which the result of a(input) is output implicitly)

Pari/GP, 26 bytes


Try it online!

Pari/GP, 28 bytes


Try it online!

Pari/GP, 15 bytes (invalid because of finite precision)


Try it online!


Mathematica, 40 bytes


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;}


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


Uses the same formula as Mego's answer.

Alternate method, 52 bytes using the PerMallows library


Actually, 18 bytes


Try it online!


;                   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:


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:

derangement formula

This formula is relatively easy to implement in Actually:

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

;             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
  • \$\begingroup\$ Keeps giving negative answers for me... \$\endgroup\$
    – Leaky Nun
    Apr 24, 2017 at 7:30
  • \$\begingroup\$ @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. \$\endgroup\$
    – user45941
    Apr 24, 2017 at 7:36
  • \$\begingroup\$ Even for n=4... \$\endgroup\$
    – Leaky Nun
    Apr 24, 2017 at 7:40
  • \$\begingroup\$ @LeakyNun Oh. I don't know why I was using floored division. Fixing it now. \$\endgroup\$
    – user45941
    Apr 24, 2017 at 7:42
  • \$\begingroup\$ 16 bytes \$\endgroup\$
    – Leaky Nun
    Apr 24, 2017 at 8:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.