# HP Forums

Full Version: (HP65) Factorial and Gamma Function
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
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
(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
(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
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
Hello Gamo,

which formula you use for your little program?
(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 xx · e–x. In your program and my first version this is done with two consecutive multiplications of xx/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
(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
Hello peacecalc

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

Gamo
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)}$.
(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
(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?
Hello friends,

a interesting remark: the coprocessor worked stack-orientated and for calculating the sum with the Bernoulli-numbers I used the horner-scheme.
(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
(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! :)
OT for the Pascal lovers (hi!): so do you love also the HPPL?

I mean the syntax is pretty similar.
Here is the formula the Stirling series.

Gamo
(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.
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :