(42S) Subfactorial

[Eddie W. Shore]

This is a request by Marko Draisma and gratitude to Mr. Draisma.

Calculating the Subfactorial

A common, and perhaps the most straight forward, formula to calculate the subfactorial is:

!n = n! × Σ((-1)^k ÷ k!, k=0 to n)

Yes, the subfactorial is written with the exclamation point first. The subfactorial finds all the possible arrangements of a set of objects where none of the objects end up in their original position.

For example, when arranging the set {1, 2, 3, 4} the subfactorial counts sets such as {2, 1, 4, 3} and {3, 4, 1, 2} but not {1, 4, 3, 2}. For the positive integers: !n < n!.

I am going to present two programs. The first will use the formula stated above.

The second uses this formula, which will not require recursion or loops:

!n = floor[ (e + 1/e) × n! ] - floor[ e × n! ]

Note: Since the N! function on the DM42 accepts only positive integers, we can use the IP (integer part) to simulate the floor function.

integer(x) = { floor(x) if x ≥ 0, ceiling(x) if x < 0

The following programs can be used on Free42, HP 42S, or Swiss Micros DM42.

Version 1: The Traditional Route

Registers used:
R01: k, counter
R02: sum register
R03: n!, later !n

Code:

01  LBL "!N"
02  STO 01
03  N!
04  STO 03
05  0
06  STO 02
07  RCL 01
08  1E3
09  ÷
10  STO 01
11  LBL 00
12  RCL 01
13  IP
14  ENTER
15  ENTER
16  -1
17  X<>Y 
18  Y↑X
19  X<>Y
20  N!
21  ÷
22  STO+ 02
23  ISG 01
24  GTO 00
25  RCL 02
26  RCL× 03
27  STO 03
28  RTN

Version 2: Closed Formula

I only put 2 in the label to distinguish the two programs.

Code:

01  LBL "!N 2"
02  N!
03  ENTER
04  ENTER
05  1
06  E↑X
07  ENTER
08  1/X
09  +
10  ×
11  IP
12  X<>Y
13  1
14  E↑X
15  ×
16  IP
17  -
18  RTN

Examples
!4 = 9
!5 = 44
!9 = 133,496

Sources:
"Calculus How To: Subfactorial" College Help Central, LLC .https://www.calculushowto.com/subfactorial/ Retrieved September 5, 2021.


Weisstein, Eric W. "Subfactorial." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/Subfactorial.html Retrieved September 5, 2021

Blog Link: http://edspi31415.blogspot.com/2021/09/s...orial.html

[John Keith]

A simpler formula for the first program is a(n) = (n-1)*(a(n-1) + a(n-2)). Maybe not as easy to implement with a 4-level stack because you have to keep the two previous values to compute the current one.

The second formula can also be simplified to round(n!/e). More information and formulas at A000166.

On Free42 or the DM42 you have over 30 digits of precision so your first program (either formula) will be exact for fairly large values of n.

[Albert Chan]

(09-06-2021 06:50 PM)John Keith Wrote:  A simpler formula for the first program is a(n) = (n-1)*(a(n-1) + a(n-2)). Maybe not as easy to implement with a 4-level stack because you have to keep the two previous values to compute the current one.

We can use https://mathworld.wolfram.com/Subfactorial.html, #5

!n = n * !(n-1) + (-1)^n

Code:

def subfac(n, r=0):     # integer n > 0
    for i in range(2,n+1): r = 1-i*r
    return abs(r)

[Albert Chan]

(09-06-2021 06:50 PM)John Keith Wrote:  The second formula can also be simplified to round(n!/e). More information and formulas at A000166.

Another perspective is with n! numerator, n!/!n is best approximation for e

With (n-1)! numerator, (n-1)!/!(n-1) is also best approximation for e

Difference of two fraction should be as tiny as possible, but not zero.

numer(n!/!n - (n-1)!/!(n-1)) = n! * !(n-1) - (n-1)! * !n = (n-1)! * (n * !(n-1) - !n)

Minimizing numerator (but not zero), last term should be ±1 (depends on parity of n)

!n = n * !(n-1) ± 1

With !1 = round(1!/e) = 0, !2 = round(2!/e) = 1, !2 = 2 * !1 + 1

!n = n * !(n-1) + (-1)^n

[Werner]

(09-07-2021 02:21 PM)Albert Chan Wrote:  

Code:

def subfac(n, r=0):     # integer n > 0
    for i in range(2,n+1): r = 1-i*r
    return abs(r)

Shouldn't that be
for i in range (2,n) ?
Cheers, Werner

[Werner]

Code:

00 { 32-Byte Prgm }
01▸LBL "!N"
02 FUNC 11
03 1ᴇ3
04 ÷
05 1
06 STO+ ST Y
07▸LBL 10
08 RCL ST Y
09 IP
10 ×
11 1
12 X<>Y
13 -
14 ISG ST Y
15 GTO 10
16 ABS
17 END

Cheers, Werner

[Werner]

If you want to make the routine behave like N!, with the Free42 extensions, include his after FUNC 11

Code:

00 REAL?
00 GTO 00
00 STR?
00 RTNERR 1 @ Alpha Data is Invalid
00 RTNERR 4 @ Invalid Type
00 LBL 00
00 ENTER
00 IP
00 X>0?
00 X#Y?
00 RTNERR 5 @ Invalid Data

Cheers, Werner

[Albert Chan]

(09-09-2021 07:24 AM)Werner Wrote:  

(09-07-2021 02:21 PM)Albert Chan Wrote:  

Code:

def subfac(n, r=0):     # integer n > 0
    for i in range(2,n+1): r = 1-i*r
    return abs(r)

Shouldn't that be
for i in range (2,n) ?

No, Python use 0-based indexing, open-intervals. To loop from 2 to n, we need range(2, n+1)

Mathematica have Subfactorial[n], but DIY is trivial, just nest (n+1) - 2 = n-1 times:

In[1]:= subfac[n_] := Nest[{1+First[#], 1-Times@@#}&, {2,0}, n-1] // Last // Abs

In[2]:= subfac[30]
Out[2]= 97581073836835777732377428235481

In[3]:= Round[30!/E]
Out[3]= 97581073836835777732377428235481

[ijabbott]

(09-06-2021 04:43 PM)Eddie W. Shore Wrote:  The subfactorial finds all the possible arrangements of a set of objects where none of the objects end up in their original position.

For example, when arranging the set {1, 2, 3, 4} the subfactorial counts sets such as {2, 1, 4, 3} and {3, 4, 1, 2} but not {1, 4, 3, 2}. For the positive integers: !n < n!.

I believe such arrangements are called "derangements" in combinatorics, so the subfactorial is the number of derangements. Great word!

[C.Ret]

(09-06-2021 06:50 PM)John Keith Wrote:  [...]The second formula can also be simplified to round(n!/e). More information and formulas at A000166.

On Free42 or the DM42 you have over 30 digits of precision so your first program (either formula) will be exact for fairly large values of n.

The round trick give me the idea for a short and fast code for my HP-15C :
   

Quote:In[2]:= subfac[30]
Out[2]= 97581073836835777732377428235481

30 f D compute !30 and f PREFIX display 97.58107385 E 30 these makes 9/10 digits are alright, fast and not bad for a truthfully Voyager oldie !

5 f D display exactly 44, let verify the statistics :
12345 - 21345 - 31245 - 41235 - 51234 *
12354 - 21354 - 31254 * 41253 * 51243 -
12435 - 21435 - 31425 - 41325 - 51324 -
12453 - 21453 * 31452 * 41352 - 51342 -
12534 - 21534 * 31524 * 41523 * 51423 *
12543 - 21543 - 31542 - 41532 * 51432 *
13245 - 23145 - 32145 - 42135 - 52134 -
13254 - 23154 * 32154 - 42153 - 52143 -
13425 - 23415 - 32415 - 42315 - 52314 -
13452 - 23451 * 32451 - 42351 - 52341 -
13524 - 23514 * 32514 - 42513 - 52413 -
13542 - 23541 - 32541 - 42531 - 52431 -
14235 - 24135 - 34125 - 43125 - 53124 *
14253 - 24153 * 34152 * 43152 * 53142 -
14325 - 24315 - 34215 - 43215 - 53214 *
14352 - 24351 - 34251 * 43251 * 53241 -
14523 - 24513 * 34512 * 43512 * 53412 *
14532 - 24531 * 34521 * 43521 * 53421 *
15234 - 25134 * 35124 * 45123 * 54123 *
15243 - 25143 - 35142 - 45132 * 54132 *
15324 - 25314 - 35214 * 45213 * 54213 *
15342 - 25341 - 35241 - 45231 * 54231 *
15423 - 25413 * 35412 * 45312 - 54312 -
15432 - 25431 * 35421 * 45321 - 54321 -

Total :
5! = 120 arrangements (- & *)
!5 = 44 derangements ( * )
How the asterisk are randomly distributed makes me believe of un vrai dérangement.

Couldn't Albert Chan please verify that with 30 elements, he exactly gets 97581073836835777732377428235481 asterisks (*) ?

[Gil]

Dividing 30! by 2.718281828... (100 digits) gives
the following digits: 97581073836835777732377428235480.9687...

The above rounded number is then correctly the one given in the previous post ending by 481.

Regards,
Gil

[Albert Chan]

(09-12-2021 12:05 AM)Gil Wrote:  Dividing 30! by 2.718281828... (100 digits) gives
the following digits: 97581073836835777732377428235480.9687...

The above rounded number is then correctly the one given in the previous post ending by 481.

Thanks, Gil

If we consider removed fractional part, (1-0.9687) = 0.0313 ≈ 0.970/31

Let P(n) = !n/n!, abs_error(P(30)) < 1/31!

Consider only n-th element, probability of derangemnt = 1-1/n
(i.e. n-th element cannot be at n-th place)

To count derangements, we have to remove over-counted cases.
Example: 12345 has all elements in the wrong spots.

Recursively apply Inclusion Exclusion Principle to remove over-counts:

P(n) = 1 - 1/1*(1 - 1/2*(1 - 1/3*( ... (1 - 1/(n-1)*(1-1/n)))))
       = 1 - 1 + 1/2! - 1/3! + 1/4! - ... + (-1)^n/n!

1/e  = 1 - 1 + 1/2! - 1/3! + 1/4! - ... + (-1)^n/n! + (-1)^(n+1)/(n+1)! + ...

With alternate signs, abs_error( P(n)-1/e ) < 1/(n+1)!

→ abs_error( !n - n!/e ) < 1/(n+1)

For n≥1, 1/(n+1) ≤ 1/2 : !n = round(n!/e)