(HP65) Factorial and Gamma Function

[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