Inspired by Egan Ford's work for the 12C and 12CP (and now his "new" 10C) I thought I would offer an admittedly leaner but less well appointed variation.
If you visit this page of Peter Luschny and scroll down to the section entitled "Which Approximation to Choose?", you find in blue pseudocode the following version of the Stieltjes continued fraction estimate of the Gamma function:
Sqrt[ 2*Pi/y ] Exp[ y (Log[y]  1) + 1/(12*y + 2/(5 y + 53/(42y))) ]
Now I have played around with this at some length in Mathematica and am convinced that simplifying the innermost coefficient of 53/42 to 5/4 doesn't compromise overall accuracy too much. So on making the substitution and doing a little rearranging that is friendly to the 12C/12CP, I get this.
Exp[ y (Log[y]  1) + 1/(12*y + 2/(5 (y + 1/(4y)))) ] Sqrt[ 2*Pi ] / Sqrt[y]
This gives the following 28 steps for the 12C/CP:
STO 0
4
*
1/x
RCL 0
+
5
*
2
x<>y
/
RCL 0
g 12x
+
1/x
RCL 0
g LN
1

RCL 0
*
+
g e^x
RCL 0
g sqrt
/
RCL 1
*
The setup is to first store the constant Sqrt [2 Pi] in register 1, which can be done manually or by program code if there are steps left. I like to take advantage of all 12 available digits on the 12CP so this is what I key in: 2.506628274 STO 1 63 EEX 11 CHS STO + 1.
That done, assuming the program is the top one in memory, in run mode I hit f CLEAR PRGM, enter my positive argument, hit R/S, voila. E.g., on my 12CP:
5.5 R/S gives 52.34277784 (actual is 52.34277778)
2.25 R/S gives 1.133003995 (actual 1.133003096)
.127 R/S gives 8.540153601 (actual 7.409558064)
10 R/S gives 362880.0000 (exact to 10 digits)
15.8 R/S gives 7.567994848E11 (exact to 10 digits)
I don't have a 12C "classic" to test, but I would not expect the full 10 digit accuracy for the larger input values since values on that calculator carry only 10 and not 12 digits on the stack and there would be more rounding error accumulation.
Folks interested in gamma/factorial approximations know that for most (the exceptions being things like the Lanczos and Spouge formulae) accuracy goes down for smaller values but can be improved by repeatedly applying a "shift and divide" trick based on the identity Gamma[x+1] == x Gamma[x]. Egan's routines for 12C/CP and 10C include this but as of yet mine does not. Also I do not yet handle negative arguments, since this requires a way to compute the sine function as required by reflection formula. Egan's 99 step 12C routine builds this in.
Also keep in mind that the little routine computes Gamma straight up, not x!. To get the same behaviour as f x! on the 15C, simply add 1 to the input argument before running this routine.
I hope this little offering is of interest. Keep in mind though it is very similar to Egan's work, it is not exactly the same, computing Gamma rather than x! and using a slightly different formula that, at least in my foolings, seems to be a little more accurate, especially with the user accessible 12digit precision of the 12CP. I hope this inspires some discussion by special function lovers.
Les