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
(10212017 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 HP34C 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 41series (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.
(10212017 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.
(10212017 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.
(10212017 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. ;)
(10212017 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
(10212017 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 shiftanddivide 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 10digit 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.
RPN65 SD iOS from CuVee Soft that emulated HP65
can not do decimal factorial.
Gamo
Hello Gamo,
which formula you use for your little program?
(10222017 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
(10222017 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. Twentyfive years ago I wrote a "turbopascal" program for the gammafct 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(1x)}\] f. e.:
\[ \Gamma(3.6) =\frac{\pi}{\sin(\pi (3.6))\cdot\Gamma(4.6)}\].
(10252017 04:06 PM)peacecalc Wrote: [ > ]thank you for your answers. Twentyfive years ago I wrote a "turbopascal" program for the gammafct with real arguments.
Ah, yes, Turbo Pascal – I loved it.
(10252017 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
COS1
*
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
COS1
*
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
(10262017 06:59 AM)Dieter Wrote: [ > ] (10252017 04:06 PM)peacecalc Wrote: [ > ]thank you for your answers. Twentyfive years ago I wrote a "turbopascal" program for the gammafct 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 stackorientated and for calculating the sum with the Bernoullinumbers I used the hornerscheme.
(10262017 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
(10262017 05:34 PM)Dieter Wrote: [ > ] (10262017 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
(10272017 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 userdefined records. Pointers I could live without.