+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (42S) Subfactorial (/thread-17439.html)
(42S) Subfactorial - Eddie W. Shore - 09-06-2021 04:43 PM
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 RTNVersion 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 RTNExamples
!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/swiss-micros-dm42-subfactorial.html
RE: (42S) Subfactorial - John Keith - 09-06-2021 06:50 PM
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.
RE: (42S) Subfactorial - Albert Chan - 09-07-2021 02:21 PM
(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)RE: (42S) Subfactorial - Albert Chan - 09-08-2021 10:26 PM
(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
RE: (42S) Subfactorial - Werner - 09-09-2021 07:24 AM
(09-07-2021 02:21 PM)Albert Chan Wrote:Shouldn't that beCode:
def subfac(n, r=0): # integer n > 0
for i in range(2,n+1): r = 1-i*r
return abs(r)
for i in range (2,n) ?
Cheers, Werner
RE: (42S) Subfactorial - Werner - 09-09-2021 07:45 AM
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 ENDCheers, Werner
RE: (42S) Subfactorial - Werner - 09-09-2021 12:31 PM
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 DataCheers, Werner
RE: (42S) Subfactorial - Albert Chan - 09-09-2021 03:33 PM
(09-09-2021 07:24 AM)Werner Wrote:(09-07-2021 02:21 PM)Albert Chan Wrote:Shouldn't that beCode:
def subfac(n, r=0): # integer n > 0
for i in range(2,n+1): r = 1-i*r
return abs(r)
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
RE: (42S) Subfactorial - ijabbott - 09-11-2021 08:24 AM
(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!
RE: (42S) Subfactorial on HP-15C - C.Ret - 09-11-2021 04:15 PM
(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 (*) ?
RE: (42S) Subfactorial - Gil - 09-12-2021 12:05 AM
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
RE: (42S) Subfactorial - Albert Chan - 09-12-2021 12:46 PM
(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)