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

[Matt Agajanian]

(09-06-2021 04:43 PM)Eddie W. Shore Wrote:  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

Hi. I’m entering the following code into a 34C. Since the 34C uses I as the loop register, I've changed register 01 to I. When I input any of the sample numbers in, the program stops at Step 25 of the program below with -9.99999999 99 in the display. I've checked, rechecked five times and I get the same underflow result.

Where am I going wrong?

001 - 25 13 11 LBL A
002 - 23 14 23 STO I
003 - 25 32 INT
004 - 25 1 x!
005 - 23 3 STO 3
006 - 0 0
007 - 23 2 STO 2
008 - 24 14 23 RCL I
009 - 33 EEX
010 - 3 3
011 - 71 ÷
012 - 23 14 23 STO I
013 - 25 13 0 LBL 0
014 - 24 14 23 RCL I
015 - 25 32 INT
016 - 31 ENTER
017 - 31 ENTER
018 - 1 1
019 - 32 CHS
020 - 21 x≷y
021 - 25 3 yˣ
022 - 21 x≷y
023 - 25 32 INT
024 - 25 1 x!
025 - 71 ÷
026 - 23 51 2 STO + 2
027 - 15 24 ISG
028 - 22 0 GTO 0
029 - 24 2 RCL 2
030 - 24 61 3 RCL × 3
031 - 23 3 STO 3
032 - 25 12 RTN

[Thomas Klemm]

Using this formula for \(!n\):

\(
!n = n! \sum_{k=0}^{n} \frac{(-1)^k}{k!}
\)

We can rewrite the sum using falling factorials \(n^{\underline {k}}=\prod _{j=0}^{k-1}(n-j)\):

\(
\begin{align}
!n
&= \sum_{k=0}^{n} (-1)^k \frac{n!}{k!} \\
\\
&= \sum_{k=0}^{n} (-1)^k n^{\underline {n-k}} \\
\\
&= \left( n^{\underline {n}} - \left( n^{\underline {n-1}} - \left( n^{\underline {n-2}} - \left( \cdots - \left( n^{\underline {2}} - \left( n^{\underline {1}}
- n^{\underline {0}} \right) \right) \right) \right) \right) \right) \\
\end{align}
\)

We start inside out, that is from right to left.
Therefore, we can use DSE instead of ISG.

Code:

00 { 21-Byte Prgm }
01▸LBL "!N"
02 1
03 ENTER
04▸LBL 00
05 RCL× ST Z
06 -
07 +/-
08 LASTX
09 DSE ST Z
10 GTO 00
11 X<>Y
12 END

As a bonus we get N! for free.

Example

7 XEQ "!N"

y: 5040
x: 1854

÷

2.71844660194

---

(09-11-2021 04:15 PM)C.Ret Wrote:  

From Subfactorial:

\(
!n = \left \lfloor \frac{n!+1}{e} \right \rfloor
\)

for \(n\ge1\)

Code:

00 { 14-Byte Prgm }
01▸LBL "!N"
02 N!
03 1
04 +
05 LASTX
06 E↑X
07 ÷
08 IP
09 END

[Thomas Klemm]

(11-19-2024 02:17 AM)Matt Agajanian Wrote:  Where am I going wrong?

I don't know. But you might try this program instead:

Code:

001 25,13,11 : h LBL A
002 23,14,23 : STO f I
003        1 : 1
004       31 : ENTER
005 25,13, 0 : h LBL 0
006 24,14,23 : RCL f I
007       61 : ×
008       41 : -
009       32 : CHS
010     25 0 : h LST x
011    15 23 : g DSE
012     22 0 : GTO 0
013       21 : x<>y
014    25 12 : h RTN

Example

7 A

1,854.0000

÷

2.7184

[Matt Agajanian]

(11-19-2024 06:24 AM)Thomas Klemm Wrote:  

(11-19-2024 02:17 AM)Matt Agajanian Wrote:  Where am I going wrong?

I don't know. But you might try this program instead:

Code:

001 25,13,11 : h LBL A
002 23,14,23 : STO f I
003        1 : 1
004       31 : ENTER
005 25,13, 0 : h LBL 0
006 24,14,23 : RCL f I
007       61 : ×
008       41 : -
009       32 : CHS
010     25 0 : h LST x
011    15 23 : g DSE
012     22 0 : GTO 0
013       21 : x<>y
014    25 12 : h RTN

Example

7 A

1,854.0000

÷

2.7184

Thanks. That works!

[Thomas Klemm]

(11-21-2024 01:31 AM)Matt Agajanian Wrote:  That works!

Great to hear.

This is for the HP-15C:

Code:

   001 { 42 21 11 } f LBL A
   002 {    44 25 } STO I
   003 {        1 } 1
   004 {       36 } ENTER
   005 { 42 21  0 } f LBL 0
   006 { 45 20 25 } RCL * I
   007 {       30 } -
   008 {       16 } CHS
   009 {    43 36 } g LSTx
   010 { 42  5 25 } f DSE I
   011 {    22  0 } GTO 0
   012 {       34 } X<=>Y
   013 {    43 32 } g RTN

Example

7 f A

1,854.0000

÷

2.7184