(HP65) Factorial and Gamma Function

[Gamo]

Just noticed that HP 65 cannot calculate decimal Factorial and HP67 app for Android also cannot do it.

Here is a handy program using Stirling's approximation for Factorial and Gamma Function.

This approximation is good for x<70

Program:

Code:


LBL A
571
ENTER
2488320
/
CHS
STO 1
139
ENTER
51840
/
CHS
STO 2
288
1/x
STO 3
12
1/x
STO 4
CLx
RTN
LBL E
ENTER
STO 5
2
x
pi
x
Sqr
RCL 5
RCL 5
2
/
Y^x
STO 6
x
RCL 5
e^x
1/x
x
STO x6
RCL 5
1/x
ENTER
ENTER
ENTER
RCL 1
x
RCL 2
+
x
RCL 3
+
x
RCL 4
+
x
1
+
RCL 6
x
RTN

Instruction:
1. Press [A] Initialize
2. Press [E] for x!

Example: [A] Initialize then 4.25 [B]
result approximation 35.21

[Dieter]

(10-21-2017 08:32 AM)Gamo Wrote:  Just noticed that HP 65 cannot calculate decimal Factorial and HP67 app for Android also cannot do it.

What you are missing is a Gamma function. The one or other HP67 simulator app may have such a function, for instance the one from CuVee software for iPhone. Which app one are you referring to? Maybe yours has such a function as well?

BTW, the first HP with Gamma (by means of the x! key) was the HP-34C from 1979. Maybe this was even the first pocket calculator with an accurate Gamma function at all (does anyone know?). Earlier models did not feature this, and even the 41-series (introduced about the same time as the 34C) had no Gamma. Possbily to keep its factorial function compatible with the 67/97's. But a separate Gamma function (like on the 42s) would have been nice.

(10-21-2017 08:32 AM)Gamo Wrote:  Here is a handy program using Stirling's approximation for Factorial and Gamma Function.

What can you say about this approximation's accuracy? It looks good for large arguments but less so for small x, e.g. x=1 results in 0,9995. If you omit the first constant –571/2488320 the average accuracy actually seems to increase.

(10-21-2017 08:32 AM)Gamo Wrote:  This approximation is good for x<70

The approximation is good for even larger arguments, the accuracy even increases. But due to the limitation of the HP65/67's working range the max. x is near 69,9575 where the result approches 1E100, so larger x will cause an overflow error.

BTW, the ENTER after LBL E can and should be omitted and instead of [1/x] [x] you may use a simple division.

(10-21-2017 08:32 AM)Gamo Wrote:  Example: [A] 0.08 (this number always show when Initialize)

That's why I prefer a CLX or CLST at the end of such initialization routines. ;-)

(10-21-2017 08:32 AM)Gamo Wrote:  then 4.25 [B]
result approximation 35.21

Here the approximation has an absolute error of ~ –2E–5. Without the first constant it is only ~ +5E–6. ;-)

Dieter

[Dieter]

(10-21-2017 01:21 PM)Dieter Wrote:  What can you say about this approximation's accuracy? It looks good for large arguments but less so for small x, e.g. x=1 results in 0,9995. If you omit the first constant –571/2488320 the average accuracy actually seems to increase.

As already mentioned, the error is larger for small arguments and smaller for large ones. With a little bit of tweaking the coefficients this can be changed to a more evenly distributed error. And finally there is the shift-and-divide method: the approximation is only used for sufficiently large x, say x>6. For smaller x, e.g. 4.25, the approximation is calculated for 6.25 and finally the result divided by (5.25*6.25).

Here is a quick and dirty version of this idea, with modified coefficients:

Code:

LBL A
8
EEX
4
1/x
CHS
STO 1
.
0
0
2
6
9
6
CHS
STO 2
2
8
8
1/x
STO 3
1
2
1/x
STO 4
CLX
RTN
LBL E
6
STO 0
ST0/0
X<>Y
LBL 1
x>y?
GTO 2
1
+
STO*0
GTO 1
LBL 2
STO 5
2
*
PI
*
SQRT
RCL 5
RCL 5
2
/
Y^X
STO 6
*
RCL 5
e^x
/
STO*6
RCL 1
RCL 5
/
RCL 2
+
RCL 5
/
RCL 3
+
RCL 5
/
RCL 4
+
RCL 5
/
1
+
RCL 0
/
RCL 6
*
RTN

If evaluated exactly (!) the largest error should be about 1...2 units in the 9th significant digit. Due to the numeric limitations of a 10-digit calculator the error can and will be slightly higher here and there.

The result for x=4,25 now is 35,21161186. The true result is ...1185.

Dieter

[Gamo]

Dieter Thank You

Your quick modification is very good with more accurate approximation.

The HP67 APP for Android simply called in the play store as HP67 if you have Android device this app is highly recommended and free except that this particular app can not do decimal factorial.

I do have the HP67 iOS app from CuVee Soft and noticed that this version can do this no problem.

RPN-65 SD iOS from CuVee Soft that emulated HP65
can not do decimal factorial.

Gamo

[peacecalc]

Hello Gamo,

which formula you use for your little program?

[Dieter]

(10-22-2017 02:40 AM)Gamo Wrote:  Your quick modification is very good with more accurate approximation.

Here is an improved version. If uses a different technique to prevent overflow during the calculation of x^x · e^–x. In your program and my first version this is done with two consecutive multiplications of x^x/2, while now this term has been rearranged to (x/e)^x:

Code:

LBL A
8
EEX
4
1/x
CHS
STO 1
.
0
0
2
6
9
6
CHS
STO 2
2
8
8
1/x
STO 3
1
2
1/x
STO 4
CLX
RTN
LBL E
6
STO 0
STO/0
X<>Y
LBL 1
x>y?
GTO 2
1
+
STO*0
GTO 1
LBL 2
ENTER
ENTER
ENTER
1
CHS
e^x
*
X<>Y
Y^X
X<>Y
2
*
PI
*
SQRT
*
RCL 0
/
RCL 1
R↑
/
RCL 2
+
R↑
/
RCL 3
+
R↑
/
RCL 4
+
R↑
/
1
+
*
RTN

The program also no longer requires R5 and R6, most of the calculation is done on the stack.

Dieter

[Dieter]

(10-22-2017 08:49 AM)peacecalc Wrote:  which formula you use for your little program?

Since this has not been answered yet: it's Stirling's appproximation. The program uses the first terms of the series given in the section "Speed of convergence and error estimates". In my modified version the n² and n³ coefficents have been replaced by optimized values.

BTW I just realize that the linked Wikipedia article has some nice other approximations that may be worth a try, e.g. the one by Nemes (2007).

Dieter

[Gamo]

Hello peacecalc

That using Stirling's approximation This formula can be fit into limited 99 steps programmable calculator.

Gamo

[peacecalc]

Hello Dieter, hello Gamo,

thank you for your answers. Twenty-five years ago I wrote a "turbo-pascal" program for the gamma-fct with real arguments. I remember this, I also used for large arguments the stirling approx (x>10) as a example for coprozesser programming. But for smaller arguments I used the method described above (divsion by integer values). For negative number I used the formula:

\[ \Gamma(x) =\frac{\pi}{\sin(\pi x)\cdot\Gamma(1-x)}\] f. e.:

\[ \Gamma(-3.6) =\frac{\pi}{\sin(\pi (-3.6))\cdot\Gamma(4.6)}\].

[Dieter]

(10-25-2017 04:06 PM)peacecalc Wrote:  thank you for your answers. Twenty-five years ago I wrote a "turbo-pascal" program for the gamma-fct with real arguments.

Ah, yes, Turbo Pascal – I loved it.

(10-25-2017 04:06 PM)peacecalc Wrote:  I remember this, I also used for large arguments the stirling approx (x>10) as a example for coprozesser programming. But for smaller arguments I used the method described above (divsion by integer values). For negative number I used the formula: (...)

Great. Here is an HP67/97 version that applies the same formula, modified for x! instead of Gamma. Also the sin(pi*x) part is calculated in a special way to avoid roundoff errors for multiples of pi, especially if x is large.

Edit: code has been replaced with a slightly improved version

Code:

LBL e
8
EEX
4
1/x
CHS
STO 1
.
0
0
2
6
9
6
CHS
STO 2
2
8
8
1/x
STO 3
1
2
1/x
STO 4
CLX
RTN
LBL E
CF 2
1
STO 0
R↓
x≠0?
x>0?
GTO 0
SF 2
CHS
ENTER
ENTER
FRAC
1
CHS
COS-1
*
SIN
1
CHS
R↑
INT
Y^X
*
PI
X<>Y
/
STO 0
R↓
1
-
LBL 0
6
X<>Y
LBL 1
x>y?
GTO 2
1
+
STO*0
GTO 1
LBL 2
ENTER
ENTER
ENTER
1
CHS
e^x
*
X<>Y
Y^X
RCL 0
/
X<>Y
2
*
PI
*
SQRT
*
RCL 1
R↑
/
RCL 2
+
R↑
/
RCL 3
+
R↑
/
RCL 4
+
R↑
/
1
+
*
F2?
1/x
RTN

Initialize with f [e].

–3,6 [E] => –0,888685714
–4,6 [E] =>   0,246857143

Edit:
If you don't mind one more second execution time, here is a version with the constants directly in the code. Except R0 no other data registers are used, and an initialisation routine is not required either.

Code:

LBL E
CF 2
1
STO 0
R↓
x≠0?
x>0?
GTO 0
SF 2
CHS
ENTER
ENTER
FRAC
1
CHS
COS-1
*
SIN
1
CHS
RUP
INT
Y^X
*
PI
X<>Y
/
STO 0
R↓
1
-
LBL 0
6
X<>Y
LBL 1
x>y?
GTO 2
1
+
STO*0
GTO 1
LBL 2
ENTER
ENTER
ENTER
1
CHS
e^x
*
X<>Y
Y^X
RCL 0
/
X<>Y
2
*
PI
*
SQRT
*
8
EEX
4
1/x
CHS
R↑
/
.
0
0
2
6
9
6
-
R↑
/
2
8
8
1/x
+
R↑
/
1
2
1/x
+
R↑
/
1
+
*
F2?
1/x
RTN

Dieter

[Massimo Gnerucci]

(10-26-2017 06:59 AM)Dieter Wrote:  

(10-25-2017 04:06 PM)peacecalc Wrote:  thank you for your answers. Twenty-five years ago I wrote a "turbo-pascal" program for the gamma-fct with real arguments.

Ah, yes, Turbo Pascal – I loved it.

Yep! What about BCD math, for instance?

[peacecalc]

Hello friends,

a interesting remark: the coprocessor worked stack-orientated and for calculating the sum with the Bernoulli-numbers I used the horner-scheme.

[Dieter]

(10-26-2017 07:11 AM)Massimo Gnerucci Wrote:  Yep! What about BCD math, for instance?

AFAIK this was only available in version 3.0. I preferred the later versions with a decent IDE, especially from 4.0 to 6.0.

But BCD math indeed is a great plus. I wish it was available in more classic programming languages. BTW, what about the HP85/86's BASIC in this regard?

Dieter

[Massimo Gnerucci]

(10-26-2017 05:34 PM)Dieter Wrote:  

(10-26-2017 07:11 AM)Massimo Gnerucci Wrote:  Yep! What about BCD math, for instance?

AFAIK this was only available in version 3.0. I preferred the later versions with a decent IDE, especially from 4.0 to 6.0.

But BCD math indeed is a great plus. I wish it was available in more classic programming languages. BTW, what about the HP85/86's BASIC in this regard?

Dieter

Oh well, I had COMP[UTATIONAL]-3 type in COBOL! :)

[pier4r]

OT for the Pascal lovers (hi!): so do you love also the HPPL?

I mean the syntax is pretty similar.

[Gamo]

Here is the formula the Stirling series.

Gamo

[toml_12953]

(10-27-2017 08:06 AM)pier4r Wrote:  OT for the Pascal lovers (hi!): so do you love also the HPPL?

I mean the syntax is pretty similar.

Yes, it is similar and yes, I love it! I just wish they'd add a few things like enumeration and user-defined records. Pointers I could live without.