(Comp.sys.handhelds)
Item: 2193 by jurjen at cwi.nl
Author: [Jurjen NE Bos]
  Subj: HP48: Fast & Accurate Complex Factorial
  Date: Mon Feb 25 1991

- load the library in your HP48
- store it in a port by typing <port> STO (for example, 0 STO)
- turn your calculator off and on (you'll lose your stack!)
- enter 1617 ATTACH in the directory where you intend to use it

Now you have remapped the FACT function to a special version that also takes
complex arguments.

Features:
-   FACT of a name or algebraic gives ! instead.  To make a function that works
    with both algebraics and complex numbers, use the user function:
    FAC @warning: typed in by hand, download with caution
    %%HP: T(3);
      \<< \-> z
         \<<
           IF { 6 7 9 } z TYPE POS
           THEN z 1 \->LIST 'FAC' APPLY
           ELSE z FACT
           END
         \>>
       \>>
    This program should be called FAC.
    This should be downloaded or typed in AFTER you installed the library.
    Sorry for the inconvenience; as far as I know, it is impossible to
    incorporate library functions into an algebraic.
-   The library is only 495.5 bytes.
-   For real arguments, the program is identical to !.
-   This program is extremely accurate, because it uses Long Real and Long
    Complex internally.  The accuracy is more than 11 digits (most of the time
    all digits are correct), even around the poles.
-   The program is fast, because it uses internal functions only.  Of course,
    a complex factorial is hard to compute, so it still takes 540 milliseconds
    to compute (5 3) FACT.  As far as I know, the maximum computation time is
    875 milliseconds. (For real arguments, it can go up to 2180 milliseconds!)
-   The algorithm used is:
        For real numbers, ! is called.
        For numbers with real part >-6: the same algorithm as HP uses.  I will
        not tell you what it is, because it is probably patented.  I will only
        tell you that I did not see it before, and that it is fast and
        accurate.  Of course, the algorithm had to be adapted for complex
        numbers.
        For number with real part <-6, I use the formula:
                         pi * x
        FACT(X) = --------------------
                  sin(pi*x) * FACT(-X)
        The sine is compute using degrees for the real part of x to preserve
        accuracy around the poles.
-   Note that -250.000000001 FACT is slower (and less accurate) than
    (-250.000000001 0) FACT.
-   Weird errors or results can occur if an error occurs during the
    computation. I have not see this happen, and I only expect it to happen for
    extremely large arguments (absolute value > 10000).

Have fun!