(42S) Subfactorial - Printable Version +- 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  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/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:   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 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 END``` Cheers, 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 Data``` Cheers, 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:   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 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)