Double factorial [wp34s]
02-12-2016, 12:03 AM
Post: #1
 Gerson W. Barbosa Senior Member Posts: 1,011 Joined: Dec 2013
Double factorial [wp34s]
Code:
 // Opcode map source: wp34s.op (local directory) // Opcode SVN version: -- unknown -- // wp34s_asm.exe: Source file(s): wp34s.dat 0001 **LBL A 0002 ENTER[^] 0003 1 0004 +/- 0005 ACOS 0006 RCL[times] Y 0007 COS 0008 DEC X 0009 4 0010 / 0011 x[<->] Y 0012 2 0013 / 0014 RCL- Y 0015 RCL L 0016 x! 0017 x[<->] Y 0018 2[^x] 0019 [times] 0020 x[<->] Y 0021 # [pi] 0022 x[<->] Y 0023 y[^x] 0024 [times] 0025 RTN 0026 **LBL B 0027 RCL+ X 0028 DEC X 0029 DEC X 0030 STO+ Y 0031 x[<->] L 0032 XEQ C 0033 RTN 0034 **LBL C 0035 STO+ X 0036 RCL- Y 0037 +/- 0038 x[<->] Y 0039 XEQ 001 0040 x[<->] Y 0041 XEQ 001 0042 / 0043 RTN 0044 **LBL 001 0045 # 1/2 0046 x[<->] Y 0047 RCL[times] Y 0048 x! 0049 x[<->] Y 0050 RCL+ L 0051 2[^x] 0052 [times] 0053 # [pi] 0054 [sqrt] 0055 / 0056 RTN 0057 **LBL D 0058 ENTER[^] 0059 INC X 0060 # 002 0061 IDIV 0062 XEQ C 0063 RTN 0064 END // 64 total instructions used. // 64 total words used. // 64 single word instructions. // 0 double word instructions.

Program A: double factorial. This uses the formula used by Wolfram Alpha:

$x!!=2^{\frac{x}{2}-\frac{\cos \left ( \pi x \right )-1}{4}}\times \pi ^{\frac{\cos \left ( \pi x \right )-1}{4}}\times\left ( \frac{x}{2} \right )!$

4 A --> 8
6 A --> 48
7.77 A --> 290.3223156096666
8.88 +/- A --> 0.01106808239125682
3.33 +/- A --> -1,065802995433050
5 +/- A --> 0.3333333333333334

Program B: rising double factorial.

9.8 ENTER 3 B --> 1595.832000000000 ; 9.8 * 11.8 * 13.8 = 1595.832
8.1 ENTER 4 B --> 13957.60410000001 ; 8.1 * 10.1 * 12.1 * 14.1 = 13957.6041

Program C: falling double factorial.

7.7 ENTER 2 C --> 43.88999999999996 ; 7.7 * 5.7 = 43.89
6.6 ENTER 3 C --> 78.93600000000006 ; 6.6 * 4.6 * 2.6 = 78.936

Program D: double factorial (integer arguments). This and programs B and C uses a more simple formula found at Wikipedia.

5 +/- D --> 0.3333333333333334
3 +/- D --> -1
4 D --> 8.000000000000001
02-12-2016, 06:27 AM
Post: #2
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
Nice Functions I considered implementing natively at one point....

cos(pi x) can be replaced by (-1)^x which should save a few steps at the start and be faster and more accurate.

The "XEQ C RTN LBL C" sequence can be simplified to LBL C too I think. Likewise, the "XEQ C RTN" at the end can be replaced by GTO C

- Pauli
02-12-2016, 07:49 AM (This post was last modified: 02-12-2016 07:52 AM by Dieter.)
Post: #3
 Dieter Senior Member Posts: 2,076 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 06:27 AM)Paul Dale Wrote:  cos(pi x) can be replaced by (-1)^x which should save a few steps at the start and be faster and more accurate.

Let's not forget CNST, a handy command that IMHO is overlooked too often. CNST 86 is pi/2, and after rearranging the x!! formula to...

$$x!!=2^{\frac{x}{2}}\cdot\left(\frac{\pi}{2}\right)^{\frac{\cos \left ( \pi x \right )-1}{4}}\cdot \left(\frac{x}{2}\right)!$$

...the A program gets quite a bit shorter:

Code:
LBL A CNST 86 RCL Y (-1)^x DEC X #004 / y^x #1/2 RCLx Z 2^x STOx Y x<> L x! x RTN

Dieter
02-12-2016, 10:45 AM
Post: #4
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
CNST 86 isn't an officially supported command and cannot be entered from the keyboard...

Pauli
02-12-2016, 11:56 AM
Post: #5
 Didier Lachieze Senior Member Posts: 952 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 10:45 AM)Paul Dale Wrote:  CNST 86 isn't an officially supported command and cannot be entered from the keyboard...

It may not be officially supported but it can be entered from the keyboard, CNST is in the P.FCN menu.
02-12-2016, 02:21 PM (This post was last modified: 02-12-2016 02:24 PM by Gerson W. Barbosa.)
Post: #6
 Gerson W. Barbosa Senior Member Posts: 1,011 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 07:49 AM)Dieter Wrote:
(02-12-2016 06:27 AM)Paul Dale Wrote:  cos(pi x) can be replaced by (-1)^x which should save a few steps at the start and be faster and more accurate.

Let's not forget CNST, a handy command that IMHO is overlooked too often. CNST 86 is pi/2, and after rearranging the x!! formula to...

$$x!!=2^{\frac{x}{2}}\cdot\left(\frac{\pi}{2}\right)^{\frac{\cos \left ( \pi x \right )-1}{4}}\cdot \left(\frac{x}{2}\right)!$$

...the A program gets quite a bit shorter:

Code:
LBL A CNST 86 RCL Y (-1)^x DEC X #004 / y^x #1/2 RCLx Z 2^x STOx Y x<> L x! x RTN

7/12th shorter, actually!

The following preverves the Y register, but I need two more steps:

Code:
 LBL A ENTER ENTER || 2^x RCL L x! * x<>y (-1)^x DEC X 4 / CNST 86 x<>y y^x * RTN

Ideally, both Y and Z stack registers should be saved. Thus I would not need the alternative definition in program 01. The implementation of the rising and falling factorials might be more interesting, but they are defined in terms of the ordinary factorial and would take just a few steps each.

Walter's Blue Book is worth looking at sometimes. If I only knew the whereabouts of mine... (hidden somewhere in the house).

Thank you all for the optimization suggestions, formula rearrangement and comments.

Gerson.
02-12-2016, 07:54 PM
Post: #7
 Dieter Senior Member Posts: 2,076 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 02:21 PM)Gerson W. Barbosa Wrote:  The following preverves the Y register, but I need two more steps:

Simply replace the initial ENTER ENTER || with #1/2 RCLx Y and it's one step less again. ;-)

(02-12-2016 02:21 PM)Gerson W. Barbosa Wrote:  Ideally, both Y and Z stack registers should be saved. Thus I would not need the alternative definition in program 01.

Hmmm... this sounds quite challenging. Which does not mean it couldn't be done. ;-)

Dieter
02-12-2016, 10:33 PM (This post was last modified: 02-13-2016 10:42 PM by Paul Dale.)
Post: #8
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 11:56 AM)Didier Lachieze Wrote:  It may not be officially supported but it can be entered from the keyboard, CNST is in the P.FCN menu.

I'd forgotten about that one. A side effect for allow indirect constants -- I don't remember why we included them but do remember doing it.

The manual doesn't number the constants, instead the has user to count them. It also stops the documentation at the final user constant. [edit: this is incorrect]

I'd assumed the various system constants after that are all subject to change if XROM needed modification.

Pauli
02-12-2016, 11:52 PM
Post: #9
 Gerson W. Barbosa Senior Member Posts: 1,011 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 07:54 PM)Dieter Wrote:
(02-12-2016 02:21 PM)Gerson W. Barbosa Wrote:  The following preverves the Y register, but I need two more steps:

Simply replace the initial ENTER ENTER || with #1/2 RCLx Y and it's one step less again. ;-)

That's it! Thanks again!

(02-12-2016 07:54 PM)Dieter Wrote:
(02-12-2016 02:21 PM)Gerson W. Barbosa Wrote:  Ideally, both Y and Z stack registers should be saved. Thus I would not need the alternative definition in program 01.

Hmmm... this sounds quite challenging. Which does not mean it couldn't be done. ;-)

Why not using register I? Its content would be changed in the next complex operation anyway. Or does this look like cheating? :-)

Code:
 0001 **LBL A 0002 z[<->] I 0003 # 1/2 0004 RCL[times] Y 0005 2[^x] 0006 RCL L 0007 x! 0008 [times] 0009 x[<->] Y 0010 (-1)[^x] 0011 DEC X 0012 # 004 0013 / 0014 CNST 86 0015 x[<->] Y 0016 y[^x] 0017 [times] 0018 RCL I 0019 [<->] YZXX 0020 RTN
02-13-2016, 12:14 AM
Post: #10
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 07:54 PM)Dieter Wrote:  Hmmm... this sounds quite challenging. Which does not mean it couldn't be done. ;-)

There are instructions included in the 34S for doing this kind of thing easily:

Code:
LBL A LocR 08 STOS .00 ... X<> .00 STO L RCLS .00 RTN

Not as good as XROM's xIN/xOUT but this version will nest and you've got the input value in a register for easy access and LastX is handled correctly. Not bad for five steps.

- Pauli
02-13-2016, 12:57 AM (This post was last modified: 02-13-2016 01:02 AM by emece67.)
Post: #11
 emece67 Senior Member Posts: 363 Joined: Feb 2015
RE: Double factorial [wp34s]
This doesn't use CNST. Only one more step than Dieter's:

Code:
LBL A (-1)^x DEC X # pi # 1/2 STOx Z STOx Z STOx Y RCLx L 2^x RCL L x! x x<> Z y^x x RTN

Corrupts all the (4 level) stack, but I'm planning to put it in XROM, so that will not be a problem.

César - Information must flow.
02-13-2016, 01:36 AM (This post was last modified: 02-13-2016 01:40 AM by Gerson W. Barbosa.)
Post: #12
 Gerson W. Barbosa Senior Member Posts: 1,011 Joined: Dec 2013
RE: Double factorial [wp34s]
The following replaces my initial programs. Thank you all for your suggestions, your pointing out my obvious mistakes and for the wp34s instruction set refresh.

Code:
  0001 **LBL B 0002 RCL+ X 0003 DEC X 0004 DEC X 0005 STO+ Y 0006 x[<->] L 0007 **LBL C 0008 STO+ X 0009 RCL- Y 0010 +/- 0011 x[<->] Y 0012 XEQ A 0013 x[<->] Y 0014 XEQ A 0015 / 0016 RTN 0017 **LBL A 0018 z[<->] I 0019 # 1/2 0020 RCL[times] Y 0021 2[^x] 0022 RCL L 0023 x! 0024 [times] 0025 x[<->] Y 0026 (-1)[^x] 0027 DEC X 0028 # 004 0029 / 0030 CNST 86 0031 x[<->] Y 0032 y[^x] 0033 [times] 0034 RCL I 0035 [<->] YZXX 0036 RTN 0037 END

Program A: double factorial.

Program B: rising double factorial.

Program C: falling double factorial.

Rising factorial:

Code:
 LBL 01 x<>y Gamma X<>Y RCL+ L Gamma x<>y / RTN

Falling factorial:

Code:
 LBL 02 x<>y x! x<>y +/- RCL+ L x! / RTN
02-13-2016, 04:46 AM
Post: #13
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-13-2016 12:57 AM)emece67 Wrote:  Corrupts all the (4 level) stack, but I'm planning to put it in XROM, so that will not be a problem.

In XROM there is no reason to not use the pi/2 constant. The name is defined symbolically and if it gets renumbered, the assembler deals with it.

Pauli
02-13-2016, 06:53 AM
Post: #14
 Dieter Senior Member Posts: 2,076 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-13-2016 04:46 AM)Paul Dale Wrote:  In XROM there is no reason to not use the pi/2 constant. The name is defined symbolically and if it gets renumbered, the assembler deals with it.

Is there a reason why these useful constants were not included in the regular CONST catalog as #pi/2, #ln2, #ln10, #sqrt2pi etc.? Is it the limited space again?

Dieter
02-13-2016, 07:05 AM
Post: #15
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-13-2016 06:53 AM)Dieter Wrote:  Is there a reason why these useful constants were not included in the regular CONST catalog as #pi/2, #ln2, #ln10, #sqrt2pi etc.? Is it the limited space again?

Ask Walter Each exposed constant does consume extra bytes so space is an issue as always. I think it is twelve bytes extra (maybe a bit more) per exposed constant -- the question, as always, is when to stop. Remembering that: "pi 2 /" is three steps.

Some of the constants are also virtually useless for the user, although after digamma was dropped there aren't as many of these.

Pauli
02-13-2016, 09:15 AM
Post: #16
 emece67 Senior Member Posts: 363 Joined: Feb 2015
RE: Double factorial [wp34s]
(02-13-2016 04:46 AM)Paul Dale Wrote:  In XROM there is no reason to not use the pi/2 constant. The name is defined symbolically and if it gets renumbered, the assembler deals with it.

Thanks for the tip. I've removed some predefined constants in my build, so I try to avoid the use of CNST. "Num [pi]/2" worked fine.

César - Information must flow.
02-13-2016, 07:00 PM
Post: #17
 Dieter Senior Member Posts: 2,076 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-13-2016 07:05 AM)Paul Dale Wrote:  the question, as always, is when to stop. Remembering that: "pi 2 /" is three steps.

True – but this applies other constants as well. Some physical constants provided by the 34s are simply products or quotients of two others. Or consider cases like Φ0 and Kj – the one is simply the reciprocal of the other, so at least one of them is not required. And what about #eE? A simple 1 ex would be fine either, wouldn't it? Even without disturbing the stack like a manual [pi] 2 [/] would do.

OK, constants like ln 2 or lg e can be obtained with two program steps, and they don't even cause stack problems (well, except losing LastX). Others do, so I would trade several of the present physical constants for a simple pi/2, pi/4 or some others.

But you said that the assembler can handle several constants, I assume this applies to cases like ln 10 or sqrt(2 pi). Could you kindly provide a list of these?

Dieter
02-13-2016, 07:50 PM
Post: #18
 Marcus von Cube Senior Member Posts: 754 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-13-2016 07:00 PM)Dieter Wrote:  But you said that the assembler can handle several constants, I assume this applies to cases like ln 10 or sqrt(2 pi). Could you kindly provide a list of these?

The wp34s.op file can tell you. Here is an excerpt:
Code:
0x2000    cmd    # 1/2 0x2001    cmd    # a 0x2002    cmd    # a[sub-0] 0x2002    alias-c    # a0 0x2003    cmd    # a[sub-m] 0x2003    alias-c    # SM_luna 0x2004    cmd    # a[terra] 0x2004    alias-c    # SM_terra 0x2005    cmd    # c 0x2006    cmd    # c[sub-1] 0x2006    alias-c    # C1 0x2007    cmd    # c[sub-2] 0x2007    alias-c    # C2 0x2008    cmd    # e 0x2009    cmd    # eE 0x200a    cmd    # F 0x200b    cmd    # F[alpha] 0x200b    alias-c    # F_alpha 0x200c    cmd    # F[delta] 0x200c    alias-c    # F_delta 0x200d    cmd    # g 0x200e    cmd    # G 0x200f    cmd    # G[sub-0] 0x200f    alias-c    # Go 0x2010    cmd    # G[sub-c] 0x2010    alias-c    # catalan 0x2011    cmd    # g[sub-e] 0x2011    alias-c    # Ge 0x2012    cmd    # GM 0x2013    cmd    # h 0x2014    cmd    # [h-bar] 0x2014    alias-c    # hon2PI 0x2015    cmd    # k 0x2016    cmd    # Kj 0x2017    cmd    # l[sub-p] 0x2017    alias-c    # PlanckL 0x2018    cmd    # m[sub-e] 0x2018    alias-c    # me 0x2019    cmd    # M[sub-m] 0x2019    alias-c    # M_luna 0x201a    cmd    # m[sub-n] 0x201a    alias-c    # mn 0x201b    cmd    # m[sub-p] 0x201b    alias-c    # mp 0x201c    cmd    # M[sub-p] 0x201c    alias-c    # PlanckM 0x201d    cmd    # m[sub-u] 0x201d    alias-c    # mu 0x201e    cmd    # m[sub-u]c[^2] 0x201e    alias-c    # muc2 0x201f    cmd    # m[sub-mu] 0x201f    alias-c    # mMu 0x2020    cmd    # M[sol] 0x2020    alias-c    # M_sol 0x2021    cmd    # M[terra] 0x2021    alias-c    # M_terra 0x2022    cmd    # N[sub-A] 0x2022    alias-c    # Na 0x2023    cmd    # NaN 0x2024    cmd    # p[sub-0] 0x2024    alias-c    # atm 0x2025    cmd    # q[sub-p] 0x2025    alias-c    # PlanckQ 0x2026    cmd    # R 0x2027    cmd    # r[sub-e] 0x2027    alias-c    # Re 0x2028    cmd    # R[sub-k] 0x2028    alias-c    # Rk 0x2029    cmd    # R[sub-m] 0x2029    alias-c    # R_luna 0x202a    cmd    # R[sub-infinity] 0x202a    alias-c    # Rinf 0x202b    cmd    # R[sol] 0x202b    alias-c    # R_sol 0x202c    cmd    # R[terra] 0x202c    alias-c    # R_terra 0x202d    cmd    # Sa 0x202e    cmd    # Sb 0x202f    cmd    # Se[^2] 0x202f    alias-c    # WGS_E2 0x2030    cmd    # Se'[^2] 0x2030    alias-c    # WGS_ES2 0x2031    cmd    # Sf[^-1] 0x2031    alias-c    # WGS_F 0x2032    cmd    # T[sub-0] 0x2032    alias-c    # t 0x2033    cmd    # t[sub-p] 0x2033    alias-c    # tp 0x2034    cmd    # T[sub-p] 0x2034    alias-c    # PlanckTh 0x2035    cmd    # V[sub-m] 0x2035    alias-c    # Vm 0x2036    cmd    # Z[sub-0] 0x2036    alias-c    # Zo 0x2037    cmd    # [alpha] 0x2037    alias-c    # alpha 0x2038    cmd    # [gamma]EM 0x2038    alias-c    # EULER 0x2039    cmd    # [gamma][sub-p] 0x2039    alias-c    # gamP 0x203a    cmd    # [epsilon][sub-0] 0x203a    alias-c    # eps0 0x203b    cmd    # [lambda][sub-c] 0x203b    alias-c    # lamC 0x203c    cmd    # [lambda][sub-c][sub-n] 0x203c    alias-c    # lamCn 0x203d    cmd    # [lambda][sub-c][sub-p] 0x203d    alias-c    # lamCp 0x203e    cmd    # [mu][sub-0] 0x203e    alias-c    # mu0 0x203f    cmd    # [mu][sub-B] 0x203f    alias-c    # muB 0x2040    cmd    # [mu][sub-e] 0x2040    alias-c    # muE 0x2041    cmd    # [mu][sub-n] 0x2041    alias-c    # mun 0x2042    cmd    # [mu][sub-p] 0x2042    alias-c    # muP 0x2043    cmd    # [mu][sub-u] 0x2043    alias-c    # mu_u 0x2044    cmd    # [mu][sub-mu] 0x2044    alias-c    # mumu 0x2045    cmd    # [sigma][sub-B] 0x2045    alias-c    # sigma 0x2046    cmd    # [PHI] 0x2046    alias-c    # PHI 0x2047    cmd    # [PHI][sub-0] 0x2047    alias-c    # phi0 0x2048    cmd    # [omega] 0x2048    alias-c    # WGS_OMEGA 0x2049    cmd    # -[infinity] 0x2049    alias-c    # NEGINF 0x204a    cmd    # [infinity] 0x204a    alias-c    # INF 0x204d    cmd    # 1/eH 0x204e    cmd    # 1/eL 0x204f    cmd    # 1/[sqrt]5 0x204f    alias-c    # RECIP_SQRT5 0x2050    cmd    # 1/[sqrt][pi] 0x2050    alias-c    # RECIP_SQRTPI 0x2051    cmd    # Chi2 0x2052    cmd    # L10[^-1] 0x2052    alias-c    # RECIPLN10 0x2053    cmd    # LN2 0x2054    cmd    # LN2[^-1] 0x2054    alias-c    # RECIPLN2 0x2055    cmd    # [pi] 0x2055    alias-c    PI 0x2056    cmd    # [pi]/2 0x2056    alias-c    # PIon2 0x2057    cmd    # [sqrt]2[pi] 0x2057    alias-c    # SQRT_2_PI 0x2058    cmd    # [integral]RgB 0x2058    alias-c    # INT_R_BOUNDS

Marcus von Cube
Wehrheim, Germany
http://www.mvcsys.de
http://wp34s.sf.net
http://mvcsys.de/doc/basic-compare.html
02-13-2016, 10:41 PM
Post: #19
 Paul Dale Senior Member Posts: 1,366 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-12-2016 10:33 PM)Paul Dale Wrote:  The manual doesn't number the constants, instead the has user to count them. It also stops the documentation at the final user constant.

Walter corrected me, I was wrong here too. The manual does list the extra constants and it does number them.

Pauli
02-14-2016, 05:44 AM (This post was last modified: 02-14-2016 05:45 AM by Gerson W. Barbosa.)
Post: #20
 Gerson W. Barbosa Senior Member Posts: 1,011 Joined: Dec 2013
RE: Double factorial [wp34s]
(02-13-2016 10:41 PM)Paul Dale Wrote:
(02-12-2016 10:33 PM)Paul Dale Wrote:  The manual doesn't number the constants, instead the has user to count them. It also stops the documentation at the final user constant.

Walter corrected me, I was wrong here too. The manual does list the extra constants and it does number them.

Speaking of factorials and constants, I've found quite by chance a remarkable approximation to a Gamma-related constant (13 SD!). It's the positive real root of the polynomial equation x^28 + x^27 + x^26 + ... + x^4 + x^3 + x^2 - 130615 = 0. I think the wp34s solver can handle that, but it might be easier to check this on the hp 50g:

Code:
 « [ 1 ] { 29 } RDM 1 CON 28 [ 0 -130615 ] REPL PROOT 6 GET C→R DROP »

The four rising and falling factorials functions have been naively implemented. They might overflow when they shouldn't, like some comb and perm functions in some old calculators.

Thanks, Walter, for helping. I haven't found my copy of the Blue Book yet, but at least I've found my spare wp34s with good batteries in it :-)

Gerson.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)