HP Forums

Full Version: (42S) Subfactorial
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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.

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

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.
(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)```
(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
(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
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
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
(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
(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!
(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 :
[attachment=9806]

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 (*) ?
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
(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)
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :