HP Forums

Full Version: Double factorial [wp34s]
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2
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
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 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
CNST 86 isn't an officially supported command and cannot be entered from the keyboard...

Pauli
(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 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 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 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 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-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
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.
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 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 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 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 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.
(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: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
(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-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.
Pages: 1 2
Reference URL's
• HP Forums: http://www.hpmuseum.org/forum/index.php
• :